Actual source code: matnull.c
petsc-3.7.5 2017-01-01
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;
95: PetscReal sN;
98: VecGetBlockSize(coords,&dim);
99: VecGetLocalSize(coords,&n);
100: VecGetSize(coords,&N);
101: n /= dim;
102: N /= dim;
103: sN = 1./PetscSqrtReal((PetscReal)N);
104: switch (dim) {
105: case 1:
106: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);
107: break;
108: case 2:
109: case 3:
110: nmodes = (dim == 2) ? 3 : 6;
111: VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);
112: VecSetSizes(vec[0],dim*n,dim*N);
113: VecSetBlockSize(vec[0],dim);
114: VecSetUp(vec[0]);
115: for (i=1; i<nmodes; i++) {VecDuplicate(vec[0],&vec[i]);}
116: for (i=0; i<nmodes; i++) {VecGetArray(vec[i],&v[i]);}
117: VecGetArrayRead(coords,&x);
118: for (i=0; i<n; i++) {
119: if (dim == 2) {
120: v[0][i*2+0] = sN;
121: v[0][i*2+1] = 0.;
122: v[1][i*2+0] = 0.;
123: v[1][i*2+1] = sN;
124: /* Rotations */
125: v[2][i*2+0] = -x[i*2+1];
126: v[2][i*2+1] = x[i*2+0];
127: } else {
128: v[0][i*3+0] = sN;
129: v[0][i*3+1] = 0.;
130: v[0][i*3+2] = 0.;
131: v[1][i*3+0] = 0.;
132: v[1][i*3+1] = sN;
133: v[1][i*3+2] = 0.;
134: v[2][i*3+0] = 0.;
135: v[2][i*3+1] = 0.;
136: v[2][i*3+2] = sN;
138: v[3][i*3+0] = x[i*3+1];
139: v[3][i*3+1] = -x[i*3+0];
140: v[3][i*3+2] = 0.;
141: v[4][i*3+0] = 0.;
142: v[4][i*3+1] = -x[i*3+2];
143: v[4][i*3+2] = x[i*3+1];
144: v[5][i*3+0] = x[i*3+2];
145: v[5][i*3+1] = 0.;
146: v[5][i*3+2] = -x[i*3+0];
147: }
148: }
149: for (i=0; i<nmodes; i++) {VecRestoreArray(vec[i],&v[i]);}
150: VecRestoreArrayRead(coords,&x);
151: for (i=dim; i<nmodes; i++) {
152: /* Orthonormalize vec[i] against vec[0:i-1] */
153: VecMDot(vec[i],i,vec,dots);
154: for (j=0; j<i; j++) dots[j] *= -1.;
155: VecMAXPY(vec[i],i,dots,vec);
156: VecNormalize(vec[i],NULL);
157: }
158: MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);
159: for (i=0; i<nmodes; i++) {VecDestroy(&vec[i]);}
160: }
161: return(0);
162: }
166: /*@C
167: MatNullSpaceView - Visualizes a null space object.
169: Collective on MatNullSpace
171: Input Parameters:
172: + matnull - the null space
173: - viewer - visualization context
175: Level: advanced
177: Fortran Note:
178: This routine is not supported in Fortran.
180: .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen()
181: @*/
182: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
183: {
185: PetscBool iascii;
189: if (!viewer) {
190: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
191: }
195: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
196: if (iascii) {
197: PetscViewerFormat format;
198: PetscInt i;
199: PetscViewerGetFormat(viewer,&format);
200: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
201: PetscViewerASCIIPushTab(viewer);
202: PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
203: if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
204: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
205: for (i=0; i<sp->n; i++) {
206: VecView(sp->vecs[i],viewer);
207: }
208: }
209: PetscViewerASCIIPopTab(viewer);
210: }
211: return(0);
212: }
216: /*@
217: MatNullSpaceCreate - Creates a data structure used to project vectors
218: out of null spaces.
220: Collective on MPI_Comm
222: Input Parameters:
223: + comm - the MPI communicator associated with the object
224: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
225: . n - number of vectors (excluding constant vector) in null space
226: - vecs - the vectors that span the null space (excluding the constant vector);
227: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
228: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
229: for them by one).
231: Output Parameter:
232: . SP - the null space context
234: Level: advanced
236: Notes: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.
238: 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
239: need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().
241: Users manual sections:
242: . Section 4.19 Solving Singular Systems
244: .keywords: PC, null space, create
246: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
247: @*/
248: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
249: {
250: MatNullSpace sp;
252: PetscInt i;
255: if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
259: #if defined(PETSC_USE_DEBUG)
260: if (n) {
261: PetscScalar *dots;
262: for (i=0; i<n; i++) {
263: PetscReal norm;
264: VecNorm(vecs[i],NORM_2,&norm);
265: if (PetscAbsReal(norm - 1.0) > 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);
266: }
267: if (has_cnst) {
268: for (i=0; i<n; i++) {
269: PetscScalar sum;
270: VecSum(vecs[i],&sum);
271: 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));
272: }
273: }
274: PetscMalloc1(n-1,&dots);
275: for (i=0; i<n-1; i++) {
276: PetscInt j;
277: VecMDot(vecs[i],n-i-1,vecs+i+1,dots);
278: for (j=0;j<n-i-1;j++) {
279: 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]));
280: }
281: }
282: PetscFree(dots);
283: }
284: #endif
286: *SP = NULL;
287: MatInitializePackage();
289: PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
291: sp->has_cnst = has_cnst;
292: sp->n = n;
293: sp->vecs = 0;
294: sp->alpha = 0;
295: sp->remove = 0;
296: sp->rmctx = 0;
298: if (n) {
299: PetscMalloc1(n,&sp->vecs);
300: PetscMalloc1(n,&sp->alpha);
301: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
302: for (i=0; i<n; i++) {
303: PetscObjectReference((PetscObject)vecs[i]);
304: sp->vecs[i] = vecs[i];
305: }
306: }
308: *SP = sp;
309: return(0);
310: }
314: /*@
315: MatNullSpaceDestroy - Destroys a data structure used to project vectors
316: out of null spaces.
318: Collective on MatNullSpace
320: Input Parameter:
321: . sp - the null space context to be destroyed
323: Level: advanced
325: .keywords: PC, null space, destroy
327: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
328: @*/
329: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
330: {
334: if (!*sp) return(0);
336: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
338: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
339: PetscFree((*sp)->alpha);
340: PetscHeaderDestroy(sp);
341: return(0);
342: }
346: /*@C
347: MatNullSpaceRemove - Removes all the components of a null space from a vector.
349: Collective on MatNullSpace
351: Input Parameters:
352: + sp - the null space context
353: - vec - the vector from which the null space is to be removed
355: Level: advanced
357: .keywords: PC, null space, remove
359: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
360: @*/
361: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
362: {
363: PetscScalar sum;
364: PetscInt i,N;
371: if (sp->has_cnst) {
372: VecGetSize(vec,&N);
373: if (N > 0) {
374: VecSum(vec,&sum);
375: sum = sum/((PetscScalar)(-1.0*N));
376: VecShift(vec,sum);
377: }
378: }
380: if (sp->n) {
381: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
382: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
383: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
384: }
386: if (sp->remove) {
387: (*sp->remove)(sp,vec,sp->rmctx);
388: }
389: return(0);
390: }
394: /*@
395: MatNullSpaceTest - Tests if the claimed null space is really a
396: null space of a matrix
398: Collective on MatNullSpace
400: Input Parameters:
401: + sp - the null space context
402: - mat - the matrix
404: Output Parameters:
405: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
407: Level: advanced
409: .keywords: PC, null space, remove
411: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
412: @*/
413: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
414: {
415: PetscScalar sum;
416: PetscReal nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
417: PetscInt j,n,N;
419: Vec l,r;
420: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
421: PetscViewer viewer;
426: n = sp->n;
427: PetscOptionsGetBool(((PetscObject)sp)->options,NULL,"-mat_null_space_test_view",&flg1,NULL);
428: PetscOptionsGetBool(((PetscObject)sp)->options,NULL,"-mat_null_space_test_view_draw",&flg2,NULL);
430: if (n) {
431: VecDuplicate(sp->vecs[0],&l);
432: } else {
433: MatCreateVecs(mat,&l,NULL);
434: }
436: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
437: if (sp->has_cnst) {
438: VecDuplicate(l,&r);
439: VecGetSize(l,&N);
440: sum = 1.0/N;
441: VecSet(l,sum);
442: MatMult(mat,l,r);
443: VecNorm(r,NORM_2,&nrm);
444: if (nrm >= tol) consistent = PETSC_FALSE;
445: if (flg1) {
446: if (consistent) {
447: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
448: } else {
449: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
450: }
451: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
452: }
453: if (!consistent && flg1) {VecView(r,viewer);}
454: if (!consistent && flg2) {VecView(r,viewer);}
455: VecDestroy(&r);
456: }
458: for (j=0; j<n; j++) {
459: (*mat->ops->mult)(mat,sp->vecs[j],l);
460: VecNorm(l,NORM_2,&nrm);
461: if (nrm >= tol) consistent = PETSC_FALSE;
462: if (flg1) {
463: if (consistent) {
464: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
465: } else {
466: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
467: consistent = PETSC_FALSE;
468: }
469: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
470: }
471: if (!consistent && flg1) {VecView(l,viewer);}
472: if (!consistent && flg2) {VecView(l,viewer);}
473: }
475: if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
476: VecDestroy(&l);
477: if (isNull) *isNull = consistent;
478: return(0);
479: }