Actual source code: matnull.c

  1: #define PETSCMAT_DLL

  3: /*
  4:     Routines to project vectors out of null spaces.
  5: */

 7:  #include include/private/matimpl.h
 8:  #include petscsys.h

 10: PetscCookie  MAT_NULLSPACE_COOKIE = 0;

 14: /*@C
 15:    MatNullSpaceSetFunction - set a function that removes a null space from a vector
 16:    out of null spaces.

 18:    Collective on MatNullSpace

 20:    Input Parameters:
 21: +  sp - the null space object
 22: .  rem - the function that removes the null space
 23: -  ctx - context for the remove function

 25:    Level: advanced

 27: .keywords: PC, null space, create

 29: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
 30: @*/
 31: PetscErrorCode  MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(Vec,void*),void *ctx)
 32: {
 34:   sp->remove = rem;
 35:   sp->rmctx  = ctx;
 36:   return(0);
 37: }

 41: /*@
 42:    MatNullSpaceCreate - Creates a data structure used to project vectors 
 43:    out of null spaces.

 45:    Collective on MPI_Comm

 47:    Input Parameters:
 48: +  comm - the MPI communicator associated with the object
 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 - the vectors that span the null space (excluding the constant vector);
 52:           these vectors must be orthonormal. These vectors are NOT copied, so do not change them
 53:           after this call. You should free the array that you pass in.

 55:    Output Parameter:
 56: .  SP - the null space context

 58:    Level: advanced

 60:   Users manual sections:
 61: .   Section 4.15 Solving Singular Systems

 63: .keywords: PC, null space, create

 65: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
 66: @*/
 67: PetscErrorCode  MatNullSpaceCreate(MPI_Comm comm,PetscTruth has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
 68: {
 69:   MatNullSpace   sp;
 71:   PetscInt       i;

 77: 
 78:   *SP = PETSC_NULL;
 79: #ifndef PETSC_USE_DYNAMIC_LIBRARIES 
 80:   MatInitializePackage(PETSC_NULL);
 81: #endif 

 83:   PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0);
 84:   PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace));

 86:   sp->has_cnst = has_cnst;
 87:   sp->n        = n;
 88:   sp->vec      = PETSC_NULL;
 89:   if (n) {
 90:     PetscMalloc(n*sizeof(Vec),&sp->vecs);
 91:     for (i=0; i<n; i++) sp->vecs[i] = vecs[i];
 92:   } else {
 93:     sp->vecs = 0;
 94:   }

 96:   for (i=0; i<n; i++) {
 97:     PetscObjectReference((PetscObject)sp->vecs[i]);
 98:   }
 99:   *SP          = sp;
100:   return(0);
101: }

105: /*@
106:    MatNullSpaceDestroy - Destroys a data structure used to project vectors 
107:    out of null spaces.

109:    Collective on MatNullSpace

111:    Input Parameter:
112: .  sp - the null space context to be destroyed

114:    Level: advanced

116: .keywords: PC, null space, destroy

118: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
119: @*/
120: PetscErrorCode  MatNullSpaceDestroy(MatNullSpace sp)
121: {

125:   if (--sp->refct > 0) return(0);

127:   if (sp->vec) {VecDestroy(sp->vec);}
128:   if (sp->vecs) {
129:     VecDestroyVecs(sp->vecs,sp->n);
130:   }
131:   PetscHeaderDestroy(sp);
132:   return(0);
133: }

137: /*@
138:    MatNullSpaceRemove - Removes all the components of a null space from a vector.

140:    Collective on MatNullSpace

142:    Input Parameters:
143: +  sp - the null space context
144: .  vec - the vector from which the null space is to be removed 
145: -  out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise
146:          the removal is done in-place (in vec)

148:    Note: The user is not responsible for the vector returned and should not destroy it.

150:    Level: advanced

152: .keywords: PC, null space, remove

154: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
155: @*/
156: PetscErrorCode  MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out)
157: {
158:   PetscScalar    sum;
160:   PetscInt       j,n = sp->n,N;
161:   Vec            l = vec;


167:   if (out) {
169:     if (!sp->vec) {
170:       VecDuplicate(vec,&sp->vec);
171:     }
172:     *out = sp->vec;
173:     VecCopy(vec,*out);
174:     l    = *out;
175:   }

177:   if (sp->has_cnst) {
178:     VecSum(l,&sum);
179:     VecGetSize(l,&N);
180:     sum  = sum/(-1.0*N);
181:     VecShift(l,sum);
182:   }

184:   for (j=0; j<n; j++) {
185:     VecDot(l,sp->vecs[j],&sum);
186:     VecAXPY(l,-sum,sp->vecs[j]);
187:   }

189:   if (sp->remove){
190:     (*sp->remove)(l,sp->rmctx);
191:   }
192:   return(0);
193: }

197: /*@
198:    MatNullSpaceTest  - Tests if the claimed null space is really a
199:      null space of a matrix

201:    Collective on MatNullSpace

203:    Input Parameters:
204: +  sp - the null space context
205: -  mat - the matrix

207:    Level: advanced

209: .keywords: PC, null space, remove

211: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
212: @*/
213: PetscErrorCode  MatNullSpaceTest(MatNullSpace sp,Mat mat)
214: {
215:   PetscScalar    sum;
216:   PetscReal      nrm;
217:   PetscInt       j,n = sp->n,N,m;
219:   Vec            l,r;
220:   MPI_Comm       comm = sp->comm;
221:   PetscTruth     flg1,flg2;
222:   PetscViewer    viewer;

225:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);
226:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);

228:   if (!sp->vec) {
229:     if (n) {
230:       VecDuplicate(sp->vecs[0],&sp->vec);
231:     } else {
232:       MatGetLocalSize(mat,&m,PETSC_NULL);
233:       VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);
234:     }
235:   }
236:   l    = sp->vec;

238:   PetscViewerASCIIGetStdout(comm,&viewer);
239:   if (sp->has_cnst) {
240:     VecDuplicate(l,&r);
241:     VecGetSize(l,&N);
242:     sum  = 1.0/N;
243:     VecSet(l,sum);
244:     MatMult(mat,l,r);
245:     VecNorm(r,NORM_2,&nrm);
246:     if (nrm < 1.e-7) {PetscPrintf(comm,"Constants are likely null vector");}
247:     else {PetscPrintf(comm,"Constants are unlikely null vector ");}
248:     PetscPrintf(comm,"|| A * 1 || = %G\n",nrm);
249:     if (nrm > 1.e-7 && flg1) {VecView(r,viewer);}
250:     if (nrm > 1.e-7 && flg2) {VecView(r,viewer);}
251:     VecDestroy(r);
252:   }

254:   for (j=0; j<n; j++) {
255:     (*mat->ops->mult)(mat,sp->vecs[j],l);
256:     VecNorm(l,NORM_2,&nrm);
257:     if (nrm < 1.e-7) {PetscPrintf(comm,"Null vector %D is likely null vector",j);}
258:     else {PetscPrintf(comm,"Null vector %D unlikely null vector ",j);}
259:     PetscPrintf(comm,"|| A * v[%D] || = %G\n",j,nrm);
260:     if (nrm > 1.e-7 && flg1) {VecView(l,viewer);}
261:     if (nrm > 1.e-7 && flg2) {VecView(l,viewer);}
262:   }

264:   if (sp->remove){
265:     SETERRQ(PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
266:   }
267: 
268:   return(0);
269: }