Actual source code: ex67.c
2: static char help[] = "Krylov methods to solve u'' = f in parallel with periodic boundary conditions,\n\
3: with a singular, inconsistent system.\n\n";
5: /*
7: This tests solving singular inconsistent systems with GMRES
9: Default: Solves a symmetric system
10: -symmetric false: Solves a non-symmetric system where nullspace(A) != nullspace(A')
12: -ksp_pc_side left or right
14: See the KSPSolve() for a discussion of when right preconditioning with nullspace(A) != nullspace(A') can fail to produce the
15: norm minimizing solution.
17: Note that though this example does solve the system with right preconditioning and nullspace(A) != nullspace(A') it does not produce the
18: norm minimizing solution, that is the computed solution is not orthogonal to the nullspace(A).
20: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
21: Include "petscksp.h" so that we can use KSP solvers. Note that this
22: file automatically includes:
23: petscsys.h - base PETSc routines petscvec.h - vectors
24: petscmat.h - matrices
25: petscis.h - index sets petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers petscpc.h - preconditioners
27: petscksp.h - linear solvers
28: */
30: #include <petscdm.h>
31: #include <petscdmda.h>
32: #include <petscsnes.h>
33: #include <petsc/private/kspimpl.h>
35: PetscBool symmetric = PETSC_TRUE;
37: PetscErrorCode FormMatrix(Mat, void *);
38: PetscErrorCode FormRightHandSide(Vec, void *);
40: int main(int argc, char **argv)
41: {
42: KSP ksp;
43: Mat J;
44: DM da;
45: Vec x, r; /* vectors */
46: PetscInt M = 10;
47: MatNullSpace constants, nconstants;
50: PetscInitialize(&argc, &argv, (char *)0, help);
51: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
52: PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Create linear solver context
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: KSPCreate(PETSC_COMM_WORLD, &ksp);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create vector data structures; set function evaluation routine
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: /*
65: Create distributed array (DMDA) to manage parallel grid and vectors
66: */
67: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, M, 1, 2, NULL, &da);
68: DMSetFromOptions(da);
69: DMSetUp(da);
71: /*
72: Extract global and local vectors from DMDA; then duplicate for remaining
73: vectors that are the same types
74: */
75: DMCreateGlobalVector(da, &x);
76: VecDuplicate(x, &r);
78: /*
79: Set function evaluation routine and vector. Whenever the nonlinear
80: solver needs to compute the nonlinear function, it will call this
81: routine.
82: - Note that the final routine argument is the user-defined
83: context that provides application-specific data for the
84: function evaluation routine.
85: */
86: FormRightHandSide(r, da);
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Create matrix data structure;
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: DMCreateMatrix(da, &J);
92: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &constants);
93: if (symmetric) {
94: MatSetOption(J, MAT_SYMMETRIC, PETSC_TRUE);
95: MatSetOption(J, MAT_SYMMETRY_ETERNAL, PETSC_TRUE);
96: } else {
97: Vec n;
98: PetscInt zero = 0;
99: PetscScalar zeros = 0.0;
100: VecDuplicate(x, &n);
101: /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
102: VecSet(n, 1.0);
103: VecSetValues(n, 1, &zero, &zeros, INSERT_VALUES);
104: VecAssemblyBegin(n);
105: VecAssemblyEnd(n);
106: VecNormalize(n, NULL);
107: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 1, &n, &nconstants);
108: MatSetTransposeNullSpace(J, nconstants);
109: MatNullSpaceDestroy(&nconstants);
110: VecDestroy(&n);
111: }
112: MatSetNullSpace(J, constants);
113: FormMatrix(J, da);
115: KSPSetOperators(ksp, J, J);
117: KSPSetFromOptions(ksp);
118: KSPSolve(ksp, r, x);
119: KSPSolveTranspose(ksp, r, x);
121: VecDestroy(&x);
122: VecDestroy(&r);
123: MatDestroy(&J);
124: MatNullSpaceDestroy(&constants);
125: KSPDestroy(&ksp);
126: DMDestroy(&da);
127: PetscFinalize();
128: return 0;
129: }
131: /*
133: This intentionally includes something in the right hand side that is not in the range of the matrix A.
134: Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
135: portion of the right hand side before solving the linear system.
136: */
137: PetscErrorCode FormRightHandSide(Vec f, void *ctx)
138: {
139: DM da = (DM)ctx;
140: PetscScalar *ff;
141: PetscInt i, M, xs, xm;
142: PetscReal h;
145: DMDAVecGetArray(da, f, &ff);
147: /*
148: Get local grid boundaries (for 1-dimensional DMDA):
149: xs, xm - starting grid index, width of local grid (no ghost points)
150: */
151: DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
152: DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
154: /*
155: Compute function over locally owned part of the grid
156: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
157: */
158: h = 1.0 / M;
159: for (i = xs; i < xs + xm; i++) ff[i] = -PetscSinReal(2.0 * PETSC_PI * i * h) + 1.0;
161: /*
162: Restore vectors
163: */
164: DMDAVecRestoreArray(da, f, &ff);
165: return 0;
166: }
167: /* ------------------------------------------------------------------- */
168: PetscErrorCode FormMatrix(Mat jac, void *ctx)
169: {
170: PetscScalar A[3];
171: PetscInt i, M, xs, xm;
172: DM da = (DM)ctx;
173: MatStencil row, cols[3];
174: PetscReal h;
177: DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
179: /*
180: Get range of locally owned matrix
181: */
182: DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
184: MatZeroEntries(jac);
185: h = 1.0 / M;
186: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
187: if (symmetric) {
188: for (i = xs; i < xs + xm; i++) {
189: row.i = i;
190: cols[0].i = i - 1;
191: cols[1].i = i;
192: cols[2].i = i + 1;
193: A[0] = A[2] = 1.0 / (h * h);
194: A[1] = -2.0 / (h * h);
195: MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES);
196: }
197: } else {
198: MatStencil *acols;
199: PetscScalar *avals;
201: /* only works for one process */
202: MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
203: row.i = 0;
204: PetscMalloc1(M, &acols);
205: PetscMalloc1(M, &avals);
206: for (i = 0; i < M; i++) {
207: acols[i].i = i;
208: avals[i] = (i % 2) ? 1 : -1;
209: }
210: MatSetValuesStencil(jac, 1, &row, M, acols, avals, ADD_VALUES);
211: PetscFree(acols);
212: PetscFree(avals);
213: row.i = 1;
214: cols[0].i = -1;
215: cols[1].i = 1;
216: cols[2].i = 1 + 1;
217: A[0] = A[2] = 1.0 / (h * h);
218: A[1] = -2.0 / (h * h);
219: MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES);
220: for (i = 2; i < xs + xm - 1; i++) {
221: row.i = i;
222: cols[0].i = i - 1;
223: cols[1].i = i;
224: cols[2].i = i + 1;
225: A[0] = A[2] = 1.0 / (h * h);
226: A[1] = -2.0 / (h * h);
227: MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES);
228: }
229: row.i = M - 1;
230: cols[0].i = M - 2;
231: cols[1].i = M - 1;
232: cols[2].i = M + 1;
233: A[0] = A[2] = 1.0 / (h * h);
234: A[1] = -2.0 / (h * h);
235: MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES);
236: }
237: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
238: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
239: return 0;
240: }
242: /*TEST
244: test:
245: suffix: nonsymmetric_left
246: args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
247: filter: sed 's/ATOL/RTOL/g'
248: requires: !single
250: test:
251: suffix: nonsymmetric_right
252: args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
253: filter: sed 's/ATOL/RTOL/g'
254: requires: !single
256: test:
257: suffix: symmetric_left
258: args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
259: requires: !single
261: test:
262: suffix: symmetric_right
263: args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
264: filter: sed 's/ATOL/RTOL/g'
265: requires: !single
267: test:
268: suffix: transpose_asm
269: args: -symmetric false -ksp_monitor -ksp_view -pc_type asm -sub_pc_type lu -sub_pc_factor_zeropivot 1.e-33 -ksp_converged_reason
270: filter: sed 's/ATOL/RTOL/g'
272: TEST*/