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)) 29 { 30 PetscFunctionBegin; 31 sp->remove = remove; 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatNullSpaceCreate" 37 /*@C 38 MatNullSpaceCreate - Creates a data structure used to project vectors 39 out of null spaces. 40 41 Collective on MPI_Comm 42 43 Input Parameters: 44 + comm - the MPI communicator associated with the object 45 . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE 46 . n - number of vectors (excluding constant vector) in null space 47 - vecs - the vectors that span the null space (excluding the constant vector); 48 these vectors must be orthonormal. These vectors are NOT copied, so do not change them 49 after this call. You should free the array that you pass in. 50 51 Output Parameter: 52 . SP - the null space context 53 54 Level: advanced 55 56 Users manual sections: 57 . sec_singular 58 59 .keywords: PC, null space, create 60 61 .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction() 62 @*/ 63 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceCreate(MPI_Comm comm,PetscTruth has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP) 64 { 65 MatNullSpace sp; 66 PetscErrorCode ierr; 67 PetscInt i; 68 69 PetscFunctionBegin; 70 ierr = PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0);CHKERRQ(ierr); 71 ierr = PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace));CHKERRQ(ierr); 72 73 sp->has_cnst = has_cnst; 74 sp->n = n; 75 sp->vec = PETSC_NULL; 76 if (n) { 77 ierr = PetscMalloc(n*sizeof(Vec),&sp->vecs);CHKERRQ(ierr); 78 for (i=0; i<n; i++) sp->vecs[i] = vecs[i]; 79 } else { 80 sp->vecs = 0; 81 } 82 83 for (i=0; i<n; i++) { 84 ierr = PetscObjectReference((PetscObject)sp->vecs[i]);CHKERRQ(ierr); 85 } 86 *SP = sp; 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "MatNullSpaceDestroy" 92 /*@ 93 MatNullSpaceDestroy - Destroys a data structure used to project vectors 94 out of null spaces. 95 96 Collective on MatNullSpace 97 98 Input Parameter: 99 . sp - the null space context to be destroyed 100 101 Level: advanced 102 103 .keywords: PC, null space, destroy 104 105 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction() 106 @*/ 107 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceDestroy(MatNullSpace sp) 108 { 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 if (--sp->refct > 0) PetscFunctionReturn(0); 113 114 if (sp->vec) {ierr = VecDestroy(sp->vec);CHKERRQ(ierr);} 115 if (sp->vecs) { 116 ierr = VecDestroyVecs(sp->vecs,sp->n);CHKERRQ(ierr); 117 } 118 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatNullSpaceRemove" 124 /*@ 125 MatNullSpaceRemove - Removes all the components of a null space from a vector. 126 127 Collective on MatNullSpace 128 129 Input Parameters: 130 + sp - the null space context 131 . vec - the vector from which the null space is to be removed 132 - out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise 133 the removal is done in-place (in vec) 134 135 136 137 Level: advanced 138 139 .keywords: PC, null space, remove 140 141 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction() 142 @*/ 143 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 144 { 145 PetscScalar sum; 146 PetscErrorCode ierr; 147 PetscInt j,n = sp->n,N; 148 Vec l = vec; 149 150 PetscFunctionBegin; 151 if (out) { 152 if (!sp->vec) { 153 ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 154 } 155 *out = sp->vec; 156 ierr = VecCopy(vec,*out);CHKERRQ(ierr); 157 l = *out; 158 } 159 160 if (sp->has_cnst) { 161 ierr = VecSum(l,&sum);CHKERRQ(ierr); 162 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 163 sum = sum/(-1.0*N); 164 ierr = VecShift(l,sum);CHKERRQ(ierr); 165 } 166 167 for (j=0; j<n; j++) { 168 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 169 sum = -sum; 170 ierr = VecAXPY(l,sum,sp->vecs[j]);CHKERRQ(ierr); 171 } 172 173 if (sp->remove){ 174 ierr = (*sp->remove)(l); 175 } 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "MatNullSpaceTest" 181 /*@ 182 MatNullSpaceTest - Tests if the claimed null space is really a 183 null space of a matrix 184 185 Collective on MatNullSpace 186 187 Input Parameters: 188 + sp - the null space context 189 - mat - the matrix 190 191 Level: advanced 192 193 .keywords: PC, null space, remove 194 195 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction() 196 @*/ 197 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceTest(MatNullSpace sp,Mat mat) 198 { 199 PetscScalar sum; 200 PetscReal nrm; 201 PetscInt j,n = sp->n,N,m; 202 PetscErrorCode ierr; 203 Vec l,r; 204 MPI_Comm comm = sp->comm; 205 PetscTruth flg1,flg2; 206 207 PetscFunctionBegin; 208 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 209 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 210 211 if (!sp->vec) { 212 if (n) { 213 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 214 } else { 215 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 216 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 217 } 218 } 219 l = sp->vec; 220 221 if (sp->has_cnst) { 222 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 223 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 224 sum = 1.0/N; 225 ierr = VecSet(l,sum);CHKERRQ(ierr); 226 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 227 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 228 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 229 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 230 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 231 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 232 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 233 ierr = VecDestroy(r);CHKERRQ(ierr); 234 } 235 236 for (j=0; j<n; j++) { 237 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 238 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 239 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %D is likely null vector",j);CHKERRQ(ierr);} 240 else {ierr = PetscPrintf(comm,"Null vector %D unlikely null vector ",j);CHKERRQ(ierr);} 241 ierr = PetscPrintf(comm,"|| A * v[%D] || = %g\n",j,nrm);CHKERRQ(ierr); 242 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 243 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 244 } 245 246 if (sp->remove){ 247 SETERRQ(PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()"); 248 } 249 250 PetscFunctionReturn(0); 251 } 252 253