Actual source code: matnull.c
petsc-3.11.4 2019-09-28
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 MatNullSpace
16: Input Parameters:
17: + sp - the null space object
18: . rem - the function that removes the null space
19: - ctx - context for the remove function
21: Level: advanced
23: .keywords: PC, null space, create
25: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
26: @*/
27: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace,Vec,void*),void *ctx)
28: {
31: sp->remove = rem;
32: sp->rmctx = ctx;
33: return(0);
34: }
36: /*@C
37: MatNullSpaceGetVecs - get vectors defining the null space
39: Not Collective
41: Input Arguments:
42: . sp - null space object
44: Output Arguments:
45: + has_cnst - PETSC_TRUE if the null space contains the constant vector, otherwise PETSC_FALSE
46: . n - number of vectors (excluding constant vector) in null space
47: - vecs - orthonormal vectors that span the null space (excluding the constant vector)
49: Level: developer
51: Notes:
52: These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller
54: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
55: @*/
56: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
57: {
61: if (has_const) *has_const = sp->has_cnst;
62: if (n) *n = sp->n;
63: if (vecs) *vecs = sp->vecs;
64: return(0);
65: }
67: /*@
68: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
70: Collective on Vec
72: Input Argument:
73: . coords - block of coordinates of each node, must have block size set
75: Output Argument:
76: . sp - the null space
78: Level: advanced
80: Notes:
81: If you are solving an elasticity problem you should likely use this, in conjunction with MatSetNearNullspace(), to provide information that
82: the PCGAMG preconditioner can use to construct a much more efficient preconditioner.
84: If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with MatSetNullspace() to
85: provide this information to the linear solver so it can handle the null space appropriately in the linear solution.
88: .seealso: MatNullSpaceCreate(), MatSetNearNullspace(), MatSetNullspace()
89: @*/
90: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
91: {
92: PetscErrorCode ierr;
93: const PetscScalar *x;
94: PetscScalar *v[6],dots[5];
95: Vec vec[6];
96: PetscInt n,N,dim,nmodes,i,j;
97: PetscReal sN;
100: VecGetBlockSize(coords,&dim);
101: VecGetLocalSize(coords,&n);
102: VecGetSize(coords,&N);
103: n /= dim;
104: N /= dim;
105: sN = 1./PetscSqrtReal((PetscReal)N);
106: switch (dim) {
107: case 1:
108: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);
109: break;
110: case 2:
111: case 3:
112: nmodes = (dim == 2) ? 3 : 6;
113: VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);
114: VecSetSizes(vec[0],dim*n,dim*N);
115: VecSetBlockSize(vec[0],dim);
116: VecSetUp(vec[0]);
117: for (i=1; i<nmodes; i++) {VecDuplicate(vec[0],&vec[i]);}
118: for (i=0; i<nmodes; i++) {VecGetArray(vec[i],&v[i]);}
119: VecGetArrayRead(coords,&x);
120: for (i=0; i<n; i++) {
121: if (dim == 2) {
122: v[0][i*2+0] = sN;
123: v[0][i*2+1] = 0.;
124: v[1][i*2+0] = 0.;
125: v[1][i*2+1] = sN;
126: /* Rotations */
127: v[2][i*2+0] = -x[i*2+1];
128: v[2][i*2+1] = x[i*2+0];
129: } else {
130: v[0][i*3+0] = sN;
131: v[0][i*3+1] = 0.;
132: v[0][i*3+2] = 0.;
133: v[1][i*3+0] = 0.;
134: v[1][i*3+1] = sN;
135: v[1][i*3+2] = 0.;
136: v[2][i*3+0] = 0.;
137: v[2][i*3+1] = 0.;
138: v[2][i*3+2] = sN;
140: v[3][i*3+0] = x[i*3+1];
141: v[3][i*3+1] = -x[i*3+0];
142: v[3][i*3+2] = 0.;
143: v[4][i*3+0] = 0.;
144: v[4][i*3+1] = -x[i*3+2];
145: v[4][i*3+2] = x[i*3+1];
146: v[5][i*3+0] = x[i*3+2];
147: v[5][i*3+1] = 0.;
148: v[5][i*3+2] = -x[i*3+0];
149: }
150: }
151: for (i=0; i<nmodes; i++) {VecRestoreArray(vec[i],&v[i]);}
152: VecRestoreArrayRead(coords,&x);
153: for (i=dim; i<nmodes; i++) {
154: /* Orthonormalize vec[i] against vec[0:i-1] */
155: VecMDot(vec[i],i,vec,dots);
156: for (j=0; j<i; j++) dots[j] *= -1.;
157: VecMAXPY(vec[i],i,dots,vec);
158: VecNormalize(vec[i],NULL);
159: }
160: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);
161: for (i=0; i<nmodes; i++) {VecDestroy(&vec[i]);}
162: }
163: return(0);
164: }
166: /*@C
167: MatNullSpaceView - Visualizes a null space object.
169: Collective on MatNullSpace
171: Input Parameters:
172: + matnull - the null space
173: - viewer - visualization context
175: Level: advanced
177: Fortran Note:
178: This routine is not supported in Fortran.
180: .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen()
181: @*/
182: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
183: {
185: PetscBool iascii;
189: if (!viewer) {
190: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
191: }
195: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
196: if (iascii) {
197: PetscViewerFormat format;
198: PetscInt i;
199: PetscViewerGetFormat(viewer,&format);
200: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
201: PetscViewerASCIIPushTab(viewer);
202: PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
203: if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
204: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
205: for (i=0; i<sp->n; i++) {
206: VecView(sp->vecs[i],viewer);
207: }
208: }
209: PetscViewerASCIIPopTab(viewer);
210: }
211: return(0);
212: }
214: /*@C
215: MatNullSpaceCreate - Creates a data structure used to project vectors
216: out of null spaces.
218: Collective on MPI_Comm
220: Input Parameters:
221: + comm - the MPI communicator associated with the object
222: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
223: . n - number of vectors (excluding constant vector) in null space
224: - vecs - the vectors that span the null space (excluding the constant vector);
225: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
226: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
227: for them by one).
229: Output Parameter:
230: . SP - the null space context
232: Level: advanced
234: Notes:
235: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.
237: 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
238: need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().
240: Users manual sections:
241: . Section 4.19 Solving Singular Systems
243: .keywords: PC, null space, create
245: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
246: @*/
247: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
248: {
249: MatNullSpace sp;
251: PetscInt i;
254: if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
258: if (n) {
259: for (i=0; i<n; i++) {
260: /* prevent the user from changes values in the vector */
261: VecLockReadPush(vecs[i]);
262: }
263: }
264: #if defined(PETSC_USE_DEBUG)
265: if (n) {
266: PetscScalar *dots;
267: for (i=0; i<n; i++) {
268: PetscReal norm;
269: VecNorm(vecs[i],NORM_2,&norm);
270: if (PetscAbsReal(norm - 1) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D must have 2-norm of 1.0, it is %g",i,(double)norm);
271: }
272: if (has_cnst) {
273: for (i=0; i<n; i++) {
274: PetscScalar sum;
275: VecSum(vecs[i],&sum);
276: if (PetscAbsScalar(sum) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D must be orthogonal to constant vector, inner product is %g",i,(double)PetscAbsScalar(sum));
277: }
278: }
279: PetscMalloc1(n-1,&dots);
280: for (i=0; i<n-1; i++) {
281: PetscInt j;
282: VecMDot(vecs[i],n-i-1,vecs+i+1,dots);
283: for (j=0;j<n-i-1;j++) {
284: if (PetscAbsScalar(dots[j]) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ3(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D must be orthogonal to vector %D, inner product is %g",i,i+j+1,(double)PetscAbsScalar(dots[j]));
285: }
286: }
287: PetscFree(dots);
288: }
289: #endif
291: *SP = NULL;
292: MatInitializePackage();
294: PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
296: sp->has_cnst = has_cnst;
297: sp->n = n;
298: sp->vecs = 0;
299: sp->alpha = 0;
300: sp->remove = 0;
301: sp->rmctx = 0;
303: if (n) {
304: PetscMalloc1(n,&sp->vecs);
305: PetscMalloc1(n,&sp->alpha);
306: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
307: for (i=0; i<n; i++) {
308: PetscObjectReference((PetscObject)vecs[i]);
309: sp->vecs[i] = vecs[i];
310: }
311: }
313: *SP = sp;
314: return(0);
315: }
317: /*@
318: MatNullSpaceDestroy - Destroys a data structure used to project vectors
319: out of null spaces.
321: Collective on MatNullSpace
323: Input Parameter:
324: . sp - the null space context to be destroyed
326: Level: advanced
328: .keywords: PC, null space, destroy
330: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
331: @*/
332: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
333: {
335: PetscInt i;
338: if (!*sp) return(0);
340: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
342: for (i=0; i < (*sp)->n; i++) {
343: VecLockReadPop((*sp)->vecs[i]);
344: }
346: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
347: PetscFree((*sp)->alpha);
348: PetscHeaderDestroy(sp);
349: return(0);
350: }
352: /*@C
353: MatNullSpaceRemove - Removes all the components of a null space from a vector.
355: Collective on MatNullSpace
357: Input Parameters:
358: + sp - the null space context (if this is NULL then no null space is removed)
359: - vec - the vector from which the null space is to be removed
361: Level: advanced
363: .keywords: PC, null space, remove
365: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
366: @*/
367: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
368: {
369: PetscScalar sum;
370: PetscInt i,N;
374: if (!sp) return(0);
378: if (sp->has_cnst) {
379: VecGetSize(vec,&N);
380: if (N > 0) {
381: VecSum(vec,&sum);
382: sum = sum/((PetscScalar)(-1.0*N));
383: VecShift(vec,sum);
384: }
385: }
387: if (sp->n) {
388: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
389: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
390: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
391: }
393: if (sp->remove) {
394: (*sp->remove)(sp,vec,sp->rmctx);
395: }
396: return(0);
397: }
399: /*@
400: MatNullSpaceTest - Tests if the claimed null space is really a
401: null space of a matrix
403: Collective on MatNullSpace
405: Input Parameters:
406: + sp - the null space context
407: - mat - the matrix
409: Output Parameters:
410: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
412: Level: advanced
414: .keywords: PC, null space, remove
416: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
417: @*/
418: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
419: {
420: PetscScalar sum;
421: PetscReal nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
422: PetscInt j,n,N;
424: Vec l,r;
425: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
426: PetscViewer viewer;
431: n = sp->n;
432: PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view",&flg1,NULL);
433: PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view_draw",&flg2,NULL);
435: if (n) {
436: VecDuplicate(sp->vecs[0],&l);
437: } else {
438: MatCreateVecs(mat,&l,NULL);
439: }
441: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
442: if (sp->has_cnst) {
443: VecDuplicate(l,&r);
444: VecGetSize(l,&N);
445: sum = 1.0/N;
446: VecSet(l,sum);
447: MatMult(mat,l,r);
448: VecNorm(r,NORM_2,&nrm);
449: if (nrm >= tol) consistent = PETSC_FALSE;
450: if (flg1) {
451: if (consistent) {
452: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
453: } else {
454: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
455: }
456: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
457: }
458: if (!consistent && flg1) {VecView(r,viewer);}
459: if (!consistent && flg2) {VecView(r,viewer);}
460: VecDestroy(&r);
461: }
463: for (j=0; j<n; j++) {
464: (*mat->ops->mult)(mat,sp->vecs[j],l);
465: VecNorm(l,NORM_2,&nrm);
466: if (nrm >= tol) consistent = PETSC_FALSE;
467: if (flg1) {
468: if (consistent) {
469: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
470: } else {
471: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
472: consistent = PETSC_FALSE;
473: }
474: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
475: }
476: if (!consistent && flg1) {VecView(l,viewer);}
477: if (!consistent && flg2) {VecView(l,viewer);}
478: }
480: if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
481: VecDestroy(&l);
482: if (isNull) *isNull = consistent;
483: return(0);
484: }