Actual source code: matnull.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  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(), KSPSetNullSpace(), 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: .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace()
 56: @*/
 57: PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs)
 58: {

 62:   if (has_const) *has_const = sp->has_cnst;
 63:   if (n) *n = sp->n;
 64:   if (vecs) *vecs = sp->vecs;
 65:   return(0);
 66: }

 70: /*@
 71:    MatNullSpaceCreateRigidBody - create rigid body modes from coordinates

 73:    Collective on Vec

 75:    Input Argument:
 76: .  coords - block of coordinates of each node, must have block size set

 78:    Output Argument:
 79: .  sp - the null space

 81:    Level: advanced

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

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

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

161: /*@C
162:    MatNullSpaceView - Visualizes a null space object.

164:    Collective on MatNullSpace

166:    Input Parameters:
167: +  matnull - the null space
168: -  viewer - visualization context

170:    Level: advanced

172:    Fortran Note:
173:    This routine is not supported in Fortran.

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

184:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)sp));

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

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

213:    Collective on MPI_Comm

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

224:    Output Parameter:
225: .  SP - the null space context

227:    Level: advanced

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

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

234:   Users manual sections:
235: .   Section 4.18 Solving Singular Systems

237: .keywords: PC, null space, create

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

248:   if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);

253:   *SP = NULL;
254:   MatInitializePackage();

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

258:   sp->has_cnst = has_cnst;
259:   sp->n        = n;
260:   sp->vecs     = 0;
261:   sp->alpha    = 0;
262:   sp->remove   = 0;
263:   sp->rmctx    = 0;

265:   if (n) {
266:     PetscMalloc1(n,&sp->vecs);
267:     PetscMalloc1(n,&sp->alpha);
268:     PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));
269:     for (i=0; i<n; i++) {
270:       PetscObjectReference((PetscObject)vecs[i]);
271:       sp->vecs[i] = vecs[i];
272:     }
273:   }

275:   *SP = sp;
276:   return(0);
277: }

281: /*@
282:    MatNullSpaceDestroy - Destroys a data structure used to project vectors
283:    out of null spaces.

285:    Collective on MatNullSpace

287:    Input Parameter:
288: .  sp - the null space context to be destroyed

290:    Level: advanced

292: .keywords: PC, null space, destroy

294: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
295: @*/
296: PetscErrorCode  MatNullSpaceDestroy(MatNullSpace *sp)
297: {

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

305:   VecDestroyVecs((*sp)->n,&(*sp)->vecs);
306:   PetscFree((*sp)->alpha);
307:   PetscHeaderDestroy(sp);
308:   return(0);
309: }

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

316:    Collective on MatNullSpace

318:    Input Parameters:
319: +  sp - the null space context
320: -  vec - the vector from which the null space is to be removed

322:    Level: advanced

324: .keywords: PC, null space, remove

326: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
327: @*/
328: PetscErrorCode  MatNullSpaceRemove(MatNullSpace sp,Vec vec)
329: {
330:   PetscScalar    sum;
331:   PetscInt       i,N;


338:   if (sp->has_cnst) {
339:     VecGetSize(vec,&N);
340:     if (N > 0) {
341:       VecSum(vec,&sum);
342:       sum  = sum/((PetscScalar)(-1.0*N));
343:       VecShift(vec,sum);
344:     }
345:   }

347:   if (sp->n) {
348:     VecMDot(vec,sp->n,sp->vecs,sp->alpha);
349:     for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
350:     VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);
351:   }

353:   if (sp->remove) {
354:     (*sp->remove)(sp,vec,sp->rmctx);
355:   }
356:   return(0);
357: }

361: /*@
362:    MatNullSpaceTest  - Tests if the claimed null space is really a
363:      null space of a matrix

365:    Collective on MatNullSpace

367:    Input Parameters:
368: +  sp - the null space context
369: -  mat - the matrix

371:    Output Parameters:
372: .  isNull - PETSC_TRUE if the nullspace is valid for this matrix

374:    Level: advanced

376: .keywords: PC, null space, remove

378: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
379: @*/
380: PetscErrorCode  MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool  *isNull)
381: {
382:   PetscScalar    sum;
383:   PetscReal      nrm;
384:   PetscInt       j,n,N;
386:   Vec            l,r;
387:   PetscBool      flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE;
388:   PetscViewer    viewer;

393:   n    = sp->n;
394:   PetscOptionsGetBool(NULL,"-mat_null_space_test_view",&flg1,NULL);
395:   PetscOptionsGetBool(NULL,"-mat_null_space_test_view_draw",&flg2,NULL);

397:   if (n) {
398:     VecDuplicate(sp->vecs[0],&l);
399:   } else {
400:     MatGetVecs(mat,&l,NULL);
401:   }

403:   PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);
404:   if (sp->has_cnst) {
405:     VecDuplicate(l,&r);
406:     VecGetSize(l,&N);
407:     sum  = 1.0/N;
408:     VecSet(l,sum);
409:     MatMult(mat,l,r);
410:     VecNorm(r,NORM_2,&nrm);
411:     if (nrm >= 1.e-7) consistent = PETSC_FALSE;
412:     if (flg1) {
413:       if (consistent) {
414:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");
415:       } else {
416:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");
417:       }
418:       PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);
419:     }
420:     if (!consistent && flg1) {VecView(r,viewer);}
421:     if (!consistent && flg2) {VecView(r,viewer);}
422:     VecDestroy(&r);
423:   }

425:   for (j=0; j<n; j++) {
426:     (*mat->ops->mult)(mat,sp->vecs[j],l);
427:     VecNorm(l,NORM_2,&nrm);
428:     if (nrm >= 1.e-7) consistent = PETSC_FALSE;
429:     if (flg1) {
430:       if (consistent) {
431:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);
432:       } else {
433:         PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);
434:         consistent = PETSC_FALSE;
435:       }
436:       PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);
437:     }
438:     if (!consistent && flg1) {VecView(l,viewer);}
439:     if (!consistent && flg2) {VecView(l,viewer);}
440:   }

442:   if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
443:   VecDestroy(&l);
444:   if (isNull) *isNull = consistent;
445:   return(0);
446: }