Actual source code: matnull.c

      1: #define PETSCMAT_DLL

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

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

     10: PetscCookie PETSCMAT_DLLEXPORT 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: -  remove - the function that removes the null space

     24: .keywords: PC, null space, create

     26: .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceCreate()
     27: @*/
     28: PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*remove)(Vec,void*),void *ctx)
     29: {
     31:   sp->remove = remove;
     32:   sp->rmctx  = ctx;
     33:   return(0);
     34: }

     38: /*@C
     39:    MatNullSpaceCreate - Creates a data structure used to project vectors
     40:    out of null spaces.

     42:    Collective on MPI_Comm

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

     52:    Output Parameter:
     53: .  SP - the null space context

     55:    Level: advanced

     57:   Users manual sections:
     58: .   Section 4.15 Solving Singular Systems

     60: .keywords: PC, null space, create

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

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

     74:   sp->has_cnst = has_cnst;
     75:   sp->n        = n;
     76:   sp->vec      = PETSC_NULL;
     77:   if (n) {
     78:     PetscMalloc(n*sizeof(Vec),&sp->vecs);
     79:     for (i=0; i<n; i++) sp->vecs[i] = vecs[i];
     80:   } else {
     81:     sp->vecs = 0;
     82:   }

     84:   for (i=0; i<n; i++) {
     85:     PetscObjectReference((PetscObject)sp->vecs[i]);
     86:   }
     87:   *SP          = sp;
     88:   return(0);
     89: }

     93: /*@
     94:    MatNullSpaceDestroy - Destroys a data structure used to project vectors
     95:    out of null spaces.

     97:    Collective on MatNullSpace

     99:    Input Parameter:
    100: .  sp - the null space context to be destroyed

    102:    Level: advanced

    104: .keywords: PC, null space, destroy

    106: .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction()
    107: @*/
    108: PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceDestroy(MatNullSpace sp)
    109: {

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

    115:   if (sp->vec) {VecDestroy(sp->vec);}
    116:   if (sp->vecs) {
    117:     VecDestroyVecs(sp->vecs,sp->n);
    118:   }
    119:   PetscHeaderDestroy(sp);
    120:   return(0);
    121: }

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

    128:    Collective on MatNullSpace

    130:    Input Parameters:
    131: +  sp - the null space context
    132: .  vec - the vector from which the null space is to be removed
    133: -  out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise
    134:          the removal is done in-place (in vec)

    138:    Level: advanced

    140: .keywords: PC, null space, remove

    142: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
    143: @*/
    144: PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out)
    145: {
    146:   PetscScalar    sum;
    148:   PetscInt       j,n = sp->n,N;
    149:   Vec            l = vec;

    152:   if (out) {
    153:     if (!sp->vec) {
    154:       VecDuplicate(vec,&sp->vec);
    155:     }
    156:     *out = sp->vec;
    157:     VecCopy(vec,*out);
    158:     l    = *out;
    159:   }

    161:   if (sp->has_cnst) {
    162:     VecSum(l,&sum);
    163:     VecGetSize(l,&N);
    164:     sum  = sum/(-1.0*N);
    165:     VecShift(l,sum);
    166:   }

    168:   for (j=0; j<n; j++) {
    169:     VecDot(l,sp->vecs[j],&sum);
    170:     sum  = -sum;
    171:     VecAXPY(l,sum,sp->vecs[j]);
    172:   }

    174:   if (sp->remove){
    175:     (*sp->remove)(l,sp->rmctx);
    176:   }
    177:   return(0);
    178: }

    182: /*@
    183:    MatNullSpaceTest  - Tests if the claimed null space is really a
    184:      null space of a matrix

    186:    Collective on MatNullSpace

    188:    Input Parameters:
    189: +  sp - the null space context
    190: -  mat - the matrix

    192:    Level: advanced

    194: .keywords: PC, null space, remove

    196: .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
    197: @*/
    198: PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceTest(MatNullSpace sp,Mat mat)
    199: {
    200:   PetscScalar    sum;
    201:   PetscReal      nrm;
    202:   PetscInt       j,n = sp->n,N,m;
    204:   Vec            l,r;
    205:   MPI_Comm       comm = sp->comm;
    206:   PetscTruth     flg1,flg2;

    209:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);
    210:   PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);

    212:   if (!sp->vec) {
    213:     if (n) {
    214:       VecDuplicate(sp->vecs[0],&sp->vec);
    215:     } else {
    216:       MatGetLocalSize(mat,&m,PETSC_NULL);
    217:       VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);
    218:     }
    219:   }
    220:   l    = sp->vec;

    222:   if (sp->has_cnst) {
    223:     VecDuplicate(l,&r);
    224:     VecGetSize(l,&N);
    225:     sum  = 1.0/N;
    226:     VecSet(l,sum);
    227:     MatMult(mat,l,r);
    228:     VecNorm(r,NORM_2,&nrm);
    229:     if (nrm < 1.e-7) {PetscPrintf(comm,"Constants are likely null vector");}
    230:     else {PetscPrintf(comm,"Constants are unlikely null vector ");}
    231:     PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);
    232:     if (nrm > 1.e-7 && flg1) {VecView(r,PETSC_VIEWER_STDOUT_(comm));}
    233:     if (nrm > 1.e-7 && flg2) {VecView(r,PETSC_VIEWER_DRAW_(comm));}
    234:     VecDestroy(r);
    235:   }

    237:   for (j=0; j<n; j++) {
    238:     (*mat->ops->mult)(mat,sp->vecs[j],l);
    239:     VecNorm(l,NORM_2,&nrm);
    240:     if (nrm < 1.e-7) {PetscPrintf(comm,"Null vector %D is likely null vector",j);}
    241:     else {PetscPrintf(comm,"Null vector %D unlikely null vector ",j);}
    242:     PetscPrintf(comm,"|| A * v[%D] || = %g\n",j,nrm);
    243:     if (nrm > 1.e-7 && flg1) {VecView(l,PETSC_VIEWER_STDOUT_(comm));}
    244:     if (nrm > 1.e-7 && flg2) {VecView(l,PETSC_VIEWER_DRAW_(comm));}
    245:   }

    247:   if (sp->remove){
    248:     SETERRQ(PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()");
    249:   }
    250:
    251:   return(0);
    252: }