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: Fortran Note:
52: One must pass in an array `vecs` that is large enough to hold all of the requested vectors
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: /*@C
164: MatNullSpaceView - Visualizes a null space object.
166: Collective; No Fortran Support
168: Input Parameters:
169: + sp - the null space
170: - viewer - visualization context
172: Level: advanced
174: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `PetscViewer`, `MatNullSpaceCreate()`, `PetscViewerASCIIOpen()`
175: @*/
176: PetscErrorCode MatNullSpaceView(MatNullSpace sp, PetscViewer viewer)
177: {
178: PetscBool iascii;
180: PetscFunctionBegin;
182: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
184: PetscCheckSameComm(sp, 1, viewer, 2);
186: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
187: if (iascii) {
188: PetscViewerFormat format;
189: PetscInt i;
190: PetscCall(PetscViewerGetFormat(viewer, &format));
191: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sp, viewer));
192: PetscCall(PetscViewerASCIIPushTab(viewer));
193: PetscCall(PetscViewerASCIIPrintf(viewer, "Contains %" PetscInt_FMT " vector%s%s\n", sp->n, sp->n == 1 ? "" : "s", sp->has_cnst ? " and the constant" : ""));
194: if (sp->remove) PetscCall(PetscViewerASCIIPrintf(viewer, "Has user-provided removal function\n"));
195: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
196: for (i = 0; i < sp->n; i++) PetscCall(VecView(sp->vecs[i], viewer));
197: }
198: PetscCall(PetscViewerASCIIPopTab(viewer));
199: }
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@C
204: MatNullSpaceCreate - Creates a `MatNullSpace` data structure used to project vectors out of null spaces.
206: Collective
208: Input Parameters:
209: + comm - the MPI communicator associated with the object
210: . has_cnst - `PETSC_TRUE` if the null space contains the constant vector; otherwise `PETSC_FALSE`
211: . n - number of vectors (excluding constant vector) in null space
212: - vecs - the vectors that span the null space (excluding the constant vector);
213: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
214: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
215: for them by one).
217: Output Parameter:
218: . SP - the null space context
220: Level: advanced
222: Notes:
223: See `MatNullSpaceSetFunction()` as an alternative way of providing the null space information instead of providing the vectors.
225: 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
226: need to pass in a function that eliminates the constant function into `MatNullSpaceSetFunction()`.
228: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpaceSetFunction()`
229: @*/
230: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm, PetscBool has_cnst, PetscInt n, const Vec vecs[], MatNullSpace *SP)
231: {
232: MatNullSpace sp;
233: PetscInt i;
235: PetscFunctionBegin;
236: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", n);
237: if (n) PetscAssertPointer(vecs, 4);
239: PetscAssertPointer(SP, 5);
240: if (n) {
241: for (i = 0; i < n; i++) {
242: /* prevent the user from changes values in the vector */
243: PetscCall(VecLockReadPush(vecs[i]));
244: }
245: }
246: if (PetscUnlikelyDebug(n)) {
247: PetscScalar *dots;
248: for (i = 0; i < n; i++) {
249: PetscReal norm;
250: PetscCall(VecNorm(vecs[i], NORM_2, &norm));
251: 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);
252: }
253: if (has_cnst) {
254: for (i = 0; i < n; i++) {
255: PetscScalar sum;
256: PetscCall(VecSum(vecs[i], &sum));
257: 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));
258: }
259: }
260: PetscCall(PetscMalloc1(n - 1, &dots));
261: for (i = 0; i < n - 1; i++) {
262: PetscInt j;
263: PetscCall(VecMDot(vecs[i], n - i - 1, vecs + i + 1, dots));
264: for (j = 0; j < n - i - 1; j++) {
265: 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]));
266: }
267: }
268: PetscCall(PetscFree(dots));
269: }
271: *SP = NULL;
272: PetscCall(MatInitializePackage());
274: PetscCall(PetscHeaderCreate(sp, MAT_NULLSPACE_CLASSID, "MatNullSpace", "Null space", "Mat", comm, MatNullSpaceDestroy, MatNullSpaceView));
276: sp->has_cnst = has_cnst;
277: sp->n = n;
278: sp->vecs = NULL;
279: sp->alpha = NULL;
280: sp->remove = NULL;
281: sp->rmctx = NULL;
283: if (n) {
284: PetscCall(PetscMalloc1(n, &sp->vecs));
285: PetscCall(PetscMalloc1(n, &sp->alpha));
286: for (i = 0; i < n; i++) {
287: PetscCall(PetscObjectReference((PetscObject)vecs[i]));
288: sp->vecs[i] = vecs[i];
289: }
290: }
292: *SP = sp;
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: MatNullSpaceDestroy - Destroys a data structure used to project vectors out of null spaces.
299: Collective
301: Input Parameter:
302: . sp - the null space context to be destroyed
304: Level: advanced
306: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceRemove()`, `MatNullSpaceSetFunction()`
307: @*/
308: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
309: {
310: PetscInt i;
312: PetscFunctionBegin;
313: if (!*sp) PetscFunctionReturn(PETSC_SUCCESS);
315: if (--((PetscObject)*sp)->refct > 0) {
316: *sp = NULL;
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: for (i = 0; i < (*sp)->n; i++) PetscCall(VecLockReadPop((*sp)->vecs[i]));
322: PetscCall(VecDestroyVecs((*sp)->n, &(*sp)->vecs));
323: PetscCall(PetscFree((*sp)->alpha));
324: PetscCall(PetscHeaderDestroy(sp));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@C
329: MatNullSpaceRemove - Removes all the components of a null space from a vector.
331: Collective
333: Input Parameters:
334: + sp - the null space context (if this is `NULL` then no null space is removed)
335: - vec - the vector from which the null space is to be removed
337: Level: advanced
339: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
340: @*/
341: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
342: {
343: PetscScalar sum;
344: PetscInt i, N;
346: PetscFunctionBegin;
347: if (!sp) PetscFunctionReturn(PETSC_SUCCESS);
351: if (sp->has_cnst) {
352: PetscCall(VecGetSize(vec, &N));
353: if (N > 0) {
354: PetscCall(VecSum(vec, &sum));
355: sum = sum / ((PetscScalar)(-1.0 * N));
356: PetscCall(VecShift(vec, sum));
357: }
358: }
360: if (sp->n) {
361: PetscCall(VecMDot(vec, sp->n, sp->vecs, sp->alpha));
362: for (i = 0; i < sp->n; i++) sp->alpha[i] = -sp->alpha[i];
363: PetscCall(VecMAXPY(vec, sp->n, sp->alpha, sp->vecs));
364: }
366: if (sp->remove) PetscCall((*sp->remove)(sp, vec, sp->rmctx));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: /*@
371: MatNullSpaceTest - Tests if the claimed null space is really a null space of a matrix
373: Collective
375: Input Parameters:
376: + sp - the null space context
377: - mat - the matrix
379: Output Parameter:
380: . isNull - `PETSC_TRUE` if the nullspace is valid for this matrix
382: Level: advanced
384: .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
385: @*/
386: PetscErrorCode MatNullSpaceTest(MatNullSpace sp, Mat mat, PetscBool *isNull)
387: {
388: PetscScalar sum;
389: PetscReal nrm, tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
390: PetscInt j, n, N;
391: Vec l, r;
392: PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE, consistent = PETSC_TRUE;
393: PetscViewer viewer;
395: PetscFunctionBegin;
398: n = sp->n;
399: PetscCall(PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view", &flg1, NULL));
400: PetscCall(PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view_draw", &flg2, NULL));
402: if (n) {
403: PetscCall(VecDuplicate(sp->vecs[0], &l));
404: } else {
405: PetscCall(MatCreateVecs(mat, &l, NULL));
406: }
408: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
409: if (sp->has_cnst) {
410: PetscCall(VecDuplicate(l, &r));
411: PetscCall(VecGetSize(l, &N));
412: sum = 1.0 / PetscSqrtReal(N);
413: PetscCall(VecSet(l, sum));
414: PetscCall(MatMult(mat, l, r));
415: PetscCall(VecNorm(r, NORM_2, &nrm));
416: if (nrm >= tol) consistent = PETSC_FALSE;
417: if (flg1) {
418: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are %s null vector ", consistent ? "likely" : "unlikely"));
419: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * 1/sqrt(N) || = %g\n", (double)nrm));
420: }
421: if (!consistent && (flg1 || flg2)) PetscCall(VecView(r, viewer));
422: PetscCall(VecDestroy(&r));
423: }
425: for (j = 0; j < n; j++) {
426: PetscUseTypeMethod(mat, mult, sp->vecs[j], l);
427: PetscCall(VecNorm(l, NORM_2, &nrm));
428: if (nrm >= tol) consistent = PETSC_FALSE;
429: if (flg1) {
430: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " is %s null vector ", j, consistent ? "likely" : "unlikely"));
431: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * v[%" PetscInt_FMT "] || = %g\n", j, (double)nrm));
432: }
433: if (!consistent && (flg1 || flg2)) PetscCall(VecView(l, viewer));
434: }
436: PetscCheck(!sp->remove, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
437: PetscCall(VecDestroy(&l));
438: if (isNull) *isNull = consistent;
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }