Actual source code: ex38.c
1: /*
3: mpiexec -n 8 ./ex38 -ksp_type fbcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -ksp_monitor -n1 64 -n2 64
5: Contributed by Jie Chen for testing flexible BiCGStab algorithm
6: */
8: static char help[] = "Solves the PDE (in 2D) -laplacian(u) + gamma x dot grad(u) + beta u = 1\n\
9: with zero Dirichlet condition. The discretization is standard centered\n\
10: difference. Input parameters include:\n\
11: -n1 : number of mesh points in 1st dimension (default 64)\n\
12: -n2 : number of mesh points in 2nd dimension (default 64)\n\
13: -h : spacing between mesh points (default 1/n1)\n\
14: -gamma : gamma (default 4/h)\n\
15: -beta : beta (default 0.01/h^2)\n\n";
17: /*
18: Include "petscksp.h" so that we can use KSP solvers. Note that this file
19: automatically includes:
20: petscsys.h - base PETSc routines petscvec.h - vectors
21: petscmat.h - matrices
22: petscis.h - index sets petscksp.h - Krylov subspace methods
23: petscviewer.h - viewers petscpc.h - preconditioners
24: */
25: #include <petscksp.h>
27: int main(int argc, char **args)
28: {
29: Vec x, b, u; /* approx solution, RHS, working vector */
30: Mat A; /* linear system matrix */
31: KSP ksp; /* linear solver context */
32: PetscInt n1, n2; /* parameters */
33: PetscReal h, gamma, beta; /* parameters */
34: PetscInt i, j, Ii, J, Istart, Iend;
35: PetscScalar v, co1, co2;
36: PetscLogStage stage;
38: PetscFunctionBeginUser;
39: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
40: n1 = 64;
41: n2 = 64;
42: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n1", &n1, NULL));
43: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n2", &n2, NULL));
45: h = 1.0 / n1;
46: gamma = 4.0;
47: beta = 0.01;
48: PetscCall(PetscOptionsGetReal(NULL, NULL, "-h", &h, NULL));
49: PetscCall(PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL));
50: PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &beta, NULL));
51: gamma = gamma / h;
52: beta = beta / (h * h);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Compute the matrix and set right-hand-side vector.
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: /*
58: Create parallel matrix, specifying only its global dimensions.
59: When using MatCreate(), the matrix format can be specified at
60: runtime. Also, the parallel partitioning of the matrix is
61: determined by PETSc at runtime.
63: Performance tuning note: For problems of substantial size,
64: preallocation of matrix memory is crucial for attaining good
65: performance. See the matrix chapter of the users manual for details.
66: */
67: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
68: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n1 * n2, n1 * n2));
69: PetscCall(MatSetFromOptions(A));
70: PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
71: PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
72: PetscCall(MatSetUp(A));
74: /*
75: Currently, all PETSc parallel matrix formats are partitioned by
76: contiguous chunks of rows across the processors. Determine which
77: rows of the matrix are locally owned.
78: */
79: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
81: /*
82: Set matrix elements for the 2-D, five-point stencil in parallel.
83: - Each processor needs to insert only elements that it owns
84: locally (but any non-local elements will be sent to the
85: appropriate processor during matrix assembly).
86: - Always specify global rows and columns of matrix entries.
87: */
88: PetscCall(PetscLogStageRegister("Assembly", &stage));
89: PetscCall(PetscLogStagePush(stage));
90: co1 = gamma * h * h / 2.0;
91: co2 = beta * h * h;
92: for (Ii = Istart; Ii < Iend; Ii++) {
93: i = Ii / n2;
94: j = Ii - i * n2;
95: if (i > 0) {
96: J = Ii - n2;
97: v = -1.0 + co1 * (PetscScalar)i;
98: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
99: }
100: if (i < n1 - 1) {
101: J = Ii + n2;
102: v = -1.0 + co1 * (PetscScalar)i;
103: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
104: }
105: if (j > 0) {
106: J = Ii - 1;
107: v = -1.0 + co1 * (PetscScalar)j;
108: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
109: }
110: if (j < n2 - 1) {
111: J = Ii + 1;
112: v = -1.0 + co1 * (PetscScalar)j;
113: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
114: }
115: v = 4.0 + co2;
116: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
117: }
119: /*
120: Assemble matrix, using the 2-step process:
121: MatAssemblyBegin(), MatAssemblyEnd()
122: Computations can be done while messages are in transition
123: by placing code between these two statements.
124: */
125: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
126: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
127: PetscCall(PetscLogStagePop());
129: /*
130: Create parallel vectors.
131: - We form 1 vector from scratch and then duplicate as needed.
132: - When using VecCreate(), VecSetSizes and VecSetFromOptions()
133: in this example, we specify only the
134: vector's global dimension; the parallel partitioning is determined
135: at runtime.
136: - When solving a linear system, the vectors and matrices MUST
137: be partitioned accordingly. PETSc automatically generates
138: appropriately partitioned matrices and vectors when MatCreate()
139: and VecCreate() are used with the same communicator.
140: - The user can alternatively specify the local vector and matrix
141: dimensions when more sophisticated partitioning is needed
142: (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
143: below).
144: */
145: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
146: PetscCall(VecSetSizes(b, PETSC_DECIDE, n1 * n2));
147: PetscCall(VecSetFromOptions(b));
148: PetscCall(VecDuplicate(b, &x));
149: PetscCall(VecDuplicate(b, &u));
151: /*
152: Set right-hand side.
153: */
154: PetscCall(VecSet(b, 1.0));
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Create the linear solver and set various options
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: /*
160: Create linear solver context
161: */
162: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
164: /*
165: Set operators. Here the matrix that defines the linear system
166: also serves as the preconditioning matrix.
167: */
168: PetscCall(KSPSetOperators(ksp, A, A));
170: /*
171: Set linear solver defaults for this problem (optional).
172: - By extracting the KSP and PC contexts from the KSP context,
173: we can then directly call any KSP and PC routines to set
174: various options.
175: */
176: PetscCall(KSPSetTolerances(ksp, 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, 200));
178: /*
179: Set runtime options, e.g.,
180: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
181: These options will override those specified above as long as
182: KSPSetFromOptions() is called _after_ any other customization
183: routines.
184: */
185: PetscCall(KSPSetFromOptions(ksp));
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Solve the linear system
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: PetscCall(KSPSolve(ksp, b, x));
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Clean up
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: /*
197: Free work space. All PETSc objects should be destroyed when they
198: are no longer needed.
199: */
200: PetscCall(KSPDestroy(&ksp));
201: PetscCall(VecDestroy(&u));
202: PetscCall(VecDestroy(&x));
203: PetscCall(VecDestroy(&b));
204: PetscCall(MatDestroy(&A));
206: /*
207: Always call PetscFinalize() before exiting a program. This routine
208: - finalizes the PETSc libraries as well as MPI
209: - provides summary and diagnostic information if certain runtime
210: options are chosen (e.g., -log_view).
211: */
212: PetscCall(PetscFinalize());
213: return 0;
214: }
216: /*TEST
218: test:
219: nsize: 8
220: args: -ksp_type fbcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 64 -n2 64
222: test:
223: suffix: 2
224: nsize: 8
225: args: -ksp_type qmrcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 64 -n2 64
226: output_file: output/ex38_1.out
228: TEST*/