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 ierr = PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0);CHKERRQ(ierr); 48 ierr = PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace));CHKERRQ(ierr); 49 50 sp->has_cnst = has_cnst; 51 sp->n = n; 52 sp->vec = PETSC_NULL; 53 if (n) { 54 ierr = PetscMalloc(n*sizeof(Vec),&sp->vecs);CHKERRQ(ierr); 55 for (i=0; i<n; i++) sp->vecs[i] = vecs[i]; 56 } else { 57 sp->vecs = 0; 58 } 59 60 for (i=0; i<n; i++) { 61 ierr = PetscObjectReference((PetscObject)sp->vecs[i]);CHKERRQ(ierr); 62 } 63 *SP = sp; 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "MatNullSpaceDestroy" 69 /*@ 70 MatNullSpaceDestroy - Destroys a data structure used to project vectors 71 out of null spaces. 72 73 Collective on MatNullSpace 74 75 Input Parameter: 76 . sp - the null space context to be destroyed 77 78 Level: advanced 79 80 .keywords: PC, null space, destroy 81 82 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove() 83 @*/ 84 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceDestroy(MatNullSpace sp) 85 { 86 PetscErrorCode ierr; 87 88 PetscFunctionBegin; 89 if (--sp->refct > 0) PetscFunctionReturn(0); 90 91 if (sp->vec) {ierr = VecDestroy(sp->vec);CHKERRQ(ierr);} 92 if (sp->vecs) { 93 ierr = VecDestroyVecs(sp->vecs,sp->n);CHKERRQ(ierr); 94 } 95 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 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 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 121 { 122 PetscScalar sum; 123 PetscErrorCode ierr; 124 PetscInt j,n = sp->n,N; 125 Vec l = vec; 126 127 PetscFunctionBegin; 128 if (out) { 129 if (!sp->vec) { 130 ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 131 } 132 *out = sp->vec; 133 ierr = VecCopy(vec,*out);CHKERRQ(ierr); 134 l = *out; 135 } 136 137 if (sp->has_cnst) { 138 ierr = VecSum(l,&sum);CHKERRQ(ierr); 139 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 140 sum = sum/(-1.0*N); 141 ierr = VecShift(l,sum);CHKERRQ(ierr); 142 } 143 144 for (j=0; j<n; j++) { 145 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 146 ierr = VecAXPY(l,-sum,sp->vecs[j]);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 PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceTest(MatNullSpace sp,Mat mat) 171 { 172 PetscScalar sum; 173 PetscReal nrm; 174 PetscInt j,n = sp->n,N,m; 175 PetscErrorCode ierr; 176 Vec l,r; 177 MPI_Comm comm = sp->comm; 178 PetscTruth flg1,flg2; 179 180 PetscFunctionBegin; 181 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 182 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 183 184 if (!sp->vec) { 185 if (n) { 186 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 187 } else { 188 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 189 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 190 } 191 } 192 l = sp->vec; 193 194 if (sp->has_cnst) { 195 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 196 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 197 sum = 1.0/N; 198 ierr = VecSet(l,sum);CHKERRQ(ierr); 199 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 200 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 201 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 202 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 203 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 204 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 205 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 206 ierr = VecDestroy(r);CHKERRQ(ierr); 207 } 208 209 for (j=0; j<n; j++) { 210 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 211 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 212 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %D is likely null vector",j);CHKERRQ(ierr);} 213 else {ierr = PetscPrintf(comm,"Null vector %D unlikely null vector ",j);CHKERRQ(ierr);} 214 ierr = PetscPrintf(comm,"|| A * v[%D] || = %g\n",j,nrm);CHKERRQ(ierr); 215 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 216 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 217 } 218 219 PetscFunctionReturn(0); 220 } 221 222