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: /*@
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()`, `MatNullSpaceRemoveFn`
23: @*/
24: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, MatNullSpaceRemoveFn *rem, PetscCtx ctx)
25: {
26: PetscFunctionBegin;
28: sp->remove = rem;
29: sp->rmctx = ctx;
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: /*@
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 - returns array of length `n` containing the 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 returned are owned by the `MatNullSpace` and should not be destroyed or freeded by the caller
51: Fortran Note:
52: Call `MatNullSpaceRestoreVecs()` when the array of `Vec` is no longer needed
54: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatGetNullSpace()`, `MatGetNearNullSpace()`
55: @*/
56: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp, PetscBool *has_const, PetscInt *n, const Vec *vecs[])
57: {
58: PetscFunctionBegin;
60: if (has_const) *has_const = sp->has_cnst;
61: if (n) *n = sp->n;
62: if (vecs) *vecs = sp->vecs;
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: /*@
67: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
69: Collective
71: Input Parameter:
72: . coords - block of coordinates of each node, must have block size set
74: Output Parameter:
75: . sp - the null space
77: Level: advanced
79: Notes:
80: If you are solving an elasticity problem you should likely use this, in conjunction with `MatSetNearNullSpace()`, to provide information that
81: the `PCGAMG` preconditioner can use to construct a much more efficient preconditioner.
83: If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with `MatSetNullSpace()` to
84: provide this information to the linear solver so it can handle the null space appropriately in the linear solution.
86: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatSetNearNullSpace()`, `MatSetNullSpace()`, `PCGAMG`
87: @*/
88: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords, MatNullSpace *sp)
89: {
90: const PetscScalar *x;
91: PetscScalar *v[6], dots[5];
92: Vec vec[6];
93: PetscInt n, N, dim, nmodes, i, j;
94: PetscReal sN;
96: PetscFunctionBegin;
97: PetscCall(VecGetBlockSize(coords, &dim));
98: PetscCall(VecGetLocalSize(coords, &n));
99: PetscCall(VecGetSize(coords, &N));
100: n /= dim;
101: N /= dim;
102: sN = 1. / PetscSqrtReal((PetscReal)N);
103: switch (dim) {
104: case 1:
105: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)coords), PETSC_TRUE, 0, NULL, sp));
106: break;
107: case 2:
108: case 3:
109: nmodes = (dim == 2) ? 3 : 6;
110: PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &vec[0]));
111: PetscCall(VecSetSizes(vec[0], dim * n, dim * N));
112: PetscCall(VecSetBlockSize(vec[0], dim));
113: PetscCall(VecSetUp(vec[0]));
114: for (i = 1; i < nmodes; i++) PetscCall(VecDuplicate(vec[0], &vec[i]));
115: for (i = 0; i < nmodes; i++) PetscCall(VecGetArray(vec[i], &v[i]));
116: PetscCall(VecGetArrayRead(coords, &x));
117: for (i = 0; i < n; i++) {
118: if (dim == 2) {
119: v[0][i * 2 + 0] = sN;
120: v[0][i * 2 + 1] = 0.;
121: v[1][i * 2 + 0] = 0.;
122: v[1][i * 2 + 1] = sN;
123: /* Rotations */
124: v[2][i * 2 + 0] = -x[i * 2 + 1];
125: v[2][i * 2 + 1] = x[i * 2 + 0];
126: } else {
127: v[0][i * 3 + 0] = sN;
128: v[0][i * 3 + 1] = 0.;
129: v[0][i * 3 + 2] = 0.;
130: v[1][i * 3 + 0] = 0.;
131: v[1][i * 3 + 1] = sN;
132: v[1][i * 3 + 2] = 0.;
133: v[2][i * 3 + 0] = 0.;
134: v[2][i * 3 + 1] = 0.;
135: v[2][i * 3 + 2] = sN;
137: v[3][i * 3 + 0] = x[i * 3 + 1];
138: v[3][i * 3 + 1] = -x[i * 3 + 0];
139: v[3][i * 3 + 2] = 0.;
140: v[4][i * 3 + 0] = 0.;
141: v[4][i * 3 + 1] = -x[i * 3 + 2];
142: v[4][i * 3 + 2] = x[i * 3 + 1];
143: v[5][i * 3 + 0] = x[i * 3 + 2];
144: v[5][i * 3 + 1] = 0.;
145: v[5][i * 3 + 2] = -x[i * 3 + 0];
146: }
147: }
148: for (i = 0; i < nmodes; i++) PetscCall(VecRestoreArray(vec[i], &v[i]));
149: PetscCall(VecRestoreArrayRead(coords, &x));
150: for (i = dim; i < nmodes; i++) {
151: /* Orthonormalize vec[i] against vec[0:i-1] */
152: PetscCall(VecMDot(vec[i], i, vec, dots));
153: for (j = 0; j < i; j++) dots[j] *= -1.;
154: PetscCall(VecMAXPY(vec[i], i, dots, vec));
155: PetscCall(VecNormalize(vec[i], NULL));
156: }
157: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)coords), PETSC_FALSE, nmodes, vec, sp));
158: for (i = 0; i < nmodes; i++) PetscCall(VecDestroy(&vec[i]));
159: }
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /*@
164: MatNullSpaceView - Visualizes a null space object.
166: Collective
168: Input Parameters:
169: + sp - the null space
170: - viewer - visualization context
172: Level: advanced
174: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `PetscViewer`, `MatNullSpaceCreate()`, `MatNullSpaceLoad()`, `PetscViewerASCIIOpen()`
175: @*/
176: PetscErrorCode MatNullSpaceView(MatNullSpace sp, PetscViewer viewer)
177: {
178: PetscBool isascii, isbinary;
180: PetscFunctionBegin;
182: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
184: PetscCheckSameComm(sp, 1, viewer, 2);
186: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
187: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
188: if (isascii) {
189: PetscViewerFormat format;
191: PetscCall(PetscViewerGetFormat(viewer, &format));
192: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sp, viewer));
193: PetscCall(PetscViewerASCIIPushTab(viewer));
194: PetscCall(PetscViewerASCIIPrintf(viewer, "Contains %" PetscInt_FMT " vector%s%s\n", sp->n, sp->n == 1 ? "" : "s", sp->has_cnst ? " and the constant" : ""));
195: if (sp->remove) PetscCall(PetscViewerASCIIPrintf(viewer, "Has user-provided removal function\n"));
196: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
197: for (PetscInt i = 0; i < sp->n; i++) PetscCall(VecView(sp->vecs[i], viewer));
198: }
199: PetscCall(PetscViewerASCIIPopTab(viewer));
200: } else if (isbinary) {
201: PetscBool skipHeader;
203: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
204: if (!skipHeader) {
205: PetscInt tr = MAT_NULLSPACE_FILE_CLASSID;
207: PetscCall(PetscViewerBinaryWrite(viewer, &tr, 1, PETSC_INT));
208: }
209: PetscCall(PetscViewerBinaryWrite(viewer, &sp->has_cnst, 1, PETSC_BOOL));
210: PetscCall(PetscViewerBinaryWrite(viewer, &sp->n, 1, PETSC_INT));
211: PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE));
212: for (PetscInt i = 0; i < sp->n; i++) PetscCall(VecView(sp->vecs[i], viewer));
213: PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader));
214: }
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@
219: MatNullSpaceLoad - Loads a `MatNullSpace` from a `PETSCVIEWERBINARY` that was saved by `MatNullSpaceView()`.
221: Collective
223: Input Parameter:
224: . viewer - the binary viewer
226: Output Parameter:
227: . sp - the null space
229: Level: advanced
231: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `PetscViewer`, `MatNullSpaceCreate()`, `MatNullSpaceView()`, `MatNullSpaceDestroy()`, `PetscViewerBinaryOpen()`, `PETSCVIEWERBINARY`
232: @*/
233: PetscErrorCode MatNullSpaceLoad(PetscViewer viewer, MatNullSpace *sp)
234: {
235: MPI_Comm comm;
236: Vec *vecs;
237: PetscBool has_cnst, isbinary, skipHeader;
238: PetscInt n;
240: PetscFunctionBegin;
242: PetscAssertPointer(sp, 2);
243: comm = PetscObjectComm((PetscObject)viewer);
244: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
245: PetscCheck(isbinary, comm, PETSC_ERR_SUP, "MatNullSpaceLoad() only supports binary viewers");
246: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
247: if (!skipHeader) {
248: PetscInt tr;
250: PetscCall(PetscViewerBinaryRead(viewer, &tr, 1, NULL, PETSC_INT));
251: PetscCheck(tr == MAT_NULLSPACE_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a MatNullSpace next in file");
252: }
253: PetscCall(PetscViewerBinaryRead(viewer, &has_cnst, 1, NULL, PETSC_BOOL));
254: PetscCall(PetscViewerBinaryRead(viewer, &n, 1, NULL, PETSC_INT));
255: PetscCheck(n >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Number of null space vectors in file (%" PetscInt_FMT ") cannot be negative", n);
256: PetscCall(PetscMalloc1(n, &vecs));
257: PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE));
258: for (PetscInt i = 0; i < n; i++) {
259: PetscCall(VecCreate(comm, &vecs[i]));
260: PetscCall(VecLoad(vecs[i], viewer));
261: }
262: PetscCall(MatNullSpaceCreate(comm, has_cnst, n, vecs, sp));
263: for (PetscInt i = 0; i < n; i++) PetscCall(VecDestroy(&vecs[i]));
264: PetscCall(PetscFree(vecs));
265: PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: MatNullSpaceCreate - Creates a `MatNullSpace` data structure used to project vectors out of null spaces.
272: Collective
274: Input Parameters:
275: + comm - the MPI communicator associated with the object
276: . has_cnst - `PETSC_TRUE` if the null space contains the constant vector; otherwise `PETSC_FALSE`
277: . n - number of vectors (excluding constant vector) in null space
278: - vecs - the vectors that span the null space (excluding the constant vector);
279: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
280: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
281: for them by one).
283: Output Parameter:
284: . SP - the null space context
286: Level: advanced
288: Notes:
289: See `MatNullSpaceSetFunction()` as an alternative way of providing the null space information instead of providing the vectors.
291: 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
292: need to pass in a function that eliminates the constant function into `MatNullSpaceSetFunction()`.
294: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpaceSetFunction()`
295: @*/
296: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm, PetscBool has_cnst, PetscInt n, const Vec vecs[], MatNullSpace *SP)
297: {
298: MatNullSpace sp;
299: PetscInt i;
301: PetscFunctionBegin;
302: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", n);
303: if (n) PetscAssertPointer(vecs, 4);
305: PetscAssertPointer(SP, 5);
306: if (n) {
307: for (i = 0; i < n; i++) {
308: /* prevent the user from changes values in the vector */
309: PetscCall(VecLockReadPush(vecs[i]));
310: }
311: }
312: if (PetscUnlikelyDebug(n)) {
313: PetscScalar *dots;
314: for (i = 0; i < n; i++) {
315: PetscReal norm;
316: PetscCall(VecNorm(vecs[i], NORM_2, &norm));
317: 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);
318: }
319: if (has_cnst) {
320: for (i = 0; i < n; i++) {
321: PetscScalar sum;
322: PetscCall(VecSum(vecs[i], &sum));
323: 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));
324: }
325: }
326: PetscCall(PetscMalloc1(n - 1, &dots));
327: for (i = 0; i < n - 1; i++) {
328: PetscInt j;
329: PetscCall(VecMDot(vecs[i], n - i - 1, vecs + i + 1, dots));
330: for (j = 0; j < n - i - 1; j++) {
331: 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]));
332: }
333: }
334: PetscCall(PetscFree(dots));
335: }
337: *SP = NULL;
338: PetscCall(MatInitializePackage());
340: PetscCall(PetscHeaderCreate(sp, MAT_NULLSPACE_CLASSID, "MatNullSpace", "Null space", "Mat", comm, MatNullSpaceDestroy, MatNullSpaceView));
342: sp->has_cnst = has_cnst;
343: sp->n = n;
344: sp->vecs = NULL;
345: sp->alpha = NULL;
346: sp->remove = NULL;
347: sp->rmctx = NULL;
349: if (n) {
350: PetscCall(PetscMalloc1(n, &sp->vecs));
351: PetscCall(PetscMalloc1(n, &sp->alpha));
352: for (i = 0; i < n; i++) {
353: PetscCall(PetscObjectReference((PetscObject)vecs[i]));
354: sp->vecs[i] = vecs[i];
355: }
356: }
358: *SP = sp;
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: /*@
363: MatNullSpaceDestroy - Destroys a data structure used to project vectors out of null spaces.
365: Collective
367: Input Parameter:
368: . sp - the null space context to be destroyed
370: Level: advanced
372: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceRemove()`, `MatNullSpaceSetFunction()`
373: @*/
374: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
375: {
376: PetscFunctionBegin;
377: if (!*sp) PetscFunctionReturn(PETSC_SUCCESS);
379: if (--((PetscObject)*sp)->refct > 0) {
380: *sp = NULL;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: for (PetscInt i = 0; i < (*sp)->n; i++) PetscCall(VecLockReadPop((*sp)->vecs[i]));
386: PetscCall(VecDestroyVecs((*sp)->n, &(*sp)->vecs));
387: PetscCall(PetscFree((*sp)->alpha));
388: PetscCall(PetscHeaderDestroy(sp));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*@
393: MatNullSpaceRemove - Removes all the components of a null space from a vector.
395: Collective
397: Input Parameters:
398: + sp - the null space context (if this is `NULL` then no null space is removed)
399: - vec - the vector from which the null space is to be removed
401: Level: advanced
403: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
404: @*/
405: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
406: {
407: PetscScalar sum;
408: PetscInt N;
410: PetscFunctionBegin;
411: if (!sp) PetscFunctionReturn(PETSC_SUCCESS);
415: if (sp->has_cnst) {
416: PetscCall(VecGetSize(vec, &N));
417: if (N > 0) {
418: PetscCall(VecSum(vec, &sum));
419: sum = sum / ((PetscScalar)(-1.0 * N));
420: PetscCall(VecShift(vec, sum));
421: }
422: }
424: if (sp->n) {
425: PetscCall(VecMDot(vec, sp->n, sp->vecs, sp->alpha));
426: for (PetscInt i = 0; i < sp->n; i++) sp->alpha[i] = -sp->alpha[i];
427: PetscCall(VecMAXPY(vec, sp->n, sp->alpha, sp->vecs));
428: }
430: if (sp->remove) PetscCall((*sp->remove)(sp, vec, sp->rmctx));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: MatNullSpaceTest - Tests if the claimed null space is really a null space of a matrix
437: Collective
439: Input Parameters:
440: + sp - the null space context
441: - mat - the matrix
443: Output Parameter:
444: . isNull - `PETSC_TRUE` if the nullspace is valid for this matrix
446: Level: advanced
448: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
449: @*/
450: PetscErrorCode MatNullSpaceTest(MatNullSpace sp, Mat mat, PetscBool *isNull)
451: {
452: PetscScalar sum;
453: PetscReal nrm, tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
454: PetscInt j, n, N;
455: Vec l, r;
456: PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE, consistent = PETSC_TRUE;
457: PetscViewer viewer;
459: PetscFunctionBegin;
462: n = sp->n;
463: PetscCall(PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view", &flg1, NULL));
464: PetscCall(PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view_draw", &flg2, NULL));
466: if (n) PetscCall(VecDuplicate(sp->vecs[0], &l));
467: else PetscCall(MatCreateVecs(mat, &l, NULL));
469: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
470: if (sp->has_cnst) {
471: PetscCall(VecDuplicate(l, &r));
472: PetscCall(VecGetSize(l, &N));
473: sum = 1.0 / PetscSqrtReal(N);
474: PetscCall(VecSet(l, sum));
475: PetscCall(MatMult(mat, l, r));
476: PetscCall(VecNorm(r, NORM_2, &nrm));
477: if (nrm >= tol) consistent = PETSC_FALSE;
478: if (flg1) {
479: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are %s null vector ", consistent ? "likely" : "unlikely"));
480: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * 1/sqrt(N) || = %g\n", (double)nrm));
481: }
482: if (!consistent && (flg1 || flg2)) PetscCall(VecView(r, viewer));
483: PetscCall(VecDestroy(&r));
484: }
486: for (j = 0; j < n; j++) {
487: PetscUseTypeMethod(mat, mult, sp->vecs[j], l);
488: PetscCall(VecNorm(l, NORM_2, &nrm));
489: if (nrm >= tol) consistent = PETSC_FALSE;
490: if (flg1) {
491: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " is %s null vector ", j, consistent ? "likely" : "unlikely"));
492: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * v[%" PetscInt_FMT "] || = %g\n", j, (double)nrm));
493: }
494: if (!consistent && (flg1 || flg2)) PetscCall(VecView(l, viewer));
495: }
497: PetscCheck(!sp->remove, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
498: PetscCall(VecDestroy(&l));
499: if (isNull) *isNull = consistent;
500: PetscFunctionReturn(PETSC_SUCCESS);
501: }