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: {
28: sp->remove = rem;
29: sp->rmctx = ctx;
30: return 0;
31: }
33: /*@C
34: MatNullSpaceGetVecs - get 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: Notes:
49: These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller
51: .seealso: 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 Vec
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: 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 MatNullSpace
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: MatNullSpaceCreate(), PetscViewerASCIIOpen()
173: @*/
174: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
175: {
176: PetscBool iascii;
179: if (!viewer) {
180: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
181: }
185: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
186: if (iascii) {
187: PetscViewerFormat format;
188: PetscInt i;
189: PetscViewerGetFormat(viewer,&format);
190: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
191: PetscViewerASCIIPushTab(viewer);
192: PetscViewerASCIIPrintf(viewer,"Contains %" PetscInt_FMT " vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
193: if (sp->remove) PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");
194: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
195: for (i=0; i<sp->n; i++) {
196: VecView(sp->vecs[i],viewer);
197: }
198: }
199: PetscViewerASCIIPopTab(viewer);
200: }
201: return 0;
202: }
204: /*@C
205: MatNullSpaceCreate - Creates a data structure used to project vectors
206: out of null spaces.
208: Collective
210: Input Parameters:
211: + comm - the MPI communicator associated with the object
212: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
213: . n - number of vectors (excluding constant vector) in null space
214: - vecs - the vectors that span the null space (excluding the constant vector);
215: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
216: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
217: for them by one).
219: Output Parameter:
220: . SP - the null space context
222: Level: advanced
224: Notes:
225: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.
227: 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
228: need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().
230: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
231: @*/
232: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
233: {
234: MatNullSpace sp;
235: PetscInt i;
241: if (n) {
242: for (i=0; i<n; i++) {
243: /* prevent the user from changes values in the vector */
244: VecLockReadPush(vecs[i]);
245: }
246: }
247: if (PetscUnlikelyDebug(n)) {
248: PetscScalar *dots;
249: for (i=0; i<n; i++) {
250: PetscReal norm;
251: VecNorm(vecs[i],NORM_2,&norm);
253: }
254: if (has_cnst) {
255: for (i=0; i<n; i++) {
256: PetscScalar sum;
257: VecSum(vecs[i],&sum);
259: }
260: }
261: PetscMalloc1(n-1,&dots);
262: for (i=0; i<n-1; i++) {
263: PetscInt j;
264: VecMDot(vecs[i],n-i-1,vecs+i+1,dots);
265: for (j=0;j<n-i-1;j++) {
267: }
268: }
269: PetscFree(dots);
270: }
272: *SP = NULL;
273: MatInitializePackage();
275: PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
277: sp->has_cnst = has_cnst;
278: sp->n = n;
279: sp->vecs = NULL;
280: sp->alpha = NULL;
281: sp->remove = NULL;
282: sp->rmctx = NULL;
284: if (n) {
285: PetscMalloc1(n,&sp->vecs);
286: PetscMalloc1(n,&sp->alpha);
287: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
288: for (i=0; i<n; i++) {
289: PetscObjectReference((PetscObject)vecs[i]);
290: sp->vecs[i] = vecs[i];
291: }
292: }
294: *SP = sp;
295: return 0;
296: }
298: /*@
299: MatNullSpaceDestroy - Destroys a data structure used to project vectors
300: out of null spaces.
302: Collective on MatNullSpace
304: Input Parameter:
305: . sp - the null space context to be destroyed
307: Level: advanced
309: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
310: @*/
311: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
312: {
313: PetscInt i;
315: if (!*sp) return 0;
317: if (--((PetscObject)(*sp))->refct > 0) {*sp = NULL; return 0;}
319: for (i=0; i < (*sp)->n; i++) {
320: VecLockReadPop((*sp)->vecs[i]);
321: }
323: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
324: PetscFree((*sp)->alpha);
325: PetscHeaderDestroy(sp);
326: return 0;
327: }
329: /*@C
330: MatNullSpaceRemove - Removes all the components of a null space from a vector.
332: Collective on MatNullSpace
334: Input Parameters:
335: + sp - the null space context (if this is NULL then no null space is removed)
336: - vec - the vector from which the null space is to be removed
338: Level: advanced
340: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
341: @*/
342: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
343: {
344: PetscScalar sum;
345: PetscInt i,N;
347: if (!sp) return 0;
351: if (sp->has_cnst) {
352: VecGetSize(vec,&N);
353: if (N > 0) {
354: VecSum(vec,&sum);
355: sum = sum/((PetscScalar)(-1.0*N));
356: VecShift(vec,sum);
357: }
358: }
360: if (sp->n) {
361: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
362: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
363: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
364: }
366: if (sp->remove) {
367: (*sp->remove)(sp,vec,sp->rmctx);
368: }
369: return 0;
370: }
372: /*@
373: MatNullSpaceTest - Tests if the claimed null space is really a
374: null space of a matrix
376: Collective on MatNullSpace
378: Input Parameters:
379: + sp - the null space context
380: - mat - the matrix
382: Output Parameters:
383: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
385: Level: advanced
387: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
388: @*/
389: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
390: {
391: PetscScalar sum;
392: PetscReal nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
393: PetscInt j,n,N;
394: Vec l,r;
395: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
396: PetscViewer viewer;
400: n = sp->n;
401: PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view",&flg1,NULL);
402: PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view_draw",&flg2,NULL);
404: if (n) {
405: VecDuplicate(sp->vecs[0],&l);
406: } else {
407: MatCreateVecs(mat,&l,NULL);
408: }
410: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
411: if (sp->has_cnst) {
412: VecDuplicate(l,&r);
413: VecGetSize(l,&N);
414: sum = 1.0/PetscSqrtReal(N);
415: VecSet(l,sum);
416: MatMult(mat,l,r);
417: VecNorm(r,NORM_2,&nrm);
418: if (nrm >= tol) consistent = PETSC_FALSE;
419: if (flg1) {
420: if (consistent) {
421: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
422: } else {
423: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
424: }
425: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
426: }
427: if (!consistent && flg1) VecView(r,viewer);
428: if (!consistent && flg2) VecView(r,viewer);
429: VecDestroy(&r);
430: }
432: for (j=0; j<n; j++) {
433: (*mat->ops->mult)(mat,sp->vecs[j],l);
434: VecNorm(l,NORM_2,&nrm);
435: if (nrm >= tol) consistent = PETSC_FALSE;
436: if (flg1) {
437: if (consistent) {
438: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %" PetscInt_FMT " is likely null vector",j);
439: } else {
440: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %" PetscInt_FMT " unlikely null vector ",j);
441: consistent = PETSC_FALSE;
442: }
443: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%" PetscInt_FMT "] || = %g\n",j,(double)nrm);
444: }
445: if (!consistent && flg1) VecView(l,viewer);
446: if (!consistent && flg2) VecView(l,viewer);
447: }
450: VecDestroy(&l);
451: if (isNull) *isNull = consistent;
452: return 0;
453: }