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