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 on sp
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 on coords
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 on sp
163: Input Parameters:
164: + matnull - the null space
165: - viewer - visualization context
167: Level: advanced
169: Fortran Note:
170: This routine is not supported in Fortran.
172: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `PetscViewerASCIIOpen()`
173: @*/
174: PetscErrorCode MatNullSpaceView(MatNullSpace sp, PetscViewer viewer)
175: {
176: PetscBool iascii;
179: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer);
183: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
184: if (iascii) {
185: PetscViewerFormat format;
186: PetscInt i;
187: PetscViewerGetFormat(viewer, &format);
188: PetscObjectPrintClassNamePrefixType((PetscObject)sp, viewer);
189: PetscViewerASCIIPushTab(viewer);
190: 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) 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++) VecView(sp->vecs[i], viewer);
194: }
195: PetscViewerASCIIPopTab(viewer);
196: }
197: return 0;
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 setting vecs.
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: `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpace`, `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;
236: if (n) {
237: for (i = 0; i < n; i++) {
238: /* prevent the user from changes values in the vector */
239: VecLockReadPush(vecs[i]);
240: }
241: }
242: if (PetscUnlikelyDebug(n)) {
243: PetscScalar *dots;
244: for (i = 0; i < n; i++) {
245: PetscReal norm;
246: VecNorm(vecs[i], NORM_2, &norm);
248: }
249: if (has_cnst) {
250: for (i = 0; i < n; i++) {
251: PetscScalar sum;
252: VecSum(vecs[i], &sum);
254: }
255: }
256: PetscMalloc1(n - 1, &dots);
257: for (i = 0; i < n - 1; i++) {
258: PetscInt j;
259: VecMDot(vecs[i], n - i - 1, vecs + i + 1, dots);
260: for (j = 0; j < n - i - 1; j++) {
262: }
263: }
264: PetscFree(dots);
265: }
267: *SP = NULL;
268: MatInitializePackage();
270: PetscHeaderCreate(sp, MAT_NULLSPACE_CLASSID, "MatNullSpace", "Null space", "Mat", comm, MatNullSpaceDestroy, MatNullSpaceView);
272: sp->has_cnst = has_cnst;
273: sp->n = n;
274: sp->vecs = NULL;
275: sp->alpha = NULL;
276: sp->remove = NULL;
277: sp->rmctx = NULL;
279: if (n) {
280: PetscMalloc1(n, &sp->vecs);
281: PetscMalloc1(n, &sp->alpha);
282: for (i = 0; i < n; i++) {
283: PetscObjectReference((PetscObject)vecs[i]);
284: sp->vecs[i] = vecs[i];
285: }
286: }
288: *SP = sp;
289: return 0;
290: }
292: /*@
293: MatNullSpaceDestroy - Destroys a data structure used to project vectors out of null spaces.
295: Collective on sp
297: Input Parameter:
298: . sp - the null space context to be destroyed
300: Level: advanced
302: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceRemove()`, `MatNullSpaceSetFunction()`
303: @*/
304: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
305: {
306: PetscInt i;
308: if (!*sp) return 0;
310: if (--((PetscObject)(*sp))->refct > 0) {
311: *sp = NULL;
312: return 0;
313: }
315: for (i = 0; i < (*sp)->n; i++) VecLockReadPop((*sp)->vecs[i]);
317: VecDestroyVecs((*sp)->n, &(*sp)->vecs);
318: PetscFree((*sp)->alpha);
319: PetscHeaderDestroy(sp);
320: return 0;
321: }
323: /*@C
324: MatNullSpaceRemove - Removes all the components of a null space from a vector.
326: Collective on sp
328: Input Parameters:
329: + sp - the null space context (if this is NULL then no null space is removed)
330: - vec - the vector from which the null space is to be removed
332: Level: advanced
334: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
335: @*/
336: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
337: {
338: PetscScalar sum;
339: PetscInt i, N;
341: if (!sp) return 0;
345: if (sp->has_cnst) {
346: VecGetSize(vec, &N);
347: if (N > 0) {
348: VecSum(vec, &sum);
349: sum = sum / ((PetscScalar)(-1.0 * N));
350: VecShift(vec, sum);
351: }
352: }
354: if (sp->n) {
355: VecMDot(vec, sp->n, sp->vecs, sp->alpha);
356: for (i = 0; i < sp->n; i++) sp->alpha[i] = -sp->alpha[i];
357: VecMAXPY(vec, sp->n, sp->alpha, sp->vecs);
358: }
360: if (sp->remove) (*sp->remove)(sp, vec, sp->rmctx);
361: return 0;
362: }
364: /*@
365: MatNullSpaceTest - Tests if the claimed null space is really a null space of a matrix
367: Collective on sp
369: Input Parameters:
370: + sp - the null space context
371: - mat - the matrix
373: Output Parameters:
374: . isNull - `PETSC_TRUE` if the nullspace is valid for this matrix
376: Level: advanced
378: .seealso: `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
379: @*/
380: PetscErrorCode MatNullSpaceTest(MatNullSpace sp, Mat mat, PetscBool *isNull)
381: {
382: PetscScalar sum;
383: PetscReal nrm, tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
384: PetscInt j, n, N;
385: Vec l, r;
386: PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE, consistent = PETSC_TRUE;
387: PetscViewer viewer;
391: n = sp->n;
392: PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view", &flg1, NULL);
393: PetscOptionsGetBool(((PetscObject)sp)->options, ((PetscObject)mat)->prefix, "-mat_null_space_test_view_draw", &flg2, NULL);
395: if (n) {
396: VecDuplicate(sp->vecs[0], &l);
397: } else {
398: MatCreateVecs(mat, &l, NULL);
399: }
401: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer);
402: if (sp->has_cnst) {
403: VecDuplicate(l, &r);
404: VecGetSize(l, &N);
405: sum = 1.0 / PetscSqrtReal(N);
406: VecSet(l, sum);
407: MatMult(mat, l, r);
408: VecNorm(r, NORM_2, &nrm);
409: if (nrm >= tol) consistent = PETSC_FALSE;
410: if (flg1) {
411: if (consistent) {
412: PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are likely null vector");
413: } else {
414: PetscPrintf(PetscObjectComm((PetscObject)sp), "Constants are unlikely null vector ");
415: }
416: PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * 1/N || = %g\n", (double)nrm);
417: }
418: if (!consistent && flg1) VecView(r, viewer);
419: if (!consistent && flg2) VecView(r, viewer);
420: VecDestroy(&r);
421: }
423: for (j = 0; j < n; j++) {
424: (*mat->ops->mult)(mat, sp->vecs[j], l);
425: VecNorm(l, NORM_2, &nrm);
426: if (nrm >= tol) consistent = PETSC_FALSE;
427: if (flg1) {
428: if (consistent) {
429: PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " is likely null vector", j);
430: } else {
431: PetscPrintf(PetscObjectComm((PetscObject)sp), "Null vector %" PetscInt_FMT " unlikely null vector ", j);
432: consistent = PETSC_FALSE;
433: }
434: PetscPrintf(PetscObjectComm((PetscObject)sp), "|| A * v[%" PetscInt_FMT "] || = %g\n", j, (double)nrm);
435: }
436: if (!consistent && flg1) VecView(l, viewer);
437: if (!consistent && flg2) VecView(l, viewer);
438: }
441: VecDestroy(&l);
442: if (isNull) *isNull = consistent;
443: return 0;
444: }