Actual source code: matnull.c

petsc-master 2018-01-20
Report Typos and Errors

  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:    Notes:  If you are solving an elasticity problems you should likely use this, in conjunction with ee MatSetNearNullspace(), to provide information that 
 81:            the PCGAMG preconditioner can use to construct a much more efficient preconditioner.

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


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

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

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

165: /*@C
166:    MatNullSpaceView - Visualizes a null space object.

168:    Collective on MatNullSpace

170:    Input Parameters:
171: +  matnull - the null space
172: -  viewer - visualization context

174:    Level: advanced

176:    Fortran Note:
177:    This routine is not supported in Fortran.

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

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

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

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

217:    Collective on MPI_Comm

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

228:    Output Parameter:
229: .  SP - the null space context

231:    Level: advanced

233:    Notes: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.

235:       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
236:        need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().

238:   Users manual sections:
239: .   Section 4.19 Solving Singular Systems

241: .keywords: PC, null space, create

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

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

289:   *SP = NULL;
290:   MatInitializePackage();

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

294:   sp->has_cnst = has_cnst;
295:   sp->n        = n;
296:   sp->vecs     = 0;
297:   sp->alpha    = 0;
298:   sp->remove   = 0;
299:   sp->rmctx    = 0;

301:   if (n) {
302:     PetscMalloc1(n,&sp->vecs);
303:     PetscMalloc1(n,&sp->alpha);
304:     PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
305:     for (i=0; i<n; i++) {
306:       PetscObjectReference((PetscObject)vecs[i]);
307:       sp->vecs[i] = vecs[i];
308:     }
309:   }

311:   *SP = sp;
312:   return(0);
313: }

315: /*@
316:    MatNullSpaceDestroy - Destroys a data structure used to project vectors
317:    out of null spaces.

319:    Collective on MatNullSpace

321:    Input Parameter:
322: .  sp - the null space context to be destroyed

324:    Level: advanced

326: .keywords: PC, null space, destroy

328: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
329: @*/
330: PetscErrorCode  MatNullSpaceDestroy(MatNullSpace *sp)
331: {
333:   PetscInt       i;

336:   if (!*sp) return(0);
338:   if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}

340:   for (i=0; i < (*sp)->n; i++) {
341:     VecLockPop((*sp)->vecs[i]);
342:   }

344:   VecDestroyVecs((*sp)->n,&(*sp)->vecs);
345:   PetscFree((*sp)->alpha);
346:   PetscHeaderDestroy(sp);
347:   return(0);
348: }

350: /*@C
351:    MatNullSpaceRemove - Removes all the components of a null space from a vector.

353:    Collective on MatNullSpace

355:    Input Parameters:
356: +  sp - the null space context (if this is NULL then no null space is removed)
357: -  vec - the vector from which the null space is to be removed

359:    Level: advanced

361: .keywords: PC, null space, remove

363: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
364: @*/
365: PetscErrorCode  MatNullSpaceRemove(MatNullSpace sp,Vec vec)
366: {
367:   PetscScalar    sum;
368:   PetscInt       i,N;

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

376:   if (sp->has_cnst) {
377:     VecGetSize(vec,&N);
378:     if (N > 0) {
379:       VecSum(vec,&sum);
380:       sum  = sum/((PetscScalar)(-1.0*N));
381:       VecShift(vec,sum);
382:     }
383:   }

385:   if (sp->n) {
386:     VecMDot(vec,sp->n,sp->vecs,sp->alpha);
387:     for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
388:     VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
389:   }

391:   if (sp->remove) {
392:     (*sp->remove)(sp,vec,sp->rmctx);
393:   }
394:   return(0);
395: }

397: /*@
398:    MatNullSpaceTest  - Tests if the claimed null space is really a
399:      null space of a matrix

401:    Collective on MatNullSpace

403:    Input Parameters:
404: +  sp - the null space context
405: -  mat - the matrix

407:    Output Parameters:
408: .  isNull - PETSC_TRUE if the nullspace is valid for this matrix

410:    Level: advanced

412: .keywords: PC, null space, remove

414: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
415: @*/
416: PetscErrorCode  MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool  *isNull)
417: {
418:   PetscScalar    sum;
419:   PetscReal      nrm,tol = 10. * PETSC_SQRT_MACHINE_EPSILON;
420:   PetscInt       j,n,N;
422:   Vec            l,r;
423:   PetscBool      flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
424:   PetscViewer    viewer;

429:   n    = sp->n;
430:   PetscOptionsGetBool(((PetscObject)sp)->options,NULL,"-mat_null_space_test_view",&flg1,NULL);
431:   PetscOptionsGetBool(((PetscObject)sp)->options,NULL,"-mat_null_space_test_view_draw",&flg2,NULL);

433:   if (n) {
434:     VecDuplicate(sp->vecs[0],&l);
435:   } else {
436:     MatCreateVecs(mat,&l,NULL);
437:   }

439:   PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
440:   if (sp->has_cnst) {
441:     VecDuplicate(l,&r);
442:     VecGetSize(l,&N);
443:     sum  = 1.0/N;
444:     VecSet(l,sum);
445:     MatMult(mat,l,r);
446:     VecNorm(r,NORM_2,&nrm);
447:     if (nrm >= tol) consistent = PETSC_FALSE;
448:     if (flg1) {
449:       if (consistent) {
450:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
451:       } else {
452:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
453:       }
454:       PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
455:     }
456:     if (!consistent && flg1) {VecView(r,viewer);}
457:     if (!consistent && flg2) {VecView(r,viewer);}
458:     VecDestroy(&r);
459:   }

461:   for (j=0; j<n; j++) {
462:     (*mat->ops->mult)(mat,sp->vecs[j],l);
463:     VecNorm(l,NORM_2,&nrm);
464:     if (nrm >= tol) consistent = PETSC_FALSE;
465:     if (flg1) {
466:       if (consistent) {
467:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
468:       } else {
469:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
470:         consistent = PETSC_FALSE;
471:       }
472:       PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
473:     }
474:     if (!consistent && flg1) {VecView(l,viewer);}
475:     if (!consistent && flg2) {VecView(l,viewer);}
476:   }

478:   if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
479:   VecDestroy(&l);
480:   if (isNull) *isNull = consistent;
481:   return(0);
482: }