Actual source code: matnull.c

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

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


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

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

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

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: }

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

218:    Collective on MPI_Comm

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

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

232:    Level: advanced

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

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

240:   Users manual sections:
241: .   Section 4.19 Solving Singular Systems

243: .keywords: PC, null space, create

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

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

291:   *SP = NULL;
292:   MatInitializePackage();

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

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

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

313:   *SP = sp;
314:   return(0);
315: }

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

321:    Collective on MatNullSpace

323:    Input Parameter:
324: .  sp - the null space context to be destroyed

326:    Level: advanced

328: .keywords: PC, null space, destroy

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

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

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

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

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

355:    Collective on MatNullSpace

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

361:    Level: advanced

363: .keywords: PC, null space, remove

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

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

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

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

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

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

403:    Collective on MatNullSpace

405:    Input Parameters:
406: +  sp - the null space context
407: -  mat - the matrix

409:    Output Parameters:
410: .  isNull - PETSC_TRUE if the nullspace is valid for this matrix

412:    Level: advanced

414: .keywords: PC, null space, remove

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

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

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

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

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

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