Actual source code: matnull.c
petsc-3.3-p7 2013-05-11
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(), KSPSetNullSpace(), 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: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
56: @*/
57: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
58: {
62: if (has_const) *has_const = sp->has_cnst;
63: if (n) *n = sp->n;
64: if (vecs) *vecs = sp->vecs;
65: return(0);
66: }
70: /*@
71: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
73: Collective on Vec
75: Input Argument:
76: . coords - block of coordinates of each node, must have block size set
78: Output Argument:
79: . sp - the null space
81: Level: advanced
83: .seealso: MatNullSpaceCreate()
84: @*/
85: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
86: {
88: const PetscScalar *x;
89: PetscScalar *v[6],dots[3];
90: Vec vec[6];
91: PetscInt n,N,dim,nmodes,i,j;
94: VecGetBlockSize(coords,&dim);
95: VecGetLocalSize(coords,&n);
96: VecGetSize(coords,&N);
97: n /= dim;
98: N /= dim;
99: switch (dim) {
100: case 1:
101: MatNullSpaceCreate(((PetscObject)coords)->comm,PETSC_TRUE,0,PETSC_NULL,sp);
102: break;
103: case 2:
104: case 3:
105: nmodes = (dim == 2) ? 3 : 6;
106: VecCreate(((PetscObject)coords)->comm,&vec[0]);
107: VecSetSizes(vec[0],dim*n,dim*N);
108: VecSetBlockSize(vec[0],dim);
109: VecSetUp(vec[0]);
110: for (i=1; i<nmodes; i++) {VecDuplicate(vec[0],&vec[i]);}
111: for (i=0; i<nmodes; i++) {VecGetArray(vec[i],&v[i]);}
112: VecGetArrayRead(coords,&x);
113: for (i=0; i<n; i++) {
114: if (dim == 2) {
115: v[0][i*2+0] = 1./N;
116: v[0][i*2+1] = 0.;
117: v[1][i*2+0] = 0.;
118: v[1][i*2+1] = 1./N;
119: /* Rotations */
120: v[2][i*2+0] = -x[i*2+1];
121: v[2][i*2+1] = x[i*2+0];
122: } else {
123: v[0][i*3+0] = 1./N;
124: v[0][i*3+1] = 0.;
125: v[0][i*3+2] = 0.;
126: v[1][i*3+0] = 0.;
127: v[1][i*3+1] = 1./N;
128: v[1][i*3+2] = 0.;
129: v[2][i*3+0] = 0.;
130: v[2][i*3+1] = 0.;
131: v[2][i*3+2] = 1./N;
133: v[3][i*3+0] = x[i*3+1];
134: v[3][i*3+1] = -x[i*3+0];
135: v[3][i*3+2] = 0.;
136: v[4][i*3+0] = 0.;
137: v[4][i*3+1] = -x[i*3+2];
138: v[4][i*3+2] = x[i*3+1];
139: v[5][i*3+0] = x[i*3+2];
140: v[5][i*3+1] = 0.;
141: v[5][i*3+2] = -x[i*3+0];
142: }
143: }
144: for (i=0; i<nmodes; i++) {VecRestoreArray(vec[i],&v[i]);}
145: VecRestoreArrayRead(coords,&x);
146: for (i=dim; i<nmodes; i++) {
147: /* Orthonormalize vec[i] against vec[0:dim] */
148: VecMDot(vec[i],i,vec,dots);
149: for (j=0; j<i; j++) dots[j] *= -1.;
150: VecMAXPY(vec[i],i,dots,vec);
151: VecNormalize(vec[i],PETSC_NULL);
152: }
153: MatNullSpaceCreate(((PetscObject)coords)->comm,PETSC_FALSE,nmodes,vec,sp);
154: for (i=0; i<nmodes; i++) {VecDestroy(&vec[i]);}
155: }
156: return(0);
157: }
161: /*@C
162: MatNullSpaceView - Visualizes a null space object.
164: Collective on MatNullSpace
166: Input Parameters:
167: + matnull - the null space
168: - viewer - visualization context
170: Level: advanced
172: Fortran Note:
173: This routine is not supported in Fortran.
175: .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen()
176: @*/
177: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
178: {
180: PetscBool iascii;
184: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(((PetscObject)sp)->comm);
188: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
189: if (iascii) {
190: PetscViewerFormat format;
191: PetscInt i;
192: PetscViewerGetFormat(viewer,&format);
193: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer,"MatNullSpace Object");
194: PetscViewerASCIIPushTab(viewer);
195: PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1?"":"s",sp->has_cnst?" and the constant":"");
196: if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
197: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
198: for (i=0; i<sp->n; i++) {
199: VecView(sp->vecs[i],viewer);
200: }
201: }
202: PetscViewerASCIIPopTab(viewer);
203: }
204: return(0);
205: }
209: /*@
210: MatNullSpaceCreate - Creates a data structure used to project vectors
211: out of null spaces.
213: Collective on MPI_Comm
215: Input Parameters:
216: + comm - the MPI communicator associated with the object
217: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
218: . n - number of vectors (excluding constant vector) in null space
219: - vecs - the vectors that span the null space (excluding the constant vector);
220: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
221: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
222: for them by one).
224: Output Parameter:
225: . SP - the null space context
227: Level: advanced
229: Notes: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.
231: 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
232: need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().
234: Users manual sections:
235: . Section 4.16 Solving Singular Systems
237: .keywords: PC, null space, create
239: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
240: @*/
241: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
242: {
243: MatNullSpace sp;
245: PetscInt i;
248: if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
252:
253: *SP = PETSC_NULL;
254: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
255: MatInitializePackage(PETSC_NULL);
256: #endif
258: PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_CLASSID,0,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
260: sp->has_cnst = has_cnst;
261: sp->n = n;
262: sp->vecs = 0;
263: sp->alpha = 0;
264: sp->vec = 0;
265: sp->remove = 0;
266: sp->rmctx = 0;
268: if (n) {
269: PetscMalloc(n*sizeof(Vec),&sp->vecs);
270: PetscMalloc(n*sizeof(PetscScalar),&sp->alpha);
271: PetscLogObjectMemory(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: VecDestroy(&(*sp)->vec);
309: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
310: PetscFree((*sp)->alpha);
311: PetscHeaderDestroy(sp);
312: return(0);
313: }
317: /*@C
318: MatNullSpaceRemove - Removes all the components of a null space from a vector.
320: Collective on MatNullSpace
322: Input Parameters:
323: + sp - the null space context
324: . vec - the vector from which the null space is to be removed
325: - out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise
326: the removal is done in-place (in vec)
328: Note: The user is not responsible for the vector returned and should not destroy it.
330: Level: advanced
332: .keywords: PC, null space, remove
334: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
335: @*/
336: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out)
337: {
338: PetscScalar sum;
339: PetscInt i,N;
346: if (out) {
348: if (!sp->vec) {
349: VecDuplicate(vec,&sp->vec);
350: PetscLogObjectParent(sp,sp->vec);
351: }
352: VecCopy(vec,sp->vec);
353: vec = *out = sp->vec;
354: }
355:
356: if (sp->has_cnst) {
357: VecGetSize(vec,&N);
358: if (N > 0) {
359: VecSum(vec,&sum);
360: sum = sum/((PetscScalar)(-1.0*N));
361: VecShift(vec,sum);
362: }
363: }
364:
365: if (sp->n) {
366: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
367: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
368: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
369: }
371: if (sp->remove){
372: (*sp->remove)(sp,vec,sp->rmctx);
373: }
374: return(0);
375: }
379: /*@
380: MatNullSpaceTest - Tests if the claimed null space is really a
381: null space of a matrix
383: Collective on MatNullSpace
385: Input Parameters:
386: + sp - the null space context
387: - mat - the matrix
389: Output Parameters:
390: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
392: Level: advanced
394: .keywords: PC, null space, remove
396: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
397: @*/
398: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
399: {
400: PetscScalar sum;
401: PetscReal nrm;
402: PetscInt j,n,N;
404: Vec l,r;
405: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
406: PetscViewer viewer;
411: n = sp->n;
412: PetscOptionsGetBool(PETSC_NULL,"-mat_null_space_test_view",&flg1,PETSC_NULL);
413: PetscOptionsGetBool(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2,PETSC_NULL);
415: if (!sp->vec) {
416: if (n) {
417: VecDuplicate(sp->vecs[0],&sp->vec);
418: } else {
419: MatGetVecs(mat,&sp->vec,PETSC_NULL);
420: }
421: }
422: l = sp->vec;
424: PetscViewerASCIIGetStdout(((PetscObject)sp)->comm,&viewer);
425: if (sp->has_cnst) {
426: VecDuplicate(l,&r);
427: VecGetSize(l,&N);
428: sum = 1.0/N;
429: VecSet(l,sum);
430: MatMult(mat,l,r);
431: VecNorm(r,NORM_2,&nrm);
432: if (nrm >= 1.e-7) consistent = PETSC_FALSE;
433: if (flg1) {
434: if (consistent) {
435: PetscPrintf(((PetscObject)sp)->comm,"Constants are likely null vector");
436: } else {
437: PetscPrintf(((PetscObject)sp)->comm,"Constants are unlikely null vector ");
438: }
439: PetscPrintf(((PetscObject)sp)->comm,"|| A * 1/N || = %G\n",nrm);
440: }
441: if (!consistent && flg1) {VecView(r,viewer);}
442: if (!consistent && flg2) {VecView(r,viewer);}
443: VecDestroy(&r);
444: }
446: for (j=0; j<n; j++) {
447: (*mat->ops->mult)(mat,sp->vecs[j],l);
448: VecNorm(l,NORM_2,&nrm);
449: if (nrm >= 1.e-7) consistent = PETSC_FALSE;
450: if (flg1) {
451: if (consistent) {
452: PetscPrintf(((PetscObject)sp)->comm,"Null vector %D is likely null vector",j);
453: } else {
454: PetscPrintf(((PetscObject)sp)->comm,"Null vector %D unlikely null vector ",j);
455: consistent = PETSC_FALSE;
456: }
457: PetscPrintf(((PetscObject)sp)->comm,"|| A * v[%D] || = %G\n",j,nrm);
458: }
459: if (!consistent && flg1) {VecView(l,viewer);}
460: if (!consistent && flg2) {VecView(l,viewer);}
461: }
463: if (sp->remove) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
464: if (isNull) *isNull = consistent;
465: return(0);
466: }