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;
49: PetscInitialize(&argc,&argv,(char*)0,help);
50: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
51: PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create linear solver context
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: KSPCreate(PETSC_COMM_WORLD,&ksp);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Create vector data structures; set function evaluation routine
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: /*
64: Create distributed array (DMDA) to manage parallel grid and vectors
65: */
66: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,2,NULL,&da);
67: DMSetFromOptions(da);
68: DMSetUp(da);
70: /*
71: Extract global and local vectors from DMDA; then duplicate for remaining
72: vectors that are the same types
73: */
74: DMCreateGlobalVector(da,&x);
75: VecDuplicate(x,&r);
77: /*
78: Set function evaluation routine and vector. Whenever the nonlinear
79: solver needs to compute the nonlinear function, it will call this
80: routine.
81: - Note that the final routine argument is the user-defined
82: context that provides application-specific data for the
83: function evaluation routine.
84: */
85: FormRightHandSide(r,da);
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Create matrix data structure;
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: DMCreateMatrix(da,&J);
91: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
92: if (symmetric) {
93: MatSetOption(J,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
94: MatSetOption(J,MAT_SYMMETRIC,PETSC_TRUE);
95: } else {
96: Vec n;
97: PetscInt zero = 0;
98: PetscScalar zeros = 0.0;
99: VecDuplicate(x,&n);
100: /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
101: VecSet(n,1.0);
102: VecSetValues(n,1,&zero,&zeros,INSERT_VALUES);
103: VecAssemblyBegin(n);
104: VecAssemblyEnd(n);
105: VecNormalize(n,NULL);
106: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&n,&nconstants);
107: MatSetTransposeNullSpace(J,nconstants);
108: MatNullSpaceDestroy(&nconstants);
109: VecDestroy(&n);
110: }
111: MatSetNullSpace(J,constants);
112: FormMatrix(J,da);
114: KSPSetOperators(ksp,J,J);
116: KSPSetFromOptions(ksp);
117: KSPSolve(ksp,r,x);
118: KSPSolveTranspose(ksp,r,x);
120: VecDestroy(&x);
121: VecDestroy(&r);
122: MatDestroy(&J);
123: MatNullSpaceDestroy(&constants);
124: KSPDestroy(&ksp);
125: DMDestroy(&da);
126: PetscFinalize();
127: return 0;
128: }
130: /*
132: This intentionally includes something in the right hand side that is not in the range of the matrix A.
133: Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
134: portion of the right hand side before solving the linear system.
135: */
136: PetscErrorCode FormRightHandSide(Vec f,void *ctx)
137: {
138: DM da = (DM) ctx;
139: PetscScalar *ff;
140: PetscInt i,M,xs,xm;
141: PetscReal h;
144: DMDAVecGetArray(da,f,&ff);
146: /*
147: Get local grid boundaries (for 1-dimensional DMDA):
148: xs, xm - starting grid index, width of local grid (no ghost points)
149: */
150: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
151: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
153: /*
154: Compute function over locally owned part of the grid
155: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
156: */
157: h = 1.0/M;
158: for (i=xs; i<xs+xm; i++) ff[i] = - PetscSinReal(2.0*PETSC_PI*i*h) + 1.0;
160: /*
161: Restore vectors
162: */
163: DMDAVecRestoreArray(da,f,&ff);
164: return 0;
165: }
166: /* ------------------------------------------------------------------- */
167: PetscErrorCode FormMatrix(Mat jac,void *ctx)
168: {
169: PetscScalar A[3];
170: PetscInt i,M,xs,xm;
171: DM da = (DM) ctx;
172: MatStencil row,cols[3];
173: PetscReal h;
176: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
178: /*
179: Get range of locally owned matrix
180: */
181: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
183: MatZeroEntries(jac);
184: h = 1.0/M;
185: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
186: if (symmetric) {
187: for (i=xs; i<xs+xm; i++) {
188: row.i = i;
189: cols[0].i = i - 1;
190: cols[1].i = i;
191: cols[2].i = i + 1;
192: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
193: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
194: }
195: } else {
196: MatStencil *acols;
197: PetscScalar *avals;
199: /* only works for one process */
200: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
201: row.i = 0;
202: PetscMalloc1(M,&acols);
203: PetscMalloc1(M,&avals);
204: for (i=0; i<M; i++) {
205: acols[i].i = i;
206: avals[i] = (i % 2) ? 1 : -1;
207: }
208: MatSetValuesStencil(jac,1,&row,M,acols,avals,ADD_VALUES);
209: PetscFree(acols);
210: PetscFree(avals);
211: row.i = 1;
212: cols[0].i = - 1;
213: cols[1].i = 1;
214: cols[2].i = 1 + 1;
215: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
216: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
217: for (i=2; i<xs+xm-1; i++) {
218: row.i = i;
219: cols[0].i = i - 1;
220: cols[1].i = i;
221: cols[2].i = i + 1;
222: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
223: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
224: }
225: row.i = M - 1 ;
226: cols[0].i = M-2;
227: cols[1].i = M-1;
228: cols[2].i = M+1;
229: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
230: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
231: }
232: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
233: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
234: return 0;
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*/