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