Actual source code: ex86.c
1: static char help[] = "Solves a one-dimensional steady upwind advection system with KSP.\n\n";
3: /*
4: Include "petscksp.h" so that we can use KSP solvers. Note that this file
5: automatically includes:
6: petscsys.h - base PETSc routines petscvec.h - vectors
7: petscmat.h - matrices petscpc.h - preconditioners
8: petscis.h - index sets
9: petscviewer.h - viewers
10: */
11: #include <petscksp.h>
13: int main(int argc, char **args)
14: {
15: Vec x, b, work_vec; /* approx solution, RHS, work vector */
16: Mat A; /* linear system matrix */
17: KSP ksp; /* linear solver context */
18: PC pc; /* preconditioner context */
19: PetscInt i, j, n = 10, col[2];
20: PetscScalar work_scalar, value[2];
21: PetscRandom r;
23: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
26: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
28: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: Compute the matrix and right-hand-side vector that define
30: the linear system, Ax = b.
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: /*
34: Create vectors. Note that we form 1 vector from scratch and
35: then duplicate as needed.
36: */
37: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
38: PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
39: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
40: PetscCall(VecSetFromOptions(x));
41: PetscCall(VecDuplicate(x, &b));
42: PetscCall(VecDuplicate(x, &work_vec));
44: /*
45: Create matrix. When using MatCreate(), the matrix format can
46: be specified at runtime.
48: Performance tuning note: For problems of substantial size,
49: preallocation of matrix memory is crucial for attaining good
50: performance. See the matrix chapter of the users manual for details.
51: */
52: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
53: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
54: PetscCall(MatSetFromOptions(A));
55: PetscCall(MatSetUp(A));
57: /*
58: Assemble matrix
59: */
60: value[0] = -1.0;
61: value[1] = 1.0;
62: for (i = 1; i < n; i++) {
63: col[0] = i - 1;
64: col[1] = i;
65: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
66: }
67: i = 0;
68: j = 0;
69: work_scalar = 1;
70: PetscCall(MatSetValues(A, 1, &i, 1, &j, &work_scalar, INSERT_VALUES));
71: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
74: /*
75: Create x, b
76: */
77: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
78: PetscCall(VecSetRandom(x, r));
79: PetscCall(VecSetRandom(work_vec, r));
80: PetscCall(MatMult(A, work_vec, b));
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Create the linear solver and set various options
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
87: /*
88: Set operators. Here the matrix that defines the linear system
89: also serves as the matrix that defines the preconditioner.
90: */
91: PetscCall(KSPSetOperators(ksp, A, A));
93: /*
94: Set linear solver defaults for this problem (optional).
95: - By extracting the KSP and PC contexts from the KSP context,
96: we can then directly call any KSP and PC routines to set
97: various options.
98: - The following four statements are optional; all of these
99: parameters could alternatively be specified at runtime via
100: KSPSetFromOptions();
101: */
102: PetscCall(KSPGetPC(ksp, &pc));
103: PetscCall(PCSetType(pc, PCBJACOBI));
104: PetscCall(KSPSetTolerances(ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
106: /*
107: Set runtime options, e.g.,
108: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
109: These options will override those specified above as long as
110: KSPSetFromOptions() is called _after_ any other customization
111: routines.
112: */
113: PetscCall(KSPSetFromOptions(ksp));
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Solve the linear system
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: PetscCall(KSPSolve(ksp, b, x));
120: /*
121: Free work space. All PETSc objects should be destroyed when they
122: are no longer needed.
123: */
124: PetscCall(VecDestroy(&x));
125: PetscCall(VecDestroy(&b));
126: PetscCall(VecDestroy(&work_vec));
127: PetscCall(MatDestroy(&A));
128: PetscCall(KSPDestroy(&ksp));
129: PetscCall(PetscRandomDestroy(&r));
131: /*
132: Always call PetscFinalize() before exiting a program. This routine
133: - finalizes the PETSc libraries as well as MPI
134: - provides summary and diagnostic information if certain runtime
135: options are chosen (e.g., -log_view).
136: */
137: PetscCall(PetscFinalize());
138: return 0;
139: }
141: /*TEST
143: testset:
144: requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
145: args: -ksp_monitor -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_grid_sweeps_down 0 -pc_hypre_boomeramg_grid_sweeps_up 1 -pc_hypre_boomeramg_grid_sweeps_coarse 2 -pc_hypre_boomeramg_max_levels 2 -ksp_rtol 1e-7 -pc_hypre_boomeramg_max_coarse_size 16 -n 33 -ksp_max_it 30 -pc_hypre_boomeramg_relax_type_all Jacobi
146: test:
147: suffix: hypre
148: test:
149: suffix: hypre_air
150: args: -pc_hypre_boomeramg_restriction_type 1 -pc_hypre_boomeramg_postrelax F -pc_hypre_boomeramg_grid_sweeps_down 1 -pc_hypre_boomeramg_prerelax C
152: TEST*/