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: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are %s null vector ", consistent ? "likely" : "unlikely"));
416: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * 1/sqrt(N) || = %g\n", (double)nrm));
417: }
418: if (!consistent && (flg1 || flg2)) PetscCall(VecView(r, viewer));
419: PetscCall(VecDestroy(&r));
420: }
422: for (j = 0; j < n; j++) {
423: PetscUseTypeMethod(mat, mult, sp->vecs[j], l);
424: PetscCall(VecNorm(l, NORM_2, &nrm));
425: if (nrm >= tol) consistent = PETSC_FALSE;
426: if (flg1) {
427: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " is %s null vector ", j, consistent ? "likely" : "unlikely"));
428: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * v[%" PetscInt_FMT "] || = %g\n", j, (double)nrm));
429: }
430: if (!consistent && (flg1 || flg2)) PetscCall(VecView(l, viewer));
431: }
433: PetscCheck(!sp->remove, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
434: PetscCall(VecDestroy(&l));
435: if (isNull) *isNull = consistent;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }