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