1 /* 2 Routines to project vectors out of null spaces. 3 */ 4 5 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 6 #include "petscsys.h" 7 8 PetscCookie MAT_NULLSPACE_COOKIE = 0; 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "MatNullSpaceCreate" 12 /*@C 13 MatNullSpaceCreate - Creates a data structure used to project vectors 14 out of null spaces. 15 16 Collective on MPI_Comm 17 18 Input Parameters: 19 + comm - the MPI communicator associated with the object 20 . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE 21 . n - number of vectors (excluding constant vector) in null space 22 - vecs - the vectors that span the null space (excluding the constant vector); 23 these vectors must be orthonormal. These vectors are NOT copied, so do not change them 24 after this call. You should free the array that you pass in. 25 26 Output Parameter: 27 . SP - the null space context 28 29 Level: advanced 30 31 Users manual sections: 32 . sec_singular 33 34 .keywords: PC, null space, create 35 36 .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace 37 @*/ 38 PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,int has_cnst,int n,const Vec vecs[],MatNullSpace *SP) 39 { 40 MatNullSpace sp; 41 PetscErrorCode ierr; 42 int 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 PetscErrorCode MatNullSpaceDestroy(MatNullSpace sp) 84 { 85 PetscErrorCode 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 PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 121 { 122 PetscScalar sum; 123 PetscErrorCode ierr; 124 int 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(&sum,l);CHKERRQ(ierr); 142 } 143 144 for (j=0; j<n; j++) { 145 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 146 sum = -sum; 147 ierr = VecAXPY(&sum,sp->vecs[j],l);CHKERRQ(ierr); 148 } 149 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "MatNullSpaceTest" 155 /*@ 156 MatNullSpaceTest - Tests if the claimed null space is really a 157 null space of a matrix 158 159 Collective on MatNullSpace 160 161 Input Parameters: 162 + sp - the null space context 163 - mat - the matrix 164 165 Level: advanced 166 167 .keywords: PC, null space, remove 168 169 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 170 @*/ 171 PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat) 172 { 173 PetscScalar sum; 174 PetscReal nrm; 175 int j,n = sp->n,N,m; 176 PetscErrorCode ierr; 177 Vec l,r; 178 MPI_Comm comm = sp->comm; 179 PetscTruth flg1,flg2; 180 181 PetscFunctionBegin; 182 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 183 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 184 185 if (!sp->vec) { 186 if (n) { 187 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 188 } else { 189 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 190 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 191 } 192 } 193 l = sp->vec; 194 195 if (sp->has_cnst) { 196 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 197 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 198 sum = 1.0/N; 199 ierr = VecSet(&sum,l);CHKERRQ(ierr); 200 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 201 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 202 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 203 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 204 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 205 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 206 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 207 ierr = VecDestroy(r);CHKERRQ(ierr); 208 } 209 210 for (j=0; j<n; j++) { 211 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 212 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 213 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %d is likely null vector",j);CHKERRQ(ierr);} 214 else {ierr = PetscPrintf(comm,"Null vector %d unlikely null vector ",j);CHKERRQ(ierr);} 215 ierr = PetscPrintf(comm,"|| A * v[%d] || = %g\n",j,nrm);CHKERRQ(ierr); 216 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 217 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 218 } 219 220 PetscFunctionReturn(0); 221 } 222 223