Actual source code: matnull.c
petsc-3.6.1 2015-08-06
2: /*
3: Routines to project vectors out of null spaces.
4: */
6: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
8: PetscClassId MAT_NULLSPACE_CLASSID;
12: /*@C
13: MatNullSpaceSetFunction - set a function that removes a null space from a vector
14: out of null spaces.
16: Logically Collective on MatNullSpace
18: Input Parameters:
19: + sp - the null space object
20: . rem - the function that removes the null space
21: - ctx - context for the remove function
23: Level: advanced
25: .keywords: PC, null space, create
27: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
28: @*/
29: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace,Vec,void*),void *ctx)
30: {
33: sp->remove = rem;
34: sp->rmctx = ctx;
35: return(0);
36: }
40: /*@C
41: MatNullSpaceGetVecs - get vectors defining the null space
43: Not Collective
45: Input Arguments:
46: . sp - null space object
48: Output Arguments:
49: + has_cnst - PETSC_TRUE if the null space contains the constant vector, otherwise PETSC_FALSE
50: . n - number of vectors (excluding constant vector) in null space
51: - vecs - orthonormal vectors that span the null space (excluding the constant vector)
53: Level: developer
55: Notes:
56: These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller
58: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
59: @*/
60: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
61: {
65: if (has_const) *has_const = sp->has_cnst;
66: if (n) *n = sp->n;
67: if (vecs) *vecs = sp->vecs;
68: return(0);
69: }
73: /*@
74: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
76: Collective on Vec
78: Input Argument:
79: . coords - block of coordinates of each node, must have block size set
81: Output Argument:
82: . sp - the null space
84: Level: advanced
86: .seealso: MatNullSpaceCreate()
87: @*/
88: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
89: {
90: PetscErrorCode ierr;
91: const PetscScalar *x;
92: PetscScalar *v[6],dots[5];
93: Vec vec[6];
94: PetscInt n,N,dim,nmodes,i,j;
97: VecGetBlockSize(coords,&dim);
98: VecGetLocalSize(coords,&n);
99: VecGetSize(coords,&N);
100: n /= dim;
101: N /= dim;
102: switch (dim) {
103: case 1:
104: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);
105: break;
106: case 2:
107: case 3:
108: nmodes = (dim == 2) ? 3 : 6;
109: VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);
110: VecSetSizes(vec[0],dim*n,dim*N);
111: VecSetBlockSize(vec[0],dim);
112: VecSetUp(vec[0]);
113: for (i=1; i<nmodes; i++) {VecDuplicate(vec[0],&vec[i]);}
114: for (i=0; i<nmodes; i++) {VecGetArray(vec[i],&v[i]);}
115: VecGetArrayRead(coords,&x);
116: for (i=0; i<n; i++) {
117: if (dim == 2) {
118: v[0][i*2+0] = 1./N;
119: v[0][i*2+1] = 0.;
120: v[1][i*2+0] = 0.;
121: v[1][i*2+1] = 1./N;
122: /* Rotations */
123: v[2][i*2+0] = -x[i*2+1];
124: v[2][i*2+1] = x[i*2+0];
125: } else {
126: v[0][i*3+0] = 1./N;
127: v[0][i*3+1] = 0.;
128: v[0][i*3+2] = 0.;
129: v[1][i*3+0] = 0.;
130: v[1][i*3+1] = 1./N;
131: v[1][i*3+2] = 0.;
132: v[2][i*3+0] = 0.;
133: v[2][i*3+1] = 0.;
134: v[2][i*3+2] = 1./N;
136: v[3][i*3+0] = x[i*3+1];
137: v[3][i*3+1] = -x[i*3+0];
138: v[3][i*3+2] = 0.;
139: v[4][i*3+0] = 0.;
140: v[4][i*3+1] = -x[i*3+2];
141: v[4][i*3+2] = x[i*3+1];
142: v[5][i*3+0] = x[i*3+2];
143: v[5][i*3+1] = 0.;
144: v[5][i*3+2] = -x[i*3+0];
145: }
146: }
147: for (i=0; i<nmodes; i++) {VecRestoreArray(vec[i],&v[i]);}
148: VecRestoreArrayRead(coords,&x);
149: for (i=dim; i<nmodes; i++) {
150: /* Orthonormalize vec[i] against vec[0:i-1] */
151: VecMDot(vec[i],i,vec,dots);
152: for (j=0; j<i; j++) dots[j] *= -1.;
153: VecMAXPY(vec[i],i,dots,vec);
154: VecNormalize(vec[i],NULL);
155: }
156: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);
157: for (i=0; i<nmodes; i++) {VecDestroy(&vec[i]);}
158: }
159: return(0);
160: }
164: /*@C
165: MatNullSpaceView - Visualizes a null space object.
167: Collective on MatNullSpace
169: Input Parameters:
170: + matnull - the null space
171: - viewer - visualization context
173: Level: advanced
175: Fortran Note:
176: This routine is not supported in Fortran.
178: .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen()
179: @*/
180: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
181: {
183: PetscBool iascii;
187: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)sp));
191: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
192: if (iascii) {
193: PetscViewerFormat format;
194: PetscInt i;
195: PetscViewerGetFormat(viewer,&format);
196: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
197: PetscViewerASCIIPushTab(viewer);
198: PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
199: if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
200: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
201: for (i=0; i<sp->n; i++) {
202: VecView(sp->vecs[i],viewer);
203: }
204: }
205: PetscViewerASCIIPopTab(viewer);
206: }
207: return(0);
208: }
212: /*@
213: MatNullSpaceCreate - Creates a data structure used to project vectors
214: out of null spaces.
216: Collective on MPI_Comm
218: Input Parameters:
219: + comm - the MPI communicator associated with the object
220: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
221: . n - number of vectors (excluding constant vector) in null space
222: - vecs - the vectors that span the null space (excluding the constant vector);
223: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
224: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
225: for them by one).
227: Output Parameter:
228: . SP - the null space context
230: Level: advanced
232: Notes: 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: . Section 4.19 Solving Singular Systems
240: .keywords: PC, null space, create
242: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
243: @*/
244: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
245: {
246: MatNullSpace sp;
248: PetscInt i;
251: if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
256: *SP = NULL;
257: MatInitializePackage();
259: PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
261: sp->has_cnst = has_cnst;
262: sp->n = n;
263: sp->vecs = 0;
264: sp->alpha = 0;
265: sp->remove = 0;
266: sp->rmctx = 0;
268: if (n) {
269: PetscMalloc1(n,&sp->vecs);
270: PetscMalloc1(n,&sp->alpha);
271: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
272: for (i=0; i<n; i++) {
273: PetscObjectReference((PetscObject)vecs[i]);
274: sp->vecs[i] = vecs[i];
275: }
276: }
278: *SP = sp;
279: return(0);
280: }
284: /*@
285: MatNullSpaceDestroy - Destroys a data structure used to project vectors
286: out of null spaces.
288: Collective on MatNullSpace
290: Input Parameter:
291: . sp - the null space context to be destroyed
293: Level: advanced
295: .keywords: PC, null space, destroy
297: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
298: @*/
299: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
300: {
304: if (!*sp) return(0);
306: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
308: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
309: PetscFree((*sp)->alpha);
310: PetscHeaderDestroy(sp);
311: return(0);
312: }
316: /*@C
317: MatNullSpaceRemove - Removes all the components of a null space from a vector.
319: Collective on MatNullSpace
321: Input Parameters:
322: + sp - the null space context
323: - vec - the vector from which the null space is to be removed
325: Level: advanced
327: .keywords: PC, null space, remove
329: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
330: @*/
331: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
332: {
333: PetscScalar sum;
334: PetscInt i,N;
341: if (sp->has_cnst) {
342: VecGetSize(vec,&N);
343: if (N > 0) {
344: VecSum(vec,&sum);
345: sum = sum/((PetscScalar)(-1.0*N));
346: VecShift(vec,sum);
347: }
348: }
350: if (sp->n) {
351: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
352: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
353: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
354: }
356: if (sp->remove) {
357: (*sp->remove)(sp,vec,sp->rmctx);
358: }
359: return(0);
360: }
364: /*@
365: MatNullSpaceTest - Tests if the claimed null space is really a
366: null space of a matrix
368: Collective on MatNullSpace
370: Input Parameters:
371: + sp - the null space context
372: - mat - the matrix
374: Output Parameters:
375: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
377: Level: advanced
379: .keywords: PC, null space, remove
381: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
382: @*/
383: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
384: {
385: PetscScalar sum;
386: PetscReal nrm;
387: PetscInt j,n,N;
389: Vec l,r;
390: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
391: PetscViewer viewer;
396: n = sp->n;
397: PetscOptionsGetBool(NULL,"-mat_null_space_test_view",&flg1,NULL);
398: PetscOptionsGetBool(NULL,"-mat_null_space_test_view_draw",&flg2,NULL);
400: if (n) {
401: VecDuplicate(sp->vecs[0],&l);
402: } else {
403: MatCreateVecs(mat,&l,NULL);
404: }
406: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
407: if (sp->has_cnst) {
408: VecDuplicate(l,&r);
409: VecGetSize(l,&N);
410: sum = 1.0/N;
411: VecSet(l,sum);
412: MatMult(mat,l,r);
413: VecNorm(r,NORM_2,&nrm);
414: if (nrm >= 1.e-7) consistent = PETSC_FALSE;
415: if (flg1) {
416: if (consistent) {
417: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
418: } else {
419: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
420: }
421: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
422: }
423: if (!consistent && flg1) {VecView(r,viewer);}
424: if (!consistent && flg2) {VecView(r,viewer);}
425: VecDestroy(&r);
426: }
428: for (j=0; j<n; j++) {
429: (*mat->ops->mult)(mat,sp->vecs[j],l);
430: VecNorm(l,NORM_2,&nrm);
431: if (nrm >= 1.e-7) consistent = PETSC_FALSE;
432: if (flg1) {
433: if (consistent) {
434: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
435: } else {
436: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
437: consistent = PETSC_FALSE;
438: }
439: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
440: }
441: if (!consistent && flg1) {VecView(l,viewer);}
442: if (!consistent && flg2) {VecView(l,viewer);}
443: }
445: if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
446: VecDestroy(&l);
447: if (isNull) *isNull = consistent;
448: return(0);
449: }