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*/