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 */ 5 6 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 7 #include "petscsys.h" 8 9 int MAT_NULLSPACE_COOKIE = 0; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatNullSpaceCreate" 13 /*@C 14 MatNullSpaceCreate - Creates a data structure used to project vectors 15 out of null spaces. 16 17 Collective on MPI_Comm 18 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. 26 27 Output Parameter: 28 . SP - the null space context 29 30 Level: advanced 31 32 Users manual sections: 33 . sec_singular 34 35 .keywords: PC, null space, create 36 37 .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace 38 @*/ 39 int MatNullSpaceCreate(MPI_Comm comm,int has_cnst,int n,const Vec vecs[],MatNullSpace *SP) 40 { 41 MatNullSpace sp; 42 int ierr,i; 43 44 PetscFunctionBegin; 45 PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0); 46 PetscLogObjectCreate(sp); 47 PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace)); 48 49 sp->has_cnst = has_cnst; 50 sp->n = n; 51 sp->vec = PETSC_NULL; 52 if (n) { 53 ierr = PetscMalloc(n*sizeof(Vec),&sp->vecs);CHKERRQ(ierr); 54 for (i=0; i<n; i++) sp->vecs[i] = vecs[i]; 55 } else { 56 sp->vecs = 0; 57 } 58 59 for (i=0; i<n; i++) { 60 ierr = PetscObjectReference((PetscObject)sp->vecs[i]);CHKERRQ(ierr); 61 } 62 *SP = sp; 63 PetscFunctionReturn(0); 64 } 65 66 #undef __FUNCT__ 67 #define __FUNCT__ "MatNullSpaceDestroy" 68 /*@ 69 MatNullSpaceDestroy - Destroys a data structure used to project vectors 70 out of null spaces. 71 72 Collective on MatNullSpace 73 74 Input Parameter: 75 . sp - the null space context to be destroyed 76 77 Level: advanced 78 79 .keywords: PC, null space, destroy 80 81 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove() 82 @*/ 83 int MatNullSpaceDestroy(MatNullSpace sp) 84 { 85 int ierr; 86 87 PetscFunctionBegin; 88 if (--sp->refct > 0) PetscFunctionReturn(0); 89 90 if (sp->vec) {ierr = VecDestroy(sp->vec);CHKERRQ(ierr);} 91 if (sp->vecs) { 92 ierr = VecDestroyVecs(sp->vecs,sp->n);CHKERRQ(ierr); 93 } 94 PetscLogObjectDestroy(sp); 95 PetscHeaderDestroy(sp); 96 PetscFunctionReturn(0); 97 } 98 99 #undef __FUNCT__ 100 #define __FUNCT__ "MatNullSpaceRemove" 101 /*@ 102 MatNullSpaceRemove - Removes all the components of a null space from a vector. 103 104 Collective on MatNullSpace 105 106 Input Parameters: 107 + sp - the null space context 108 . vec - the vector from which the null space is to be removed 109 - out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise 110 the removal is done in-place (in vec) 111 112 113 114 Level: advanced 115 116 .keywords: PC, null space, remove 117 118 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 119 @*/ 120 int MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 121 { 122 PetscScalar sum; 123 int j,n = sp->n,N,ierr; 124 Vec l = vec; 125 126 PetscFunctionBegin; 127 if (out) { 128 if (!sp->vec) { 129 ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 130 } 131 *out = sp->vec; 132 ierr = VecCopy(vec,*out);CHKERRQ(ierr); 133 l = *out; 134 } 135 136 if (sp->has_cnst) { 137 ierr = VecSum(l,&sum);CHKERRQ(ierr); 138 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 139 sum = sum/(-1.0*N); 140 ierr = VecShift(&sum,l);CHKERRQ(ierr); 141 } 142 143 for (j=0; j<n; j++) { 144 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 145 sum = -sum; 146 ierr = VecAXPY(&sum,sp->vecs[j],l);CHKERRQ(ierr); 147 } 148 149 PetscFunctionReturn(0); 150 } 151 152 #undef __FUNCT__ 153 #define __FUNCT__ "MatNullSpaceTest" 154 /*@ 155 MatNullSpaceTest - Tests if the claimed null space is really a 156 null space of a matrix 157 158 Collective on MatNullSpace 159 160 Input Parameters: 161 + sp - the null space context 162 - mat - the matrix 163 164 Level: advanced 165 166 .keywords: PC, null space, remove 167 168 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 169 @*/ 170 int MatNullSpaceTest(MatNullSpace sp,Mat mat) 171 { 172 PetscScalar sum; 173 PetscReal nrm; 174 int j,n = sp->n,N,ierr,m; 175 Vec l,r; 176 MPI_Comm comm = sp->comm; 177 PetscTruth flg1,flg2; 178 179 PetscFunctionBegin; 180 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 181 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 182 183 if (!sp->vec) { 184 if (n) { 185 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 186 } else { 187 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 188 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 189 } 190 } 191 l = sp->vec; 192 193 if (sp->has_cnst) { 194 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 195 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 196 sum = 1.0/N; 197 ierr = VecSet(&sum,l);CHKERRQ(ierr); 198 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 199 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 200 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 201 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 202 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 203 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 204 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 205 ierr = VecDestroy(r);CHKERRQ(ierr); 206 } 207 208 for (j=0; j<n; j++) { 209 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 210 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 211 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %d is likely null vector",j);CHKERRQ(ierr);} 212 else {ierr = PetscPrintf(comm,"Null vector %d unlikely null vector ",j);CHKERRQ(ierr);} 213 ierr = PetscPrintf(comm,"|| A * v[%d] || = %g\n",j,nrm);CHKERRQ(ierr); 214 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 215 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 216 } 217 218 PetscFunctionReturn(0); 219 } 220 221