1: /*$Id: matnull.c,v 1.40 2001/09/07 20:09:09 bsmith Exp $*/
  2: /*
  3:     Routines to project vectors out of null spaces.
  4: */

 6:  #include src/mat/matimpl.h
 7:  #include petscsys.h

  9: int MAT_NULLSPACE_COOKIE;

 11: #undef __FUNCT__  
 13: /*@C
 14:    MatNullSpaceCreate - Creates a data structure used to project vectors 
 15:    out of null spaces.

 17:    Collective on MPI_Comm

 19:    Input Parameters:
 20: +  comm - the MPI communicator associated with the object
 21: .  has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
 22: .  n - number of vectors (excluding constant vector) in null space
 23: -  vecs - the vectors that span the null space (excluding the constant vector);
 24:           these vectors must be orthonormal. These vectors are NOT copied, so do not change them
 25:           after this call. You should free the array that you pass in.

 27:    Output Parameter:
 28: .  SP - the null space context

 30:    Level: advanced

 32:   Users manual sections:
 33: .   Section 4.15 Solving Singular Systems

 35: .keywords: PC, null space, create

 37: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), PCNullSpaceAttach(), MatNullSpace
 38: @*/
 39: int MatNullSpaceCreate(MPI_Comm comm,int has_cnst,int n,Vec *vecs,MatNullSpace *SP)
 40: {
 41:   MatNullSpace sp;
 42:   int          ierr,i;

 45:   PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0);
 46:   PetscLogObjectCreate(sp);
 47:   PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace));

 49:   sp->has_cnst = has_cnst;
 50:   sp->n        = n;
 51:   sp->vec      = PETSC_NULL;
 52:   if (n) {
 53:     PetscMalloc(n*sizeof(Vec),&sp->vecs);
 54:     for (i=0; i<n; i++) sp->vecs[i] = vecs[i];
 55:   } else {
 56:     sp->vecs = 0;
 57:   }

 59:   for (i=0; i<n; i++) {
 60:     PetscObjectReference((PetscObject)sp->vecs[i]);
 61:   }
 62:   *SP          = sp;
 63:   return(0);
 64: }

 66: #undef __FUNCT__  
 68: /*@
 69:    MatNullSpaceDestroy - Destroys a data structure used to project vectors 
 70:    out of null spaces.

 72:    Collective on MatNullSpace

 74:    Input Parameter:
 75: .  sp - the null space context to be destroyed

 77:    Level: advanced

 79: .keywords: PC, null space, destroy

 81: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove()
 82: @*/
 83: int MatNullSpaceDestroy(MatNullSpace sp)
 84: {

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

 90:   if (sp->vec) {VecDestroy(sp->vec);}
 91:   if (sp->vecs) {
 92:     VecDestroyVecs(sp->vecs,sp->n);
 93:   }
 94:   PetscLogObjectDestroy(sp);
 95:   PetscHeaderDestroy(sp);
 96:   return(0);
 97: }

 99: #undef __FUNCT__  
101: /*@
102:    MatNullSpaceRemove - Removes all the components of a null space from a vector.

104:    Collective on MatNullSpace

106:    Input Parameters:
107: +  sp - the null space context
108: -  vec - the vector from which the null space is to be removed 

110:    Level: advanced

112: .keywords: PC, null space, remove

114: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy()
115: @*/
116: int MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out)
117: {
118:   PetscScalar sum;
119:   int         j,n = sp->n,N,ierr;
120:   Vec         l = vec;

123:   if (out) {
124:     if (!sp->vec) {
125:       VecDuplicate(vec,&sp->vec);
126:     }
127:     *out = sp->vec;
128:     VecCopy(vec,*out);
129:     l    = *out;
130:   }

132:   if (sp->has_cnst) {
133:     VecSum(l,&sum);
134:     VecGetSize(l,&N);
135:     sum  = sum/(-1.0*N);
136:     VecShift(&sum,l);
137:   }

139:   for (j=0; j<n; j++) {
140:     VecDot(l,sp->vecs[j],&sum);
141:     sum  = -sum;
142:     VecAXPY(&sum,sp->vecs[j],l);
143:   }
144: 
145:   return(0);
146: }

148: #undef __FUNCT__  
150: /*@
151:    MatNullSpaceTest  - Tests if the claimed null space is really a
152:      null space of a matrix

154:    Collective on MatNullSpace

156:    Input Parameters:
157: +  sp - the null space context
158: -  mat - the matrix

160:    Level: advanced

162: .keywords: PC, null space, remove

164: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy()
165: @*/
166: int MatNullSpaceTest(MatNullSpace sp,Mat mat)
167: {
168:   PetscScalar  sum;
169:   PetscReal    nrm;
170:   int          j,n = sp->n,N,ierr,m;
171:   Vec          l,r;
172:   MPI_Comm     comm = sp->comm;
173:   PetscTruth   flg1,flg2;

176:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);
177:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);

179:   if (!sp->vec) {
180:     if (n) {
181:       VecDuplicate(sp->vecs[0],&sp->vec);
182:     } else {
183:       MatGetLocalSize(mat,&m,PETSC_NULL);
184:       VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);
185:     }
186:   }
187:   l    = sp->vec;

189:   if (sp->has_cnst) {
190:     VecDuplicate(l,&r);
191:     VecGetSize(l,&N);
192:     sum  = 1.0/N;
193:     VecSet(&sum,l);
194:     MatMult(mat,l,r);
195:     VecNorm(r,NORM_2,&nrm);
196:     if (nrm < 1.e-7) {PetscPrintf(comm,"Constants are likely null vector");}
197:     else {PetscPrintf(comm,"Constants are unlikely null vector ");}
198:     PetscPrintf(comm,"|| A * 1 || = %gn",nrm);
199:     if (nrm > 1.e-7 && flg1) {VecView(r,PETSC_VIEWER_STDOUT_(comm));}
200:     if (nrm > 1.e-7 && flg2) {VecView(r,PETSC_VIEWER_DRAW_(comm));}
201:     VecDestroy(r);
202:   }

204:   for (j=0; j<n; j++) {
205:     (*mat->ops->mult)(mat,sp->vecs[j],l);
206:     VecNorm(l,NORM_2,&nrm);
207:     if (nrm < 1.e-7) {PetscPrintf(comm,"Null vector %d is likely null vector",j);}
208:     else {PetscPrintf(comm,"Null vector %d unlikely null vector ",j);}
209:     PetscPrintf(comm,"|| A * v[%d] || = %gn",j,nrm);
210:     if (nrm > 1.e-7 && flg1) {VecView(l,PETSC_VIEWER_STDOUT_(comm));}
211:     if (nrm > 1.e-7 && flg2) {VecView(l,PETSC_VIEWER_DRAW_(comm));}
212:   }
213: 
214:   return(0);
215: }