Actual source code: matnull.c

  1: /*
  2:     Routines to project vectors out of null spaces.
  3: */

  5: #include <petsc/private/matimpl.h>

  7: PetscClassId MAT_NULLSPACE_CLASSID;

  9: /*@C
 10:   MatNullSpaceSetFunction - set a function that removes a null space from a vector
 11:   out of null spaces.

 13:   Logically Collective

 15:   Input Parameters:
 16: + sp  - the `MatNullSpace` null space object
 17: . rem - the function that removes the null space
 18: - ctx - context for the remove function

 20:   Level: advanced

 22: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpaceCreate()`
 23: @*/
 24: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace, Vec, void *), void *ctx)
 25: {
 26:   PetscFunctionBegin;
 28:   sp->remove = rem;
 29:   sp->rmctx  = ctx;
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: /*@C
 34:   MatNullSpaceGetVecs - get the vectors defining the null space

 36:   Not Collective

 38:   Input Parameter:
 39: . sp - null space object

 41:   Output Parameters:
 42: + has_const - `PETSC_TRUE` if the null space contains the constant vector, otherwise `PETSC_FALSE`
 43: . n         - number of vectors (excluding constant vector) in the null space
 44: - vecs      - orthonormal vectors that span the null space (excluding the constant vector), `NULL` if `n` is 0

 46:   Level: developer

 48:   Note:
 49:   These vectors and the array are owned by the `MatNullSpace` and should not be destroyed or freeded by the caller

 51: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatGetNullSpace()`, `MatGetNearNullSpace()`
 52: @*/
 53: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp, PetscBool *has_const, PetscInt *n, const Vec **vecs)
 54: {
 55:   PetscFunctionBegin;
 57:   if (has_const) *has_const = sp->has_cnst;
 58:   if (n) *n = sp->n;
 59:   if (vecs) *vecs = sp->vecs;
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: /*@
 64:   MatNullSpaceCreateRigidBody - create rigid body modes from coordinates

 66:   Collective

 68:   Input Parameter:
 69: . coords - block of coordinates of each node, must have block size set

 71:   Output Parameter:
 72: . sp - the null space

 74:   Level: advanced

 76:   Notes:
 77:   If you are solving an elasticity problem you should likely use this, in conjunction with `MatSetNearNullSpace()`, to provide information that
 78:   the `PCGAMG` preconditioner can use to construct a much more efficient preconditioner.

 80:   If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with `MatSetNullSpace()` to
 81:   provide this information to the linear solver so it can handle the null space appropriately in the linear solution.

 83: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatSetNearNullSpace()`, `MatSetNullSpace()`, `PCGAMG`
 84: @*/
 85: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords, MatNullSpace *sp)
 86: {
 87:   const PetscScalar *x;
 88:   PetscScalar       *v[6], dots[5];
 89:   Vec                vec[6];
 90:   PetscInt           n, N, dim, nmodes, i, j;
 91:   PetscReal          sN;

 93:   PetscFunctionBegin;
 94:   PetscCall(VecGetBlockSize(coords, &dim));
 95:   PetscCall(VecGetLocalSize(coords, &n));
 96:   PetscCall(VecGetSize(coords, &N));
 97:   n /= dim;
 98:   N /= dim;
 99:   sN = 1. / PetscSqrtReal((PetscReal)N);
100:   switch (dim) {
101:   case 1:
102:     PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)coords), PETSC_TRUE, 0, NULL, sp));
103:     break;
104:   case 2:
105:   case 3:
106:     nmodes = (dim == 2) ? 3 : 6;
107:     PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &vec[0]));
108:     PetscCall(VecSetSizes(vec[0], dim * n, dim * N));
109:     PetscCall(VecSetBlockSize(vec[0], dim));
110:     PetscCall(VecSetUp(vec[0]));
111:     for (i = 1; i < nmodes; i++) PetscCall(VecDuplicate(vec[0], &vec[i]));
112:     for (i = 0; i < nmodes; i++) PetscCall(VecGetArray(vec[i], &v[i]));
113:     PetscCall(VecGetArrayRead(coords, &x));
114:     for (i = 0; i < n; i++) {
115:       if (dim == 2) {
116:         v[0][i * 2 + 0] = sN;
117:         v[0][i * 2 + 1] = 0.;
118:         v[1][i * 2 + 0] = 0.;
119:         v[1][i * 2 + 1] = sN;
120:         /* Rotations */
121:         v[2][i * 2 + 0] = -x[i * 2 + 1];
122:         v[2][i * 2 + 1] = x[i * 2 + 0];
123:       } else {
124:         v[0][i * 3 + 0] = sN;
125:         v[0][i * 3 + 1] = 0.;
126:         v[0][i * 3 + 2] = 0.;
127:         v[1][i * 3 + 0] = 0.;
128:         v[1][i * 3 + 1] = sN;
129:         v[1][i * 3 + 2] = 0.;
130:         v[2][i * 3 + 0] = 0.;
131:         v[2][i * 3 + 1] = 0.;
132:         v[2][i * 3 + 2] = sN;

134:         v[3][i * 3 + 0] = x[i * 3 + 1];
135:         v[3][i * 3 + 1] = -x[i * 3 + 0];
136:         v[3][i * 3 + 2] = 0.;
137:         v[4][i * 3 + 0] = 0.;
138:         v[4][i * 3 + 1] = -x[i * 3 + 2];
139:         v[4][i * 3 + 2] = x[i * 3 + 1];
140:         v[5][i * 3 + 0] = x[i * 3 + 2];
141:         v[5][i * 3 + 1] = 0.;
142:         v[5][i * 3 + 2] = -x[i * 3 + 0];
143:       }
144:     }
145:     for (i = 0; i < nmodes; i++) PetscCall(VecRestoreArray(vec[i], &v[i]));
146:     PetscCall(VecRestoreArrayRead(coords, &x));
147:     for (i = dim; i < nmodes; i++) {
148:       /* Orthonormalize vec[i] against vec[0:i-1] */
149:       PetscCall(VecMDot(vec[i], i, vec, dots));
150:       for (j = 0; j < i; j++) dots[j] *= -1.;
151:       PetscCall(VecMAXPY(vec[i], i, dots, vec));
152:       PetscCall(VecNormalize(vec[i], NULL));
153:     }
154:     PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)coords), PETSC_FALSE, nmodes, vec, sp));
155:     for (i = 0; i < nmodes; i++) PetscCall(VecDestroy(&vec[i]));
156:   }
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: /*@C
161:   MatNullSpaceView - Visualizes a null space object.

163:   Collective; No Fortran Support

165:   Input Parameters:
166: + sp     - the null space
167: - viewer - visualization context

169:   Level: advanced

171: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `PetscViewer`, `MatNullSpaceCreate()`, `PetscViewerASCIIOpen()`
172: @*/
173: PetscErrorCode MatNullSpaceView(MatNullSpace sp, PetscViewer viewer)
174: {
175:   PetscBool iascii;

177:   PetscFunctionBegin;
179:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
181:   PetscCheckSameComm(sp, 1, viewer, 2);

183:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
184:   if (iascii) {
185:     PetscViewerFormat format;
186:     PetscInt          i;
187:     PetscCall(PetscViewerGetFormat(viewer, &format));
188:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sp, viewer));
189:     PetscCall(PetscViewerASCIIPushTab(viewer));
190:     PetscCall(PetscViewerASCIIPrintf(viewer, "Contains %" PetscInt_FMT " vector%s%s\n", sp->n, sp->n == 1 ? "" : "s", sp->has_cnst ? " and the constant" : ""));
191:     if (sp->remove) PetscCall(PetscViewerASCIIPrintf(viewer, "Has user-provided removal function\n"));
192:     if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
193:       for (i = 0; i < sp->n; i++) PetscCall(VecView(sp->vecs[i], viewer));
194:     }
195:     PetscCall(PetscViewerASCIIPopTab(viewer));
196:   }
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: /*@C
201:   MatNullSpaceCreate - Creates a `MatNullSpace` data structure used to project vectors out of null spaces.

203:   Collective

205:   Input Parameters:
206: + comm     - the MPI communicator associated with the object
207: . has_cnst - `PETSC_TRUE` if the null space contains the constant vector; otherwise `PETSC_FALSE`
208: . n        - number of vectors (excluding constant vector) in null space
209: - vecs     - the vectors that span the null space (excluding the constant vector);
210:           these vectors must be orthonormal. These vectors are NOT copied, so do not change them
211:           after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
212:           for them by one).

214:   Output Parameter:
215: . SP - the null space context

217:   Level: advanced

219:   Notes:
220:   See `MatNullSpaceSetFunction()` as an alternative way of providing the null space information instead of providing the vectors.

222:   If has_cnst is `PETSC_TRUE` you do not need to pass a constant vector in as a fourth argument to this routine, nor do you
223:   need to pass in a function that eliminates the constant function into `MatNullSpaceSetFunction()`.

225: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpaceSetFunction()`
226: @*/
227: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm, PetscBool has_cnst, PetscInt n, const Vec vecs[], MatNullSpace *SP)
228: {
229:   MatNullSpace sp;
230:   PetscInt     i;

232:   PetscFunctionBegin;
233:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", n);
234:   if (n) PetscAssertPointer(vecs, 4);
236:   PetscAssertPointer(SP, 5);
237:   if (n) {
238:     for (i = 0; i < n; i++) {
239:       /* prevent the user from changes values in the vector */
240:       PetscCall(VecLockReadPush(vecs[i]));
241:     }
242:   }
243:   if (PetscUnlikelyDebug(n)) {
244:     PetscScalar *dots;
245:     for (i = 0; i < n; i++) {
246:       PetscReal norm;
247:       PetscCall(VecNorm(vecs[i], NORM_2, &norm));
248:       PetscCheck(PetscAbsReal(norm - 1) <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)vecs[i]), PETSC_ERR_ARG_WRONG, "Vector %" PetscInt_FMT " must have 2-norm of 1.0, it is %g", i, (double)norm);
249:     }
250:     if (has_cnst) {
251:       for (i = 0; i < n; i++) {
252:         PetscScalar sum;
253:         PetscCall(VecSum(vecs[i], &sum));
254:         PetscCheck(PetscAbsScalar(sum) <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)vecs[i]), PETSC_ERR_ARG_WRONG, "Vector %" PetscInt_FMT " must be orthogonal to constant vector, inner product is %g", i, (double)PetscAbsScalar(sum));
255:       }
256:     }
257:     PetscCall(PetscMalloc1(n - 1, &dots));
258:     for (i = 0; i < n - 1; i++) {
259:       PetscInt j;
260:       PetscCall(VecMDot(vecs[i], n - i - 1, vecs + i + 1, dots));
261:       for (j = 0; j < n - i - 1; j++) {
262:         PetscCheck(PetscAbsScalar(dots[j]) <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)vecs[i]), PETSC_ERR_ARG_WRONG, "Vector %" PetscInt_FMT " must be orthogonal to vector %" PetscInt_FMT ", inner product is %g", i, i + j + 1, (double)PetscAbsScalar(dots[j]));
263:       }
264:     }
265:     PetscCall(PetscFree(dots));
266:   }

268:   *SP = NULL;
269:   PetscCall(MatInitializePackage());

271:   PetscCall(PetscHeaderCreate(sp, MAT_NULLSPACE_CLASSID, "MatNullSpace", "Null space", "Mat", comm, MatNullSpaceDestroy, MatNullSpaceView));

273:   sp->has_cnst = has_cnst;
274:   sp->n        = n;
275:   sp->vecs     = NULL;
276:   sp->alpha    = NULL;
277:   sp->remove   = NULL;
278:   sp->rmctx    = NULL;

280:   if (n) {
281:     PetscCall(PetscMalloc1(n, &sp->vecs));
282:     PetscCall(PetscMalloc1(n, &sp->alpha));
283:     for (i = 0; i < n; i++) {
284:       PetscCall(PetscObjectReference((PetscObject)vecs[i]));
285:       sp->vecs[i] = vecs[i];
286:     }
287:   }

289:   *SP = sp;
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

293: /*@
294:   MatNullSpaceDestroy - Destroys a data structure used to project vectors out of null spaces.

296:   Collective

298:   Input Parameter:
299: . sp - the null space context to be destroyed

301:   Level: advanced

303: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceRemove()`, `MatNullSpaceSetFunction()`
304: @*/
305: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
306: {
307:   PetscInt i;

309:   PetscFunctionBegin;
310:   if (!*sp) PetscFunctionReturn(PETSC_SUCCESS);
312:   if (--((PetscObject)(*sp))->refct > 0) {
313:     *sp = NULL;
314:     PetscFunctionReturn(PETSC_SUCCESS);
315:   }

317:   for (i = 0; i < (*sp)->n; i++) PetscCall(VecLockReadPop((*sp)->vecs[i]));

319:   PetscCall(VecDestroyVecs((*sp)->n, &(*sp)->vecs));
320:   PetscCall(PetscFree((*sp)->alpha));
321:   PetscCall(PetscHeaderDestroy(sp));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: /*@C
326:   MatNullSpaceRemove - Removes all the components of a null space from a vector.

328:   Collective

330:   Input Parameters:
331: + sp  - the null space context (if this is `NULL` then no null space is removed)
332: - vec - the vector from which the null space is to be removed

334:   Level: advanced

336: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
337: @*/
338: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
339: {
340:   PetscScalar sum;
341:   PetscInt    i, N;

343:   PetscFunctionBegin;
344:   if (!sp) PetscFunctionReturn(PETSC_SUCCESS);

348:   if (sp->has_cnst) {
349:     PetscCall(VecGetSize(vec, &N));
350:     if (N > 0) {
351:       PetscCall(VecSum(vec, &sum));
352:       sum = sum / ((PetscScalar)(-1.0 * N));
353:       PetscCall(VecShift(vec, sum));
354:     }
355:   }

357:   if (sp->n) {
358:     PetscCall(VecMDot(vec, sp->n, sp->vecs, sp->alpha));
359:     for (i = 0; i < sp->n; i++) sp->alpha[i] = -sp->alpha[i];
360:     PetscCall(VecMAXPY(vec, sp->n, sp->alpha, sp->vecs));
361:   }

363:   if (sp->remove) PetscCall((*sp->remove)(sp, vec, sp->rmctx));
364:   PetscFunctionReturn(PETSC_SUCCESS);
365: }

367: /*@
368:   MatNullSpaceTest  - Tests if the claimed null space is really a null space of a matrix

370:   Collective

372:   Input Parameters:
373: + sp  - the null space context
374: - mat - the matrix

376:   Output Parameter:
377: . isNull - `PETSC_TRUE` if the nullspace is valid for this matrix

379:   Level: advanced

381: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
382: @*/
383: PetscErrorCode MatNullSpaceTest(MatNullSpace sp, Mat mat, PetscBool *isNull)
384: {
385:   PetscScalar sum;
386:   PetscReal   nrm, tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
387:   PetscInt    j, n, N;
388:   Vec         l, r;
389:   PetscBool   flg1 = PETSC_FALSE, flg2 = PETSC_FALSE, consistent = PETSC_TRUE;
390:   PetscViewer viewer;

392:   PetscFunctionBegin;
395:   n = sp->n;
396:   PetscCall(PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view", &flg1, NULL));
397:   PetscCall(PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view_draw", &flg2, NULL));

399:   if (n) {
400:     PetscCall(VecDuplicate(sp->vecs[0], &l));
401:   } else {
402:     PetscCall(MatCreateVecs(mat, &l, NULL));
403:   }

405:   PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
406:   if (sp->has_cnst) {
407:     PetscCall(VecDuplicate(l, &r));
408:     PetscCall(VecGetSize(l, &N));
409:     sum = 1.0 / PetscSqrtReal(N);
410:     PetscCall(VecSet(l, sum));
411:     PetscCall(MatMult(mat, l, r));
412:     PetscCall(VecNorm(r, NORM_2, &nrm));
413:     if (nrm >= tol) consistent = PETSC_FALSE;
414:     if (flg1) {
415:       if (consistent) {
416:         PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are likely null vector"));
417:       } else {
418:         PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are unlikely null vector "));
419:       }
420:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * 1/N || = %g\n", (double)nrm));
421:     }
422:     if (!consistent && flg1) PetscCall(VecView(r, viewer));
423:     if (!consistent && flg2) PetscCall(VecView(r, viewer));
424:     PetscCall(VecDestroy(&r));
425:   }

427:   for (j = 0; j < n; j++) {
428:     PetscCall((*mat->ops->mult)(mat, sp->vecs[j], l));
429:     PetscCall(VecNorm(l, NORM_2, &nrm));
430:     if (nrm >= tol) consistent = PETSC_FALSE;
431:     if (flg1) {
432:       if (consistent) {
433:         PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " is likely null vector", j));
434:       } else {
435:         PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " unlikely null vector ", j));
436:         consistent = PETSC_FALSE;
437:       }
438:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * v[%" PetscInt_FMT "] || = %g\n", j, (double)nrm));
439:     }
440:     if (!consistent && flg1) PetscCall(VecView(l, viewer));
441:     if (!consistent && flg2) PetscCall(VecView(l, viewer));
442:   }

444:   PetscCheck(!sp->remove, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
445:   PetscCall(VecDestroy(&l));
446:   if (isNull) *isNull = consistent;
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }