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 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: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
24: @*/
25: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace,Vec,void*),void *ctx)
26: {
29: sp->remove = rem;
30: sp->rmctx = ctx;
31: return(0);
32: }
34: /*@C
35: MatNullSpaceGetVecs - get vectors defining the null space
37: Not Collective
39: Input Parameter:
40: . sp - null space object
42: Output Parameters:
43: + has_cnst - PETSC_TRUE if the null space contains the constant vector, otherwise PETSC_FALSE
44: . n - number of vectors (excluding constant vector) in null space
45: - vecs - orthonormal vectors that span the null space (excluding the constant vector)
47: Level: developer
49: Notes:
50: These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller
52: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
53: @*/
54: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
55: {
59: if (has_const) *has_const = sp->has_cnst;
60: if (n) *n = sp->n;
61: if (vecs) *vecs = sp->vecs;
62: return(0);
63: }
65: /*@
66: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
68: Collective on Vec
70: Input Parameter:
71: . coords - block of coordinates of each node, must have block size set
73: Output Parameter:
74: . sp - the null space
76: Level: advanced
78: Notes:
79: If you are solving an elasticity problem you should likely use this, in conjunction with MatSetNearNullspace(), to provide information that
80: the PCGAMG preconditioner can use to construct a much more efficient preconditioner.
82: If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with MatSetNullspace() to
83: provide this information to the linear solver so it can handle the null space appropriately in the linear solution.
85: .seealso: MatNullSpaceCreate(), MatSetNearNullspace(), MatSetNullspace()
86: @*/
87: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
88: {
89: PetscErrorCode ierr;
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;
97: VecGetBlockSize(coords,&dim);
98: VecGetLocalSize(coords,&n);
99: VecGetSize(coords,&N);
100: n /= dim;
101: N /= dim;
102: sN = 1./PetscSqrtReal((PetscReal)N);
103: switch (dim) {
104: case 1:
105: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);
106: break;
107: case 2:
108: case 3:
109: nmodes = (dim == 2) ? 3 : 6;
110: VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);
111: VecSetSizes(vec[0],dim*n,dim*N);
112: VecSetBlockSize(vec[0],dim);
113: VecSetUp(vec[0]);
114: for (i=1; i<nmodes; i++) {VecDuplicate(vec[0],&vec[i]);}
115: for (i=0; i<nmodes; i++) {VecGetArray(vec[i],&v[i]);}
116: 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;
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++) {VecRestoreArray(vec[i],&v[i]);}
149: VecRestoreArrayRead(coords,&x);
150: for (i=dim; i<nmodes; i++) {
151: /* Orthonormalize vec[i] against vec[0:i-1] */
152: VecMDot(vec[i],i,vec,dots);
153: for (j=0; j<i; j++) dots[j] *= -1.;
154: VecMAXPY(vec[i],i,dots,vec);
155: VecNormalize(vec[i],NULL);
156: }
157: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);
158: for (i=0; i<nmodes; i++) {VecDestroy(&vec[i]);}
159: }
160: return(0);
161: }
163: /*@C
164: MatNullSpaceView - Visualizes a null space object.
166: Collective on MatNullSpace
168: Input Parameters:
169: + matnull - the null space
170: - viewer - visualization context
172: Level: advanced
174: Fortran Note:
175: This routine is not supported in Fortran.
177: .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen()
178: @*/
179: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
180: {
182: PetscBool iascii;
186: if (!viewer) {
187: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
188: }
192: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
193: if (iascii) {
194: PetscViewerFormat format;
195: PetscInt i;
196: PetscViewerGetFormat(viewer,&format);
197: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
198: PetscViewerASCIIPushTab(viewer);
199: PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
200: if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
201: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
202: for (i=0; i<sp->n; i++) {
203: VecView(sp->vecs[i],viewer);
204: }
205: }
206: PetscViewerASCIIPopTab(viewer);
207: }
208: return(0);
209: }
211: /*@C
212: MatNullSpaceCreate - Creates a data structure used to project vectors
213: out of null spaces.
215: Collective
217: Input Parameters:
218: + comm - the MPI communicator associated with the object
219: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
220: . n - number of vectors (excluding constant vector) in null space
221: - vecs - the vectors that span the null space (excluding the constant vector);
222: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
223: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
224: for them by one).
226: Output Parameter:
227: . SP - the null space context
229: Level: advanced
231: Notes:
232: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.
234: 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
235: need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().
237: Users manual sections:
238: . sec_singular
240: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
241: @*/
242: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
243: {
244: MatNullSpace sp;
246: PetscInt i;
249: if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
253: if (n) {
254: for (i=0; i<n; i++) {
255: /* prevent the user from changes values in the vector */
256: VecLockReadPush(vecs[i]);
257: }
258: }
259: if (PetscUnlikelyDebug(n)) {
260: PetscScalar *dots;
261: for (i=0; i<n; i++) {
262: PetscReal norm;
263: VecNorm(vecs[i],NORM_2,&norm);
264: 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);
265: }
266: if (has_cnst) {
267: for (i=0; i<n; i++) {
268: PetscScalar sum;
269: VecSum(vecs[i],&sum);
270: 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));
271: }
272: }
273: PetscMalloc1(n-1,&dots);
274: for (i=0; i<n-1; i++) {
275: PetscInt j;
276: VecMDot(vecs[i],n-i-1,vecs+i+1,dots);
277: for (j=0;j<n-i-1;j++) {
278: 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]));
279: }
280: }
281: PetscFree(dots);
282: }
284: *SP = NULL;
285: MatInitializePackage();
287: PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
289: sp->has_cnst = has_cnst;
290: sp->n = n;
291: sp->vecs = NULL;
292: sp->alpha = NULL;
293: sp->remove = NULL;
294: sp->rmctx = NULL;
296: if (n) {
297: PetscMalloc1(n,&sp->vecs);
298: PetscMalloc1(n,&sp->alpha);
299: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
300: for (i=0; i<n; i++) {
301: PetscObjectReference((PetscObject)vecs[i]);
302: sp->vecs[i] = vecs[i];
303: }
304: }
306: *SP = sp;
307: return(0);
308: }
310: /*@
311: MatNullSpaceDestroy - Destroys a data structure used to project vectors
312: out of null spaces.
314: Collective on MatNullSpace
316: Input Parameter:
317: . sp - the null space context to be destroyed
319: Level: advanced
321: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
322: @*/
323: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
324: {
326: PetscInt i;
329: if (!*sp) return(0);
331: if (--((PetscObject)(*sp))->refct > 0) {*sp = NULL; return(0);}
333: for (i=0; i < (*sp)->n; i++) {
334: VecLockReadPop((*sp)->vecs[i]);
335: }
337: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
338: PetscFree((*sp)->alpha);
339: PetscHeaderDestroy(sp);
340: return(0);
341: }
343: /*@C
344: MatNullSpaceRemove - Removes all the components of a null space from a vector.
346: Collective on MatNullSpace
348: Input Parameters:
349: + sp - the null space context (if this is NULL then no null space is removed)
350: - vec - the vector from which the null space is to be removed
352: Level: advanced
354: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
355: @*/
356: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
357: {
358: PetscScalar sum;
359: PetscInt i,N;
363: if (!sp) return(0);
367: if (sp->has_cnst) {
368: VecGetSize(vec,&N);
369: if (N > 0) {
370: VecSum(vec,&sum);
371: sum = sum/((PetscScalar)(-1.0*N));
372: VecShift(vec,sum);
373: }
374: }
376: if (sp->n) {
377: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
378: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
379: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
380: }
382: if (sp->remove) {
383: (*sp->remove)(sp,vec,sp->rmctx);
384: }
385: return(0);
386: }
388: /*@
389: MatNullSpaceTest - Tests if the claimed null space is really a
390: null space of a matrix
392: Collective on MatNullSpace
394: Input Parameters:
395: + sp - the null space context
396: - mat - the matrix
398: Output Parameters:
399: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
401: Level: advanced
403: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
404: @*/
405: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
406: {
407: PetscScalar sum;
408: PetscReal nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
409: PetscInt j,n,N;
411: Vec l,r;
412: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
413: PetscViewer viewer;
418: n = sp->n;
419: PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view",&flg1,NULL);
420: PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view_draw",&flg2,NULL);
422: if (n) {
423: VecDuplicate(sp->vecs[0],&l);
424: } else {
425: MatCreateVecs(mat,&l,NULL);
426: }
428: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
429: if (sp->has_cnst) {
430: VecDuplicate(l,&r);
431: VecGetSize(l,&N);
432: sum = 1.0/PetscSqrtReal(N);
433: VecSet(l,sum);
434: MatMult(mat,l,r);
435: VecNorm(r,NORM_2,&nrm);
436: if (nrm >= tol) consistent = PETSC_FALSE;
437: if (flg1) {
438: if (consistent) {
439: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
440: } else {
441: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
442: }
443: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
444: }
445: if (!consistent && flg1) {VecView(r,viewer);}
446: if (!consistent && flg2) {VecView(r,viewer);}
447: VecDestroy(&r);
448: }
450: for (j=0; j<n; j++) {
451: (*mat->ops->mult)(mat,sp->vecs[j],l);
452: VecNorm(l,NORM_2,&nrm);
453: if (nrm >= tol) consistent = PETSC_FALSE;
454: if (flg1) {
455: if (consistent) {
456: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
457: } else {
458: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
459: consistent = PETSC_FALSE;
460: }
461: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
462: }
463: if (!consistent && flg1) {VecView(l,viewer);}
464: if (!consistent && flg2) {VecView(l,viewer);}
465: }
467: if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
468: VecDestroy(&l);
469: if (isNull) *isNull = consistent;
470: return(0);
471: }