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: `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpace`, `MatNullSpaceCreate()`
24: @*/
25: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace, Vec, void *), void *ctx)
26: {
28: sp->remove = rem;
29: sp->rmctx = ctx;
30: return 0;
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_cnst - `PETSC_TRUE` if the null space contains the constant vector, otherwise `PETSC_FALSE`
43: . n - number of vectors (excluding constant vector) in null space
44: - vecs - orthonormal vectors that span the null space (excluding the constant vector)
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: `MatNullSpace`, `MatNullSpaceCreate()`, `MatGetNullSpace()`, `MatGetNearNullSpace()`
52: @*/
53: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp, PetscBool *has_const, PetscInt *n, const Vec **vecs)
54: {
56: if (has_const) *has_const = sp->has_cnst;
57: if (n) *n = sp->n;
58: if (vecs) *vecs = sp->vecs;
59: return 0;
60: }
62: /*@
63: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
65: Collective
67: Input Parameter:
68: . coords - block of coordinates of each node, must have block size set
70: Output Parameter:
71: . sp - the null space
73: Level: advanced
75: Notes:
76: If you are solving an elasticity problem you should likely use this, in conjunction with `MatSetNearNullSpace()`, to provide information that
77: the `PCGAMG` preconditioner can use to construct a much more efficient preconditioner.
79: If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with `MatSetNullSpace()` to
80: provide this information to the linear solver so it can handle the null space appropriately in the linear solution.
82: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatSetNearNullSpace()`, `MatSetNullSpace()`
83: @*/
84: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords, MatNullSpace *sp)
85: {
86: const PetscScalar *x;
87: PetscScalar *v[6], dots[5];
88: Vec vec[6];
89: PetscInt n, N, dim, nmodes, i, j;
90: PetscReal sN;
92: VecGetBlockSize(coords, &dim);
93: VecGetLocalSize(coords, &n);
94: VecGetSize(coords, &N);
95: n /= dim;
96: N /= dim;
97: sN = 1. / PetscSqrtReal((PetscReal)N);
98: switch (dim) {
99: case 1:
100: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords), PETSC_TRUE, 0, NULL, sp);
101: break;
102: case 2:
103: case 3:
104: nmodes = (dim == 2) ? 3 : 6;
105: VecCreate(PetscObjectComm((PetscObject)coords), &vec[0]);
106: VecSetSizes(vec[0], dim * n, dim * N);
107: VecSetBlockSize(vec[0], dim);
108: VecSetUp(vec[0]);
109: for (i = 1; i < nmodes; i++) VecDuplicate(vec[0], &vec[i]);
110: for (i = 0; i < nmodes; i++) VecGetArray(vec[i], &v[i]);
111: VecGetArrayRead(coords, &x);
112: for (i = 0; i < n; i++) {
113: if (dim == 2) {
114: v[0][i * 2 + 0] = sN;
115: v[0][i * 2 + 1] = 0.;
116: v[1][i * 2 + 0] = 0.;
117: v[1][i * 2 + 1] = sN;
118: /* Rotations */
119: v[2][i * 2 + 0] = -x[i * 2 + 1];
120: v[2][i * 2 + 1] = x[i * 2 + 0];
121: } else {
122: v[0][i * 3 + 0] = sN;
123: v[0][i * 3 + 1] = 0.;
124: v[0][i * 3 + 2] = 0.;
125: v[1][i * 3 + 0] = 0.;
126: v[1][i * 3 + 1] = sN;
127: v[1][i * 3 + 2] = 0.;
128: v[2][i * 3 + 0] = 0.;
129: v[2][i * 3 + 1] = 0.;
130: v[2][i * 3 + 2] = sN;
132: v[3][i * 3 + 0] = x[i * 3 + 1];
133: v[3][i * 3 + 1] = -x[i * 3 + 0];
134: v[3][i * 3 + 2] = 0.;
135: v[4][i * 3 + 0] = 0.;
136: v[4][i * 3 + 1] = -x[i * 3 + 2];
137: v[4][i * 3 + 2] = x[i * 3 + 1];
138: v[5][i * 3 + 0] = x[i * 3 + 2];
139: v[5][i * 3 + 1] = 0.;
140: v[5][i * 3 + 2] = -x[i * 3 + 0];
141: }
142: }
143: for (i = 0; i < nmodes; i++) VecRestoreArray(vec[i], &v[i]);
144: VecRestoreArrayRead(coords, &x);
145: for (i = dim; i < nmodes; i++) {
146: /* Orthonormalize vec[i] against vec[0:i-1] */
147: VecMDot(vec[i], i, vec, dots);
148: for (j = 0; j < i; j++) dots[j] *= -1.;
149: VecMAXPY(vec[i], i, dots, vec);
150: VecNormalize(vec[i], NULL);
151: }
152: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords), PETSC_FALSE, nmodes, vec, sp);
153: for (i = 0; i < nmodes; i++) VecDestroy(&vec[i]);
154: }
155: return 0;
156: }
158: /*@C
159: MatNullSpaceView - Visualizes a null space object.
161: Collective; No Fortran Support
163: Input Parameters:
164: + matnull - the null space
165: - viewer - visualization context
167: Level: advanced
169: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `PetscViewerASCIIOpen()`
170: @*/
171: PetscErrorCode MatNullSpaceView(MatNullSpace sp, PetscViewer viewer)
172: {
173: PetscBool iascii;
176: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer);
180: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
181: if (iascii) {
182: PetscViewerFormat format;
183: PetscInt i;
184: PetscViewerGetFormat(viewer, &format);
185: PetscObjectPrintClassNamePrefixType((PetscObject)sp, viewer);
186: PetscViewerASCIIPushTab(viewer);
187: PetscViewerASCIIPrintf(viewer, "Contains %" PetscInt_FMT " vector%s%s\n", sp->n, sp->n == 1 ? "" : "s", sp->has_cnst ? " and the constant" : "");
188: if (sp->remove) PetscViewerASCIIPrintf(viewer, "Has user-provided removal function\n");
189: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
190: for (i = 0; i < sp->n; i++) VecView(sp->vecs[i], viewer);
191: }
192: PetscViewerASCIIPopTab(viewer);
193: }
194: return 0;
195: }
197: /*@C
198: MatNullSpaceCreate - Creates a `MatNullSpace` data structure used to project vectors out of null spaces.
200: Collective
202: Input Parameters:
203: + comm - the MPI communicator associated with the object
204: . has_cnst - `PETSC_TRUE` if the null space contains the constant vector; otherwise `PETSC_FALSE`
205: . n - number of vectors (excluding constant vector) in null space
206: - vecs - the vectors that span the null space (excluding the constant vector);
207: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
208: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
209: for them by one).
211: Output Parameter:
212: . SP - the null space context
214: Level: advanced
216: Notes:
217: See `MatNullSpaceSetFunction()` as an alternative way of providing the null space information instead of setting vecs.
219: 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
220: need to pass in a function that eliminates the constant function into `MatNullSpaceSetFunction()`.
222: .seealso: `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpace`, `MatNullSpaceSetFunction()`
223: @*/
224: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm, PetscBool has_cnst, PetscInt n, const Vec vecs[], MatNullSpace *SP)
225: {
226: MatNullSpace sp;
227: PetscInt i;
233: if (n) {
234: for (i = 0; i < n; i++) {
235: /* prevent the user from changes values in the vector */
236: VecLockReadPush(vecs[i]);
237: }
238: }
239: if (PetscUnlikelyDebug(n)) {
240: PetscScalar *dots;
241: for (i = 0; i < n; i++) {
242: PetscReal norm;
243: VecNorm(vecs[i], NORM_2, &norm);
245: }
246: if (has_cnst) {
247: for (i = 0; i < n; i++) {
248: PetscScalar sum;
249: VecSum(vecs[i], &sum);
251: }
252: }
253: PetscMalloc1(n - 1, &dots);
254: for (i = 0; i < n - 1; i++) {
255: PetscInt j;
256: VecMDot(vecs[i], n - i - 1, vecs + i + 1, dots);
257: for (j = 0; j < n - i - 1; j++) {
259: }
260: }
261: PetscFree(dots);
262: }
264: *SP = NULL;
265: MatInitializePackage();
267: PetscHeaderCreate(sp, MAT_NULLSPACE_CLASSID, "MatNullSpace", "Null space", "Mat", comm, MatNullSpaceDestroy, MatNullSpaceView);
269: sp->has_cnst = has_cnst;
270: sp->n = n;
271: sp->vecs = NULL;
272: sp->alpha = NULL;
273: sp->remove = NULL;
274: sp->rmctx = NULL;
276: if (n) {
277: PetscMalloc1(n, &sp->vecs);
278: PetscMalloc1(n, &sp->alpha);
279: for (i = 0; i < n; i++) {
280: PetscObjectReference((PetscObject)vecs[i]);
281: sp->vecs[i] = vecs[i];
282: }
283: }
285: *SP = sp;
286: return 0;
287: }
289: /*@
290: MatNullSpaceDestroy - Destroys a data structure used to project vectors out of null spaces.
292: Collective
294: Input Parameter:
295: . sp - the null space context to be destroyed
297: Level: advanced
299: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceRemove()`, `MatNullSpaceSetFunction()`
300: @*/
301: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
302: {
303: PetscInt i;
305: if (!*sp) return 0;
307: if (--((PetscObject)(*sp))->refct > 0) {
308: *sp = NULL;
309: return 0;
310: }
312: for (i = 0; i < (*sp)->n; i++) VecLockReadPop((*sp)->vecs[i]);
314: VecDestroyVecs((*sp)->n, &(*sp)->vecs);
315: PetscFree((*sp)->alpha);
316: PetscHeaderDestroy(sp);
317: return 0;
318: }
320: /*@C
321: MatNullSpaceRemove - Removes all the components of a null space from a vector.
323: Collective
325: Input Parameters:
326: + sp - the null space context (if this is NULL then no null space is removed)
327: - vec - the vector from which the null space is to be removed
329: Level: advanced
331: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
332: @*/
333: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
334: {
335: PetscScalar sum;
336: PetscInt i, N;
338: if (!sp) return 0;
342: if (sp->has_cnst) {
343: VecGetSize(vec, &N);
344: if (N > 0) {
345: VecSum(vec, &sum);
346: sum = sum / ((PetscScalar)(-1.0 * N));
347: VecShift(vec, sum);
348: }
349: }
351: if (sp->n) {
352: VecMDot(vec, sp->n, sp->vecs, sp->alpha);
353: for (i = 0; i < sp->n; i++) sp->alpha[i] = -sp->alpha[i];
354: VecMAXPY(vec, sp->n, sp->alpha, sp->vecs);
355: }
357: if (sp->remove) (*sp->remove)(sp, vec, sp->rmctx);
358: return 0;
359: }
361: /*@
362: MatNullSpaceTest - Tests if the claimed null space is really a null space of a matrix
364: Collective
366: Input Parameters:
367: + sp - the null space context
368: - mat - the matrix
370: Output Parameters:
371: . isNull - `PETSC_TRUE` if the nullspace is valid for this matrix
373: Level: advanced
375: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
376: @*/
377: PetscErrorCode MatNullSpaceTest(MatNullSpace sp, Mat mat, PetscBool *isNull)
378: {
379: PetscScalar sum;
380: PetscReal nrm, tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
381: PetscInt j, n, N;
382: Vec l, r;
383: PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE, consistent = PETSC_TRUE;
384: PetscViewer viewer;
388: n = sp->n;
389: PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view", &flg1, NULL);
390: PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view_draw", &flg2, NULL);
392: if (n) {
393: VecDuplicate(sp->vecs[0], &l);
394: } else {
395: MatCreateVecs(mat, &l, NULL);
396: }
398: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer);
399: if (sp->has_cnst) {
400: VecDuplicate(l, &r);
401: VecGetSize(l, &N);
402: sum = 1.0 / PetscSqrtReal(N);
403: VecSet(l, sum);
404: MatMult(mat, l, r);
405: VecNorm(r, NORM_2, &nrm);
406: if (nrm >= tol) consistent = PETSC_FALSE;
407: if (flg1) {
408: if (consistent) {
409: PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are likely null vector");
410: } else {
411: PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are unlikely null vector ");
412: }
413: PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * 1/N || = %g\n", (double)nrm);
414: }
415: if (!consistent && flg1) VecView(r, viewer);
416: if (!consistent && flg2) VecView(r, viewer);
417: VecDestroy(&r);
418: }
420: for (j = 0; j < n; j++) {
421: (*mat->ops->mult)(mat, sp->vecs[j], l);
422: VecNorm(l, NORM_2, &nrm);
423: if (nrm >= tol) consistent = PETSC_FALSE;
424: if (flg1) {
425: if (consistent) {
426: PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " is likely null vector", j);
427: } else {
428: PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " unlikely null vector ", j);
429: consistent = PETSC_FALSE;
430: }
431: PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * v[%" PetscInt_FMT "] || = %g\n", j, (double)nrm);
432: }
433: if (!consistent && flg1) VecView(l, viewer);
434: if (!consistent && flg2) VecView(l, viewer);
435: }
438: VecDestroy(&l);
439: if (isNull) *isNull = consistent;
440: return 0;
441: }