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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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