Actual source code: matnull.c


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

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

  8: PetscClassId MAT_NULLSPACE_CLASSID;

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

 14:    Logically Collective

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

 21:    Level: advanced

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

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

 37:    Not Collective

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

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

 47:    Level: developer

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

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

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

 67:    Collective

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

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

 75:    Level: advanced

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

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

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

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

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

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

164:    Collective; No Fortran Support

166:    Input Parameters:
167: +  matnull - the null space
168: -  viewer - visualization context

170:    Level: advanced

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

178:   PetscFunctionBegin;
180:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
182:   PetscCheckSameComm(sp, 1, viewer, 2);

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

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

204:    Collective

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

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

218:    Level: advanced

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

223:     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
224:     need to pass in a function that eliminates the constant function into `MatNullSpaceSetFunction()`.

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

233:   PetscFunctionBegin;
234:   PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", n);
238:   if (n) {
239:     for (i = 0; i < n; i++) {
240:       /* prevent the user from changes values in the vector */
241:       PetscCall(VecLockReadPush(vecs[i]));
242:     }
243:   }
244:   if (PetscUnlikelyDebug(n)) {
245:     PetscScalar *dots;
246:     for (i = 0; i < n; i++) {
247:       PetscReal norm;
248:       PetscCall(VecNorm(vecs[i], NORM_2, &norm));
249:       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);
250:     }
251:     if (has_cnst) {
252:       for (i = 0; i < n; i++) {
253:         PetscScalar sum;
254:         PetscCall(VecSum(vecs[i], &sum));
255:         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));
256:       }
257:     }
258:     PetscCall(PetscMalloc1(n - 1, &dots));
259:     for (i = 0; i < n - 1; i++) {
260:       PetscInt j;
261:       PetscCall(VecMDot(vecs[i], n - i - 1, vecs + i + 1, dots));
262:       for (j = 0; j < n - i - 1; j++) {
263:         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]));
264:       }
265:     }
266:     PetscCall(PetscFree(dots));
267:   }

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

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

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

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

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

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

297:    Collective

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

302:    Level: advanced

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

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

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

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

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

329:    Collective

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

335:    Level: advanced

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

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

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

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

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

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

371:    Collective

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

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

380:    Level: advanced

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

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

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

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

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

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