1 #define PETSCMAT_DLL 2 3 /* 4 Routines to project vectors out of null spaces. 5 */ 6 7 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 8 #include "petscsys.h" 9 10 PetscCookie PETSCMAT_DLLEXPORT MAT_NULLSPACE_COOKIE = 0; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatNullSpaceSetFunction" 14 /*@C 15 MatNullSpaceSetFunction - set a function that removes a null space from a vector 16 out of null spaces. 17 18 Collective on MatNullSpace 19 20 Input Parameters: 21 + sp - the null space object 22 - remove - the function that removes the null space 23 24 .keywords: PC, null space, create 25 26 .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceCreate() 27 @*/ 28 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*remove)(Vec,void*),void *ctx) 29 { 30 PetscFunctionBegin; 31 sp->remove = remove; 32 sp->rmctx = ctx; 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatNullSpaceCreate" 38 /*@C 39 MatNullSpaceCreate - Creates a data structure used to project vectors 40 out of null spaces. 41 42 Collective on MPI_Comm 43 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. 51 52 Output Parameter: 53 . SP - the null space context 54 55 Level: advanced 56 57 Users manual sections: 58 . sec_singular 59 60 .keywords: PC, null space, create 61 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; 67 PetscErrorCode ierr; 68 PetscInt i; 69 70 PetscFunctionBegin; 71 ierr = PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0);CHKERRQ(ierr); 72 ierr = PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace));CHKERRQ(ierr); 73 74 sp->has_cnst = has_cnst; 75 sp->n = n; 76 sp->vec = PETSC_NULL; 77 if (n) { 78 ierr = PetscMalloc(n*sizeof(Vec),&sp->vecs);CHKERRQ(ierr); 79 for (i=0; i<n; i++) sp->vecs[i] = vecs[i]; 80 } else { 81 sp->vecs = 0; 82 } 83 84 for (i=0; i<n; i++) { 85 ierr = PetscObjectReference((PetscObject)sp->vecs[i]);CHKERRQ(ierr); 86 } 87 *SP = sp; 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNCT__ 92 #define __FUNCT__ "MatNullSpaceDestroy" 93 /*@ 94 MatNullSpaceDestroy - Destroys a data structure used to project vectors 95 out of null spaces. 96 97 Collective on MatNullSpace 98 99 Input Parameter: 100 . sp - the null space context to be destroyed 101 102 Level: advanced 103 104 .keywords: PC, null space, destroy 105 106 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction() 107 @*/ 108 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceDestroy(MatNullSpace sp) 109 { 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 if (--sp->refct > 0) PetscFunctionReturn(0); 114 115 if (sp->vec) {ierr = VecDestroy(sp->vec);CHKERRQ(ierr);} 116 if (sp->vecs) { 117 ierr = VecDestroyVecs(sp->vecs,sp->n);CHKERRQ(ierr); 118 } 119 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 #undef __FUNCT__ 124 #define __FUNCT__ "MatNullSpaceRemove" 125 /*@ 126 MatNullSpaceRemove - Removes all the components of a null space from a vector. 127 128 Collective on MatNullSpace 129 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) 135 136 137 138 Level: advanced 139 140 .keywords: PC, null space, remove 141 142 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction() 143 @*/ 144 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 145 { 146 PetscScalar sum; 147 PetscErrorCode ierr; 148 PetscInt j,n = sp->n,N; 149 Vec l = vec; 150 151 PetscFunctionBegin; 152 if (out) { 153 if (!sp->vec) { 154 ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 155 } 156 *out = sp->vec; 157 ierr = VecCopy(vec,*out);CHKERRQ(ierr); 158 l = *out; 159 } 160 161 if (sp->has_cnst) { 162 ierr = VecSum(l,&sum);CHKERRQ(ierr); 163 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 164 sum = sum/(-1.0*N); 165 ierr = VecShift(l,sum);CHKERRQ(ierr); 166 } 167 168 for (j=0; j<n; j++) { 169 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 170 sum = -sum; 171 ierr = VecAXPY(l,sum,sp->vecs[j]);CHKERRQ(ierr); 172 } 173 174 if (sp->remove){ 175 ierr = (*sp->remove)(l,sp->rmctx); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "MatNullSpaceTest" 182 /*@ 183 MatNullSpaceTest - Tests if the claimed null space is really a 184 null space of a matrix 185 186 Collective on MatNullSpace 187 188 Input Parameters: 189 + sp - the null space context 190 - mat - the matrix 191 192 Level: advanced 193 194 .keywords: PC, null space, remove 195 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; 203 PetscErrorCode ierr; 204 Vec l,r; 205 MPI_Comm comm = sp->comm; 206 PetscTruth flg1,flg2; 207 208 PetscFunctionBegin; 209 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 210 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 211 212 if (!sp->vec) { 213 if (n) { 214 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 215 } else { 216 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 217 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 218 } 219 } 220 l = sp->vec; 221 222 if (sp->has_cnst) { 223 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 224 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 225 sum = 1.0/N; 226 ierr = VecSet(l,sum);CHKERRQ(ierr); 227 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 228 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 229 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 230 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 231 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 232 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 233 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 234 ierr = VecDestroy(r);CHKERRQ(ierr); 235 } 236 237 for (j=0; j<n; j++) { 238 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 239 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 240 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %D is likely null vector",j);CHKERRQ(ierr);} 241 else {ierr = PetscPrintf(comm,"Null vector %D unlikely null vector ",j);CHKERRQ(ierr);} 242 ierr = PetscPrintf(comm,"|| A * v[%D] || = %g\n",j,nrm);CHKERRQ(ierr); 243 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 244 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 245 } 246 247 if (sp->remove){ 248 SETERRQ(PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()"); 249 } 250 251 PetscFunctionReturn(0); 252 } 253 254