Actual source code: matnull.c
petsc-3.5.4 2015-05-23
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: {
87: PetscErrorCode ierr;
88: const PetscScalar *x;
89: PetscScalar *v[6],dots[5];
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(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);
102: break;
103: case 2:
104: case 3:
105: nmodes = (dim == 2) ? 3 : 6;
106: VecCreate(PetscObjectComm((PetscObject)coords),&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:i-1] */
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],NULL);
152: }
153: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),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_(PetscObjectComm((PetscObject)sp));
188: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
189: if (iascii) {
190: PetscViewerFormat format;
191: PetscInt i;
192: PetscViewerGetFormat(viewer,&format);
193: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
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.18 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);
253: *SP = NULL;
254: MatInitializePackage();
256: PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
258: sp->has_cnst = has_cnst;
259: sp->n = n;
260: sp->vecs = 0;
261: sp->alpha = 0;
262: sp->remove = 0;
263: sp->rmctx = 0;
265: if (n) {
266: PetscMalloc1(n,&sp->vecs);
267: PetscMalloc1(n,&sp->alpha);
268: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
269: for (i=0; i<n; i++) {
270: PetscObjectReference((PetscObject)vecs[i]);
271: sp->vecs[i] = vecs[i];
272: }
273: }
275: *SP = sp;
276: return(0);
277: }
281: /*@
282: MatNullSpaceDestroy - Destroys a data structure used to project vectors
283: out of null spaces.
285: Collective on MatNullSpace
287: Input Parameter:
288: . sp - the null space context to be destroyed
290: Level: advanced
292: .keywords: PC, null space, destroy
294: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
295: @*/
296: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
297: {
301: if (!*sp) return(0);
303: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
305: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
306: PetscFree((*sp)->alpha);
307: PetscHeaderDestroy(sp);
308: return(0);
309: }
313: /*@C
314: MatNullSpaceRemove - Removes all the components of a null space from a vector.
316: Collective on MatNullSpace
318: Input Parameters:
319: + sp - the null space context
320: - vec - the vector from which the null space is to be removed
322: Level: advanced
324: .keywords: PC, null space, remove
326: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
327: @*/
328: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
329: {
330: PetscScalar sum;
331: PetscInt i,N;
338: if (sp->has_cnst) {
339: VecGetSize(vec,&N);
340: if (N > 0) {
341: VecSum(vec,&sum);
342: sum = sum/((PetscScalar)(-1.0*N));
343: VecShift(vec,sum);
344: }
345: }
347: if (sp->n) {
348: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
349: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
350: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
351: }
353: if (sp->remove) {
354: (*sp->remove)(sp,vec,sp->rmctx);
355: }
356: return(0);
357: }
361: /*@
362: MatNullSpaceTest - Tests if the claimed null space is really a
363: null space of a matrix
365: Collective on MatNullSpace
367: Input Parameters:
368: + sp - the null space context
369: - mat - the matrix
371: Output Parameters:
372: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
374: Level: advanced
376: .keywords: PC, null space, remove
378: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
379: @*/
380: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
381: {
382: PetscScalar sum;
383: PetscReal nrm;
384: PetscInt j,n,N;
386: Vec l,r;
387: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
388: PetscViewer viewer;
393: n = sp->n;
394: PetscOptionsGetBool(NULL,"-mat_null_space_test_view",&flg1,NULL);
395: PetscOptionsGetBool(NULL,"-mat_null_space_test_view_draw",&flg2,NULL);
397: if (n) {
398: VecDuplicate(sp->vecs[0],&l);
399: } else {
400: MatGetVecs(mat,&l,NULL);
401: }
403: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
404: if (sp->has_cnst) {
405: VecDuplicate(l,&r);
406: VecGetSize(l,&N);
407: sum = 1.0/N;
408: VecSet(l,sum);
409: MatMult(mat,l,r);
410: VecNorm(r,NORM_2,&nrm);
411: if (nrm >= 1.e-7) consistent = PETSC_FALSE;
412: if (flg1) {
413: if (consistent) {
414: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
415: } else {
416: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
417: }
418: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
419: }
420: if (!consistent && flg1) {VecView(r,viewer);}
421: if (!consistent && flg2) {VecView(r,viewer);}
422: VecDestroy(&r);
423: }
425: for (j=0; j<n; j++) {
426: (*mat->ops->mult)(mat,sp->vecs[j],l);
427: VecNorm(l,NORM_2,&nrm);
428: if (nrm >= 1.e-7) consistent = PETSC_FALSE;
429: if (flg1) {
430: if (consistent) {
431: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
432: } else {
433: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
434: consistent = PETSC_FALSE;
435: }
436: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
437: }
438: if (!consistent && flg1) {VecView(l,viewer);}
439: if (!consistent && flg2) {VecView(l,viewer);}
440: }
442: if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
443: VecDestroy(&l);
444: if (isNull) *isNull = consistent;
445: return(0);
446: }