Actual source code: matnull.c
petsc-3.8.4 2018-03-24
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: .keywords: PC, null space, create
25: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
26: @*/
27: PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace,Vec,void*),void *ctx)
28: {
31: sp->remove = rem;
32: sp->rmctx = ctx;
33: return(0);
34: }
36: /*@C
37: MatNullSpaceGetVecs - get vectors defining the null space
39: Not Collective
41: Input Arguments:
42: . sp - null space object
44: Output Arguments:
45: + has_cnst - PETSC_TRUE if the null space contains the constant vector, otherwise PETSC_FALSE
46: . n - number of vectors (excluding constant vector) in null space
47: - vecs - orthonormal vectors that span the null space (excluding the constant vector)
49: Level: developer
51: Notes:
52: These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller
54: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
55: @*/
56: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
57: {
61: if (has_const) *has_const = sp->has_cnst;
62: if (n) *n = sp->n;
63: if (vecs) *vecs = sp->vecs;
64: return(0);
65: }
67: /*@
68: MatNullSpaceCreateRigidBody - create rigid body modes from coordinates
70: Collective on Vec
72: Input Argument:
73: . coords - block of coordinates of each node, must have block size set
75: Output Argument:
76: . sp - the null space
78: Level: advanced
80: .seealso: MatNullSpaceCreate()
81: @*/
82: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
83: {
84: PetscErrorCode ierr;
85: const PetscScalar *x;
86: PetscScalar *v[6],dots[5];
87: Vec vec[6];
88: PetscInt n,N,dim,nmodes,i,j;
89: 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: {
177: PetscBool iascii;
181: if (!viewer) {
182: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
183: }
187: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
188: if (iascii) {
189: PetscViewerFormat format;
190: PetscInt i;
191: PetscViewerGetFormat(viewer,&format);
192: PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
193: PetscViewerASCIIPushTab(viewer);
194: PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
195: if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
196: if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
197: for (i=0; i<sp->n; i++) {
198: VecView(sp->vecs[i],viewer);
199: }
200: }
201: PetscViewerASCIIPopTab(viewer);
202: }
203: return(0);
204: }
206: /*@C
207: MatNullSpaceCreate - Creates a data structure used to project vectors
208: out of null spaces.
210: Collective on MPI_Comm
212: Input Parameters:
213: + comm - the MPI communicator associated with the object
214: . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
215: . n - number of vectors (excluding constant vector) in null space
216: - vecs - the vectors that span the null space (excluding the constant vector);
217: these vectors must be orthonormal. These vectors are NOT copied, so do not change them
218: after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
219: for them by one).
221: Output Parameter:
222: . SP - the null space context
224: Level: advanced
226: Notes: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.
228: 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
229: need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().
231: Users manual sections:
232: . Section 4.19 Solving Singular Systems
234: .keywords: PC, null space, create
236: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
237: @*/
238: PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
239: {
240: MatNullSpace sp;
242: PetscInt i;
245: if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
249: if (n) {
250: for (i=0; i<n; i++) {
251: /* prevent the user from changes values in the vector */
252: VecLockPush(vecs[i]);
253: }
254: }
255: #if defined(PETSC_USE_DEBUG)
256: if (n) {
257: PetscScalar *dots;
258: for (i=0; i<n; i++) {
259: PetscReal norm;
260: VecNorm(vecs[i],NORM_2,&norm);
261: 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);
262: }
263: if (has_cnst) {
264: for (i=0; i<n; i++) {
265: PetscScalar sum;
266: VecSum(vecs[i],&sum);
267: 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));
268: }
269: }
270: PetscMalloc1(n-1,&dots);
271: for (i=0; i<n-1; i++) {
272: PetscInt j;
273: VecMDot(vecs[i],n-i-1,vecs+i+1,dots);
274: for (j=0;j<n-i-1;j++) {
275: 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]));
276: }
277: }
278: PetscFree(dots);
279: }
280: #endif
282: *SP = NULL;
283: MatInitializePackage();
285: PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);
287: sp->has_cnst = has_cnst;
288: sp->n = n;
289: sp->vecs = 0;
290: sp->alpha = 0;
291: sp->remove = 0;
292: sp->rmctx = 0;
294: if (n) {
295: PetscMalloc1(n,&sp->vecs);
296: PetscMalloc1(n,&sp->alpha);
297: PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
298: for (i=0; i<n; i++) {
299: PetscObjectReference((PetscObject)vecs[i]);
300: sp->vecs[i] = vecs[i];
301: }
302: }
304: *SP = sp;
305: return(0);
306: }
308: /*@
309: MatNullSpaceDestroy - Destroys a data structure used to project vectors
310: out of null spaces.
312: Collective on MatNullSpace
314: Input Parameter:
315: . sp - the null space context to be destroyed
317: Level: advanced
319: .keywords: PC, null space, destroy
321: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
322: @*/
323: PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp)
324: {
326: PetscInt i;
329: if (!*sp) return(0);
331: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
333: for (i=0; i < (*sp)->n; i++) {
334: VecLockPop((*sp)->vecs[i]);
335: }
337: VecDestroyVecs((*sp)->n,&(*sp)->vecs);
338: PetscFree((*sp)->alpha);
339: PetscHeaderDestroy(sp);
340: return(0);
341: }
343: /*@C
344: MatNullSpaceRemove - Removes all the components of a null space from a vector.
346: Collective on MatNullSpace
348: Input Parameters:
349: + sp - the null space context
350: - vec - the vector from which the null space is to be removed
352: Level: advanced
354: .keywords: PC, null space, remove
356: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
357: @*/
358: PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec)
359: {
360: PetscScalar sum;
361: PetscInt i,N;
368: if (sp->has_cnst) {
369: VecGetSize(vec,&N);
370: if (N > 0) {
371: VecSum(vec,&sum);
372: sum = sum/((PetscScalar)(-1.0*N));
373: VecShift(vec,sum);
374: }
375: }
377: if (sp->n) {
378: VecMDot(vec,sp->n,sp->vecs,sp->alpha);
379: for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
380: VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
381: }
383: if (sp->remove) {
384: (*sp->remove)(sp,vec,sp->rmctx);
385: }
386: return(0);
387: }
389: /*@
390: MatNullSpaceTest - Tests if the claimed null space is really a
391: null space of a matrix
393: Collective on MatNullSpace
395: Input Parameters:
396: + sp - the null space context
397: - mat - the matrix
399: Output Parameters:
400: . isNull - PETSC_TRUE if the nullspace is valid for this matrix
402: Level: advanced
404: .keywords: PC, null space, remove
406: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
407: @*/
408: PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull)
409: {
410: PetscScalar sum;
411: PetscReal nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
412: PetscInt j,n,N;
414: Vec l,r;
415: PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
416: PetscViewer viewer;
421: n = sp->n;
422: PetscOptionsGetBool(((PetscObject)sp)->options,NULL,"-mat_null_space_test_view",&flg1,NULL);
423: PetscOptionsGetBool(((PetscObject)sp)->options,NULL,"-mat_null_space_test_view_draw",&flg2,NULL);
425: if (n) {
426: VecDuplicate(sp->vecs[0],&l);
427: } else {
428: MatCreateVecs(mat,&l,NULL);
429: }
431: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
432: if (sp->has_cnst) {
433: VecDuplicate(l,&r);
434: VecGetSize(l,&N);
435: sum = 1.0/N;
436: VecSet(l,sum);
437: MatMult(mat,l,r);
438: VecNorm(r,NORM_2,&nrm);
439: if (nrm >= tol) consistent = PETSC_FALSE;
440: if (flg1) {
441: if (consistent) {
442: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
443: } else {
444: PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
445: }
446: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
447: }
448: if (!consistent && flg1) {VecView(r,viewer);}
449: if (!consistent && flg2) {VecView(r,viewer);}
450: VecDestroy(&r);
451: }
453: for (j=0; j<n; j++) {
454: (*mat->ops->mult)(mat,sp->vecs[j],l);
455: VecNorm(l,NORM_2,&nrm);
456: if (nrm >= tol) consistent = PETSC_FALSE;
457: if (flg1) {
458: if (consistent) {
459: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
460: } else {
461: PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
462: consistent = PETSC_FALSE;
463: }
464: PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
465: }
466: if (!consistent && flg1) {VecView(l,viewer);}
467: if (!consistent && flg2) {VecView(l,viewer);}
468: }
470: if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
471: VecDestroy(&l);
472: if (isNull) *isNull = consistent;
473: return(0);
474: }