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 */
 53:   PetscInt       M = 10;
 54:   MatNullSpace   constants,nconstants;

 56:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 57:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 58:   PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:      Create linear solver context
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   KSPCreate(PETSC_COMM_WORLD,&ksp);

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:      Create vector data structures; set function evaluation routine
 68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 70:   /*
 71:      Create distributed array (DMDA) to manage parallel grid and vectors
 72:   */
 73:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,2,NULL,&da);
 74:   DMSetFromOptions(da);
 75:   DMSetUp(da);

 77:   /*
 78:      Extract global and local vectors from DMDA; then duplicate for remaining
 79:      vectors that are the same types
 80:   */
 81:   DMCreateGlobalVector(da,&x);
 82:   VecDuplicate(x,&r);

 84:   /*
 85:      Set function evaluation routine and vector.  Whenever the nonlinear
 86:      solver needs to compute the nonlinear function, it will call this
 87:      routine.
 88:       - Note that the final routine argument is the user-defined
 89:         context that provides application-specific data for the
 90:         function evaluation routine.
 91:   */
 92:   FormRightHandSide(r,da);

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Create matrix data structure;
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 97:   DMCreateMatrix(da,&J);
 98:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
 99:   if (symmetric) {
100:     MatSetOption(J,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
101:     MatSetOption(J,MAT_SYMMETRIC,PETSC_TRUE);
102:   } else {
103:     Vec         n;
104:     PetscInt    zero = 0;
105:     PetscScalar zeros = 0.0;
106:     VecDuplicate(x,&n);
107:     /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
108:     VecSet(n,1.0);
109:     VecSetValues(n,1,&zero,&zeros,INSERT_VALUES);
110:     VecAssemblyBegin(n);
111:     VecAssemblyEnd(n);
112:     VecNormalize(n,NULL);
113:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&n,&nconstants);
114:     MatSetTransposeNullSpace(J,nconstants);
115:     MatNullSpaceDestroy(&nconstants);
116:     VecDestroy(&n);
117:   }
118:   MatSetNullSpace(J,constants);
119:   FormMatrix(J,da);

121:   KSPSetOperators(ksp,J,J);

123:   KSPSetFromOptions(ksp);
124:   KSPSolve(ksp,r,x);
125:   KSPSolveTranspose(ksp,r,x);

127:   VecDestroy(&x);
128:   VecDestroy(&r);
129:   MatDestroy(&J);
130:   MatNullSpaceDestroy(&constants);
131:   KSPDestroy(&ksp);
132:   DMDestroy(&da);
133:   PetscFinalize();
134:   return ierr;
135: }

137: /*

139:      This intentionally includes something in the right hand side that is not in the range of the matrix A.
140:      Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
141:      portion of the right hand side before solving the linear system.
142: */
143: PetscErrorCode FormRightHandSide(Vec f,void *ctx)
144: {
145:   DM             da    = (DM) ctx;
146:   PetscScalar    *ff;
148:   PetscInt       i,M,xs,xm;
149:   PetscReal      h;

152:   DMDAVecGetArray(da,f,&ff);

154:   /*
155:      Get local grid boundaries (for 1-dimensional DMDA):
156:        xs, xm  - starting grid index, width of local grid (no ghost points)
157:   */
158:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
159:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

161:   /*
162:      Compute function over locally owned part of the grid
163:      Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
164:   */
165:   h = 1.0/M;
166:   for (i=xs; i<xs+xm; i++) ff[i] = - PetscSinReal(2.0*PETSC_PI*i*h) + 1.0;

168:   /*
169:      Restore vectors
170:   */
171:   DMDAVecRestoreArray(da,f,&ff);
172:   return(0);
173: }
174: /* ------------------------------------------------------------------- */
175: PetscErrorCode FormMatrix(Mat jac,void *ctx)
176: {
177:   PetscScalar    A[3];
179:   PetscInt       i,M,xs,xm;
180:   DM             da = (DM) ctx;
181:   MatStencil     row,cols[3];
182:   PetscReal      h;

185:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

187:   /*
188:     Get range of locally owned matrix
189:   */
190:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

192:   MatZeroEntries(jac);
193:   h = 1.0/M;
194:   /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
195:   if (symmetric) {
196:     for (i=xs; i<xs+xm; i++) {
197:       row.i = i;
198:       cols[0].i = i - 1;
199:       cols[1].i = i;
200:       cols[2].i = i + 1;
201:       A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
202:       MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
203:     }
204:   } else {
205:     MatStencil  *acols;
206:     PetscScalar *avals;

208:     /* only works for one process */
209:     MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
210:     row.i = 0;
211:     PetscMalloc1(M,&acols);
212:     PetscMalloc1(M,&avals);
213:     for (i=0; i<M; i++) {
214:       acols[i].i = i;
215:       avals[i]  = (i % 2) ? 1 : -1;
216:     }
217:     MatSetValuesStencil(jac,1,&row,M,acols,avals,ADD_VALUES);
218:     PetscFree(acols);
219:     PetscFree(avals);
220:     row.i = 1;
221:     cols[0].i = - 1;
222:     cols[1].i = 1;
223:     cols[2].i = 1 + 1;
224:     A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
225:     MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
226:     for (i=2; i<xs+xm-1; i++) {
227:       row.i = i;
228:       cols[0].i = i - 1;
229:       cols[1].i = i;
230:       cols[2].i = i + 1;
231:       A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
232:       MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
233:     }
234:     row.i = M - 1 ;
235:     cols[0].i = M-2;
236:     cols[1].i = M-1;
237:     cols[2].i = M+1;
238:     A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
239:     MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
240:   }
241:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
242:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
243:   return(0);
244: }

246: /*TEST

248:    test:
249:       suffix: nonsymmetric_left
250:       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
251:       filter: sed 's/ATOL/RTOL/g'
252:       requires: !single

254:    test:
255:       suffix: nonsymmetric_right
256:       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
257:       filter: sed 's/ATOL/RTOL/g'
258:       requires: !single

260:    test:
261:       suffix: symmetric_left
262:       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
263:       requires: !single

265:    test:
266:       suffix: symmetric_right
267:       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
268:       filter: sed 's/ATOL/RTOL/g'
269:       requires: !single

271:    test:
272:       suffix: transpose_asm
273:       args: -symmetric false -ksp_monitor -ksp_view -pc_type asm -sub_pc_type lu -sub_pc_factor_zeropivot 1.e-33 -ksp_converged_reason
274:       filter: sed 's/ATOL/RTOL/g'

276: TEST*/