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: }