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;

 50:   PetscInitialize(&argc, &argv, (char *)0, help);
 51:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 52:   PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL);

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:      Create linear solver context
 56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 58:   KSPCreate(PETSC_COMM_WORLD, &ksp);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:      Create vector data structures; set function evaluation routine
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   /*
 65:      Create distributed array (DMDA) to manage parallel grid and vectors
 66:   */
 67:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, M, 1, 2, NULL, &da);
 68:   DMSetFromOptions(da);
 69:   DMSetUp(da);

 71:   /*
 72:      Extract global and local vectors from DMDA; then duplicate for remaining
 73:      vectors that are the same types
 74:   */
 75:   DMCreateGlobalVector(da, &x);
 76:   VecDuplicate(x, &r);

 78:   /*
 79:      Set function evaluation routine and vector.  Whenever the nonlinear
 80:      solver needs to compute the nonlinear function, it will call this
 81:      routine.
 82:       - Note that the final routine argument is the user-defined
 83:         context that provides application-specific data for the
 84:         function evaluation routine.
 85:   */
 86:   FormRightHandSide(r, da);

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

115:   KSPSetOperators(ksp, J, J);

117:   KSPSetFromOptions(ksp);
118:   KSPSolve(ksp, r, x);
119:   KSPSolveTranspose(ksp, r, x);

121:   VecDestroy(&x);
122:   VecDestroy(&r);
123:   MatDestroy(&J);
124:   MatNullSpaceDestroy(&constants);
125:   KSPDestroy(&ksp);
126:   DMDestroy(&da);
127:   PetscFinalize();
128:   return 0;
129: }

131: /*

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

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

147:   /*
148:      Get local grid boundaries (for 1-dimensional DMDA):
149:        xs, xm  - starting grid index, width of local grid (no ghost points)
150:   */
151:   DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
152:   DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

154:   /*
155:      Compute function over locally owned part of the grid
156:      Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
157:   */
158:   h = 1.0 / M;
159:   for (i = xs; i < xs + xm; i++) ff[i] = -PetscSinReal(2.0 * PETSC_PI * i * h) + 1.0;

161:   /*
162:      Restore vectors
163:   */
164:   DMDAVecRestoreArray(da, f, &ff);
165:   return 0;
166: }
167: /* ------------------------------------------------------------------- */
168: PetscErrorCode FormMatrix(Mat jac, void *ctx)
169: {
170:   PetscScalar A[3];
171:   PetscInt    i, M, xs, xm;
172:   DM          da = (DM)ctx;
173:   MatStencil  row, cols[3];
174:   PetscReal   h;

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

179:   /*
180:     Get range of locally owned matrix
181:   */
182:   DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

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

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

242: /*TEST

244:    test:
245:       suffix: nonsymmetric_left
246:       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
247:       filter: sed 's/ATOL/RTOL/g'
248:       requires: !single

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

256:    test:
257:       suffix: symmetric_left
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 left
259:       requires: !single

261:    test:
262:       suffix: symmetric_right
263:       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
264:       filter: sed 's/ATOL/RTOL/g'
265:       requires: !single

267:    test:
268:       suffix: transpose_asm
269:       args: -symmetric false -ksp_monitor -ksp_view -pc_type asm -sub_pc_type lu -sub_pc_factor_zeropivot 1.e-33 -ksp_converged_reason
270:       filter: sed 's/ATOL/RTOL/g'

272: TEST*/