1 /*
2 Routines to project vectors out of null spaces.
3 */
4
5 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
6
7 PetscClassId MAT_NULLSPACE_CLASSID;
8
9 /*@C
10 MatNullSpaceSetFunction - set a function that removes a null space from a vector
11 out of null spaces.
12
13 Logically Collective
14
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
19
20 Level: advanced
21
22 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpaceCreate()`, `MatNullSpaceRemoveFn`
23 @*/
MatNullSpaceSetFunction(MatNullSpace sp,MatNullSpaceRemoveFn * rem,PetscCtx ctx)24 PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, MatNullSpaceRemoveFn *rem, PetscCtx ctx)
25 {
26 PetscFunctionBegin;
27 PetscValidHeaderSpecific(sp, MAT_NULLSPACE_CLASSID, 1);
28 sp->remove = rem;
29 sp->rmctx = ctx;
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
33 /*@C
34 MatNullSpaceGetVecs - get the vectors defining the null space
35
36 Not Collective
37
38 Input Parameter:
39 . sp - null space object
40
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 - returns array of length `n` containing the orthonormal vectors that span the null space (excluding the constant vector), `NULL` if `n` is 0
45
46 Level: developer
47
48 Note:
49 These vectors and the array returned are owned by the `MatNullSpace` and should not be destroyed or freeded by the caller
50
51 Fortran Note:
52 Call `MatNullSpaceRestoreVecs()` when the array of `Vec` is no longer needed
53
54 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatGetNullSpace()`, `MatGetNearNullSpace()`
55 @*/
MatNullSpaceGetVecs(MatNullSpace sp,PetscBool * has_const,PetscInt * n,const Vec * vecs[])56 PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp, PetscBool *has_const, PetscInt *n, const Vec *vecs[])
57 {
58 PetscFunctionBegin;
59 PetscValidHeaderSpecific(sp, MAT_NULLSPACE_CLASSID, 1);
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 }
65
66 /*@
67 MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
68
69 Collective
70
71 Input Parameter:
72 . coords - block of coordinates of each node, must have block size set
73
74 Output Parameter:
75 . sp - the null space
76
77 Level: advanced
78
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.
82
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.
85
86 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatSetNearNullSpace()`, `MatSetNullSpace()`, `PCGAMG`
87 @*/
MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace * sp)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;
95
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;
136
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 }
162
163 /*@
164 MatNullSpaceView - Visualizes a null space object.
165
166 Collective
167
168 Input Parameters:
169 + sp - the null space
170 - viewer - visualization context
171
172 Level: advanced
173
174 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `PetscViewer`, `MatNullSpaceCreate()`, `PetscViewerASCIIOpen()`
175 @*/
MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)176 PetscErrorCode MatNullSpaceView(MatNullSpace sp, PetscViewer viewer)
177 {
178 PetscBool isascii;
179
180 PetscFunctionBegin;
181 PetscValidHeaderSpecific(sp, MAT_NULLSPACE_CLASSID, 1);
182 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp), &viewer));
183 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
184 PetscCheckSameComm(sp, 1, viewer, 2);
185
186 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
187 if (isascii) {
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 }
202
203 /*@
204 MatNullSpaceCreate - Creates a `MatNullSpace` data structure used to project vectors out of null spaces.
205
206 Collective
207
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).
216
217 Output Parameter:
218 . SP - the null space context
219
220 Level: advanced
221
222 Notes:
223 See `MatNullSpaceSetFunction()` as an alternative way of providing the null space information instead of providing the vectors.
224
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()`.
227
228 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceDestroy()`, `MatNullSpaceRemove()`, `MatSetNullSpace()`, `MatNullSpaceSetFunction()`
229 @*/
MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace * SP)230 PetscErrorCode MatNullSpaceCreate(MPI_Comm comm, PetscBool has_cnst, PetscInt n, const Vec vecs[], MatNullSpace *SP)
231 {
232 MatNullSpace sp;
233 PetscInt i;
234
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);
238 for (i = 0; i < n; i++) PetscValidHeaderSpecific(vecs[i], VEC_CLASSID, 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 }
270
271 *SP = NULL;
272 PetscCall(MatInitializePackage());
273
274 PetscCall(PetscHeaderCreate(sp, MAT_NULLSPACE_CLASSID, "MatNullSpace", "Null space", "Mat", comm, MatNullSpaceDestroy, MatNullSpaceView));
275
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;
282
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 }
291
292 *SP = sp;
293 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
296 /*@
297 MatNullSpaceDestroy - Destroys a data structure used to project vectors out of null spaces.
298
299 Collective
300
301 Input Parameter:
302 . sp - the null space context to be destroyed
303
304 Level: advanced
305
306 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceRemove()`, `MatNullSpaceSetFunction()`
307 @*/
MatNullSpaceDestroy(MatNullSpace * sp)308 PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
309 {
310 PetscInt i;
311
312 PetscFunctionBegin;
313 if (!*sp) PetscFunctionReturn(PETSC_SUCCESS);
314 PetscValidHeaderSpecific(*sp, MAT_NULLSPACE_CLASSID, 1);
315 if (--((PetscObject)*sp)->refct > 0) {
316 *sp = NULL;
317 PetscFunctionReturn(PETSC_SUCCESS);
318 }
319
320 for (i = 0; i < (*sp)->n; i++) PetscCall(VecLockReadPop((*sp)->vecs[i]));
321
322 PetscCall(VecDestroyVecs((*sp)->n, &(*sp)->vecs));
323 PetscCall(PetscFree((*sp)->alpha));
324 PetscCall(PetscHeaderDestroy(sp));
325 PetscFunctionReturn(PETSC_SUCCESS);
326 }
327
328 /*@
329 MatNullSpaceRemove - Removes all the components of a null space from a vector.
330
331 Collective
332
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
336
337 Level: advanced
338
339 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
340 @*/
MatNullSpaceRemove(MatNullSpace sp,Vec vec)341 PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
342 {
343 PetscScalar sum;
344 PetscInt i, N;
345
346 PetscFunctionBegin;
347 if (!sp) PetscFunctionReturn(PETSC_SUCCESS);
348 PetscValidHeaderSpecific(sp, MAT_NULLSPACE_CLASSID, 1);
349 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
350
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 }
359
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 }
365
366 if (sp->remove) PetscCall((*sp->remove)(sp, vec, sp->rmctx));
367 PetscFunctionReturn(PETSC_SUCCESS);
368 }
369
370 /*@
371 MatNullSpaceTest - Tests if the claimed null space is really a null space of a matrix
372
373 Collective
374
375 Input Parameters:
376 + sp - the null space context
377 - mat - the matrix
378
379 Output Parameter:
380 . isNull - `PETSC_TRUE` if the nullspace is valid for this matrix
381
382 Level: advanced
383
384 .seealso: [](ch_matrices), `Mat`, `MatNullSpace`, `MatNullSpaceCreate()`, `MatNullSpaceDestroy()`, `MatNullSpaceSetFunction()`
385 @*/
MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool * isNull)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;
394
395 PetscFunctionBegin;
396 PetscValidHeaderSpecific(sp, MAT_NULLSPACE_CLASSID, 1);
397 PetscValidHeaderSpecific(mat, MAT_CLASSID, 2);
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));
401
402 if (n) {
403 PetscCall(VecDuplicate(sp->vecs[0], &l));
404 } else {
405 PetscCall(MatCreateVecs(mat, &l, NULL));
406 }
407
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 }
424
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 }
435
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 }
441