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: #if defined(PETSC_USE_LOG)
37: PetscLogStage stage;
38: #endif
41: PetscInitialize(&argc, &args, (char *)0, help);
42: n1 = 64;
43: n2 = 64;
44: PetscOptionsGetInt(NULL, NULL, "-n1", &n1, NULL);
45: PetscOptionsGetInt(NULL, NULL, "-n2", &n2, NULL);
47: h = 1.0 / n1;
48: gamma = 4.0;
49: beta = 0.01;
50: PetscOptionsGetReal(NULL, NULL, "-h", &h, NULL);
51: PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL);
52: PetscOptionsGetReal(NULL, NULL, "-beta", &beta, NULL);
53: gamma = gamma / h;
54: beta = beta / (h * h);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Compute the matrix and set right-hand-side vector.
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: /*
60: Create parallel matrix, specifying only its global dimensions.
61: When using MatCreate(), the matrix format can be specified at
62: runtime. Also, the parallel partitioning of the matrix is
63: determined by PETSc at runtime.
65: Performance tuning note: For problems of substantial size,
66: preallocation of matrix memory is crucial for attaining good
67: performance. See the matrix chapter of the users manual for details.
68: */
69: MatCreate(PETSC_COMM_WORLD, &A);
70: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n1 * n2, n1 * n2);
71: MatSetFromOptions(A);
72: MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
73: MatSeqAIJSetPreallocation(A, 5, NULL);
74: MatSetUp(A);
76: /*
77: Currently, all PETSc parallel matrix formats are partitioned by
78: contiguous chunks of rows across the processors. Determine which
79: rows of the matrix are locally owned.
80: */
81: MatGetOwnershipRange(A, &Istart, &Iend);
83: /*
84: Set matrix elements for the 2-D, five-point stencil in parallel.
85: - Each processor needs to insert only elements that it owns
86: locally (but any non-local elements will be sent to the
87: appropriate processor during matrix assembly).
88: - Always specify global rows and columns of matrix entries.
89: */
90: PetscLogStageRegister("Assembly", &stage);
91: PetscLogStagePush(stage);
92: co1 = gamma * h * h / 2.0;
93: co2 = beta * h * h;
94: for (Ii = Istart; Ii < Iend; Ii++) {
95: i = Ii / n2;
96: j = Ii - i * n2;
97: if (i > 0) {
98: J = Ii - n2;
99: v = -1.0 + co1 * (PetscScalar)i;
100: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
101: }
102: if (i < n1 - 1) {
103: J = Ii + n2;
104: v = -1.0 + co1 * (PetscScalar)i;
105: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
106: }
107: if (j > 0) {
108: J = Ii - 1;
109: v = -1.0 + co1 * (PetscScalar)j;
110: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
111: }
112: if (j < n2 - 1) {
113: J = Ii + 1;
114: v = -1.0 + co1 * (PetscScalar)j;
115: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
116: }
117: v = 4.0 + co2;
118: MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
119: }
121: /*
122: Assemble matrix, using the 2-step process:
123: MatAssemblyBegin(), MatAssemblyEnd()
124: Computations can be done while messages are in transition
125: by placing code between these two statements.
126: */
127: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
129: PetscLogStagePop();
131: /*
132: Create parallel vectors.
133: - We form 1 vector from scratch and then duplicate as needed.
134: - When using VecCreate(), VecSetSizes and VecSetFromOptions()
135: in this example, we specify only the
136: vector's global dimension; the parallel partitioning is determined
137: at runtime.
138: - When solving a linear system, the vectors and matrices MUST
139: be partitioned accordingly. PETSc automatically generates
140: appropriately partitioned matrices and vectors when MatCreate()
141: and VecCreate() are used with the same communicator.
142: - The user can alternatively specify the local vector and matrix
143: dimensions when more sophisticated partitioning is needed
144: (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
145: below).
146: */
147: VecCreate(PETSC_COMM_WORLD, &b);
148: VecSetSizes(b, PETSC_DECIDE, n1 * n2);
149: VecSetFromOptions(b);
150: VecDuplicate(b, &x);
151: VecDuplicate(b, &u);
153: /*
154: Set right-hand side.
155: */
156: VecSet(b, 1.0);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Create the linear solver and set various options
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: /*
162: Create linear solver context
163: */
164: KSPCreate(PETSC_COMM_WORLD, &ksp);
166: /*
167: Set operators. Here the matrix that defines the linear system
168: also serves as the preconditioning matrix.
169: */
170: KSPSetOperators(ksp, A, A);
172: /*
173: Set linear solver defaults for this problem (optional).
174: - By extracting the KSP and PC contexts from the KSP context,
175: we can then directly call any KSP and PC routines to set
176: various options.
177: */
178: KSPSetTolerances(ksp, 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, 200);
180: /*
181: Set runtime options, e.g.,
182: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
183: These options will override those specified above as long as
184: KSPSetFromOptions() is called _after_ any other customization
185: routines.
186: */
187: KSPSetFromOptions(ksp);
189: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: Solve the linear system
191: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193: KSPSolve(ksp, b, x);
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Clean up
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: /*
199: Free work space. All PETSc objects should be destroyed when they
200: are no longer needed.
201: */
202: KSPDestroy(&ksp);
203: VecDestroy(&u);
204: VecDestroy(&x);
205: VecDestroy(&b);
206: MatDestroy(&A);
208: /*
209: Always call PetscFinalize() before exiting a program. This routine
210: - finalizes the PETSc libraries as well as MPI
211: - provides summary and diagnostic information if certain runtime
212: options are chosen (e.g., -log_view).
213: */
214: PetscFinalize();
215: return 0;
216: }
218: /*TEST
220: test:
221: nsize: 8
222: 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
224: test:
225: suffix: 2
226: nsize: 8
227: 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
228: output_file: output/ex38_1.out
230: TEST*/