Actual source code: matnull.c


  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: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
 24: @*/
 25: PetscErrorCode  MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace,Vec,void*),void *ctx)
 26: {
 29:   sp->remove = rem;
 30:   sp->rmctx  = ctx;
 31:   return(0);
 32: }

 34: /*@C
 35:    MatNullSpaceGetVecs - get vectors defining the null space

 37:    Not Collective

 39:    Input Parameter:
 40: .  sp - null space object

 42:    Output Parameters:
 43: +  has_cnst - PETSC_TRUE if the null space contains the constant vector, otherwise PETSC_FALSE
 44: .  n - number of vectors (excluding constant vector) in null space
 45: -  vecs - orthonormal vectors that span the null space (excluding the constant vector)

 47:    Level: developer

 49:    Notes:
 50:       These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller

 52: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
 53: @*/
 54: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
 55: {

 59:   if (has_const) *has_const = sp->has_cnst;
 60:   if (n) *n = sp->n;
 61:   if (vecs) *vecs = sp->vecs;
 62:   return(0);
 63: }

 65: /*@
 66:    MatNullSpaceCreateRigidBody - create rigid body modes from coordinates

 68:    Collective on Vec

 70:    Input Parameter:
 71: .  coords - block of coordinates of each node, must have block size set

 73:    Output Parameter:
 74: .  sp - the null space

 76:    Level: advanced

 78:    Notes:
 79:      If you are solving an elasticity problem you should likely use this, in conjunction with MatSetNearNullspace(), to provide information that
 80:      the PCGAMG preconditioner can use to construct a much more efficient preconditioner.

 82:      If you are solving an elasticity problem with pure Neumann boundary conditions you can use this in conjunction with MatSetNullspace() to
 83:      provide this information to the linear solver so it can handle the null space appropriately in the linear solution.

 85: .seealso: MatNullSpaceCreate(), MatSetNearNullspace(), MatSetNullspace()
 86: @*/
 87: PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp)
 88: {
 89:   PetscErrorCode    ierr;
 90:   const PetscScalar *x;
 91:   PetscScalar       *v[6],dots[5];
 92:   Vec               vec[6];
 93:   PetscInt          n,N,dim,nmodes,i,j;
 94:   PetscReal         sN;

 97:   VecGetBlockSize(coords,&dim);
 98:   VecGetLocalSize(coords,&n);
 99:   VecGetSize(coords,&N);
100:   n   /= dim;
101:   N   /= dim;
102:   sN = 1./PetscSqrtReal((PetscReal)N);
103:   switch (dim) {
104:   case 1:
105:     MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);
106:     break;
107:   case 2:
108:   case 3:
109:     nmodes = (dim == 2) ? 3 : 6;
110:     VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);
111:     VecSetSizes(vec[0],dim*n,dim*N);
112:     VecSetBlockSize(vec[0],dim);
113:     VecSetUp(vec[0]);
114:     for (i=1; i<nmodes; i++) {VecDuplicate(vec[0],&vec[i]);}
115:     for (i=0; i<nmodes; i++) {VecGetArray(vec[i],&v[i]);}
116:     VecGetArrayRead(coords,&x);
117:     for (i=0; i<n; i++) {
118:       if (dim == 2) {
119:         v[0][i*2+0] = sN;
120:         v[0][i*2+1] = 0.;
121:         v[1][i*2+0] = 0.;
122:         v[1][i*2+1] = sN;
123:         /* Rotations */
124:         v[2][i*2+0] = -x[i*2+1];
125:         v[2][i*2+1] = x[i*2+0];
126:       } else {
127:         v[0][i*3+0] = sN;
128:         v[0][i*3+1] = 0.;
129:         v[0][i*3+2] = 0.;
130:         v[1][i*3+0] = 0.;
131:         v[1][i*3+1] = sN;
132:         v[1][i*3+2] = 0.;
133:         v[2][i*3+0] = 0.;
134:         v[2][i*3+1] = 0.;
135:         v[2][i*3+2] = sN;

137:         v[3][i*3+0] = x[i*3+1];
138:         v[3][i*3+1] = -x[i*3+0];
139:         v[3][i*3+2] = 0.;
140:         v[4][i*3+0] = 0.;
141:         v[4][i*3+1] = -x[i*3+2];
142:         v[4][i*3+2] = x[i*3+1];
143:         v[5][i*3+0] = x[i*3+2];
144:         v[5][i*3+1] = 0.;
145:         v[5][i*3+2] = -x[i*3+0];
146:       }
147:     }
148:     for (i=0; i<nmodes; i++) {VecRestoreArray(vec[i],&v[i]);}
149:     VecRestoreArrayRead(coords,&x);
150:     for (i=dim; i<nmodes; i++) {
151:       /* Orthonormalize vec[i] against vec[0:i-1] */
152:       VecMDot(vec[i],i,vec,dots);
153:       for (j=0; j<i; j++) dots[j] *= -1.;
154:       VecMAXPY(vec[i],i,dots,vec);
155:       VecNormalize(vec[i],NULL);
156:     }
157:     MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);
158:     for (i=0; i<nmodes; i++) {VecDestroy(&vec[i]);}
159:   }
160:   return(0);
161: }

163: /*@C
164:    MatNullSpaceView - Visualizes a null space object.

166:    Collective on MatNullSpace

168:    Input Parameters:
169: +  matnull - the null space
170: -  viewer - visualization context

172:    Level: advanced

174:    Fortran Note:
175:    This routine is not supported in Fortran.

177: .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen()
178: @*/
179: PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer)
180: {
182:   PetscBool      iascii;

186:   if (!viewer) {
187:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
188:   }

192:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
193:   if (iascii) {
194:     PetscViewerFormat format;
195:     PetscInt          i;
196:     PetscViewerGetFormat(viewer,&format);
197:     PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);
198:     PetscViewerASCIIPushTab(viewer);
199:     PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");
200:     if (sp->remove) {PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");}
201:     if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) {
202:       for (i=0; i<sp->n; i++) {
203:         VecView(sp->vecs[i],viewer);
204:       }
205:     }
206:     PetscViewerASCIIPopTab(viewer);
207:   }
208:   return(0);
209: }

211: /*@C
212:    MatNullSpaceCreate - Creates a data structure used to project vectors
213:    out of null spaces.

215:    Collective

217:    Input Parameters:
218: +  comm - the MPI communicator associated with the object
219: .  has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
220: .  n - number of vectors (excluding constant vector) in null space
221: -  vecs - the vectors that span the null space (excluding the constant vector);
222:           these vectors must be orthonormal. These vectors are NOT copied, so do not change them
223:           after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
224:           for them by one).

226:    Output Parameter:
227: .  SP - the null space context

229:    Level: advanced

231:    Notes:
232:     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: .   sec_singular

240: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
241: @*/
242: PetscErrorCode  MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
243: {
244:   MatNullSpace   sp;
246:   PetscInt       i;

249:   if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
253:   if (n) {
254:     for (i=0; i<n; i++) {
255:       /* prevent the user from changes values in the vector */
256:       VecLockReadPush(vecs[i]);
257:     }
258:   }
259:   if (PetscUnlikelyDebug(n)) {
260:     PetscScalar *dots;
261:     for (i=0; i<n; i++) {
262:       PetscReal norm;
263:       VecNorm(vecs[i],NORM_2,&norm);
264:       if (PetscAbsReal(norm - 1) > 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);
265:     }
266:     if (has_cnst) {
267:       for (i=0; i<n; i++) {
268:         PetscScalar sum;
269:         VecSum(vecs[i],&sum);
270:         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));
271:       }
272:     }
273:     PetscMalloc1(n-1,&dots);
274:     for (i=0; i<n-1; i++) {
275:       PetscInt j;
276:       VecMDot(vecs[i],n-i-1,vecs+i+1,dots);
277:       for (j=0;j<n-i-1;j++) {
278:         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]));
279:       }
280:     }
281:     PetscFree(dots);
282:   }

284:   *SP = NULL;
285:   MatInitializePackage();

287:   PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);

289:   sp->has_cnst = has_cnst;
290:   sp->n        = n;
291:   sp->vecs     = NULL;
292:   sp->alpha    = NULL;
293:   sp->remove   = NULL;
294:   sp->rmctx    = NULL;

296:   if (n) {
297:     PetscMalloc1(n,&sp->vecs);
298:     PetscMalloc1(n,&sp->alpha);
299:     PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
300:     for (i=0; i<n; i++) {
301:       PetscObjectReference((PetscObject)vecs[i]);
302:       sp->vecs[i] = vecs[i];
303:     }
304:   }

306:   *SP = sp;
307:   return(0);
308: }

310: /*@
311:    MatNullSpaceDestroy - Destroys a data structure used to project vectors
312:    out of null spaces.

314:    Collective on MatNullSpace

316:    Input Parameter:
317: .  sp - the null space context to be destroyed

319:    Level: advanced

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 = NULL; return(0);}

333:   for (i=0; i < (*sp)->n; i++) {
334:     VecLockReadPop((*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 (if this is NULL then no null space is removed)
350: -  vec - the vector from which the null space is to be removed

352:    Level: advanced

354: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
355: @*/
356: PetscErrorCode  MatNullSpaceRemove(MatNullSpace sp,Vec vec)
357: {
358:   PetscScalar    sum;
359:   PetscInt       i,N;

363:   if (!sp) return(0);

367:   if (sp->has_cnst) {
368:     VecGetSize(vec,&N);
369:     if (N > 0) {
370:       VecSum(vec,&sum);
371:       sum  = sum/((PetscScalar)(-1.0*N));
372:       VecShift(vec,sum);
373:     }
374:   }

376:   if (sp->n) {
377:     VecMDot(vec,sp->n,sp->vecs,sp->alpha);
378:     for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
379:     VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
380:   }

382:   if (sp->remove) {
383:     (*sp->remove)(sp,vec,sp->rmctx);
384:   }
385:   return(0);
386: }

388: /*@
389:    MatNullSpaceTest  - Tests if the claimed null space is really a
390:      null space of a matrix

392:    Collective on MatNullSpace

394:    Input Parameters:
395: +  sp - the null space context
396: -  mat - the matrix

398:    Output Parameters:
399: .  isNull - PETSC_TRUE if the nullspace is valid for this matrix

401:    Level: advanced

403: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
404: @*/
405: PetscErrorCode  MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool  *isNull)
406: {
407:   PetscScalar    sum;
408:   PetscReal      nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
409:   PetscInt       j,n,N;
411:   Vec            l,r;
412:   PetscBool      flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
413:   PetscViewer    viewer;

418:   n    = sp->n;
419:   PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view",&flg1,NULL);
420:   PetscOptionsGetBool(((PetscObject)sp)->options,((PetscObject)mat)->prefix,"-mat_null_space_test_view_draw",&flg2,NULL);

422:   if (n) {
423:     VecDuplicate(sp->vecs[0],&l);
424:   } else {
425:     MatCreateVecs(mat,&l,NULL);
426:   }

428:   PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
429:   if (sp->has_cnst) {
430:     VecDuplicate(l,&r);
431:     VecGetSize(l,&N);
432:     sum  = 1.0/PetscSqrtReal(N);
433:     VecSet(l,sum);
434:     MatMult(mat,l,r);
435:     VecNorm(r,NORM_2,&nrm);
436:     if (nrm >= tol) consistent = PETSC_FALSE;
437:     if (flg1) {
438:       if (consistent) {
439:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
440:       } else {
441:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
442:       }
443:       PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
444:     }
445:     if (!consistent && flg1) {VecView(r,viewer);}
446:     if (!consistent && flg2) {VecView(r,viewer);}
447:     VecDestroy(&r);
448:   }

450:   for (j=0; j<n; j++) {
451:     (*mat->ops->mult)(mat,sp->vecs[j],l);
452:     VecNorm(l,NORM_2,&nrm);
453:     if (nrm >= tol) consistent = PETSC_FALSE;
454:     if (flg1) {
455:       if (consistent) {
456:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
457:       } else {
458:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
459:         consistent = PETSC_FALSE;
460:       }
461:       PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
462:     }
463:     if (!consistent && flg1) {VecView(l,viewer);}
464:     if (!consistent && flg2) {VecView(l,viewer);}
465:   }

467:   if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
468:   VecDestroy(&l);
469:   if (isNull) *isNull = consistent;
470:   return(0);
471: }