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