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 int 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 int MatNullSpaceCreate(MPI_Comm comm,int has_cnst,int n,const Vec vecs[],MatNullSpace *SP) 39 { 40 MatNullSpace sp; 41 int ierr,i; 42 43 PetscFunctionBegin; 44 PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0); 45 PetscLogObjectCreate(sp); 46 PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace)); 47 48 sp->has_cnst = has_cnst; 49 sp->n = n; 50 sp->vec = PETSC_NULL; 51 if (n) { 52 ierr = PetscMalloc(n*sizeof(Vec),&sp->vecs);CHKERRQ(ierr); 53 for (i=0; i<n; i++) sp->vecs[i] = vecs[i]; 54 } else { 55 sp->vecs = 0; 56 } 57 58 for (i=0; i<n; i++) { 59 ierr = PetscObjectReference((PetscObject)sp->vecs[i]);CHKERRQ(ierr); 60 } 61 *SP = sp; 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "MatNullSpaceDestroy" 67 /*@ 68 MatNullSpaceDestroy - Destroys a data structure used to project vectors 69 out of null spaces. 70 71 Collective on MatNullSpace 72 73 Input Parameter: 74 . sp - the null space context to be destroyed 75 76 Level: advanced 77 78 .keywords: PC, null space, destroy 79 80 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove() 81 @*/ 82 int MatNullSpaceDestroy(MatNullSpace sp) 83 { 84 int ierr; 85 86 PetscFunctionBegin; 87 if (--sp->refct > 0) PetscFunctionReturn(0); 88 89 if (sp->vec) {ierr = VecDestroy(sp->vec);CHKERRQ(ierr);} 90 if (sp->vecs) { 91 ierr = VecDestroyVecs(sp->vecs,sp->n);CHKERRQ(ierr); 92 } 93 PetscLogObjectDestroy(sp); 94 PetscHeaderDestroy(sp); 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatNullSpaceRemove" 100 /*@ 101 MatNullSpaceRemove - Removes all the components of a null space from a vector. 102 103 Collective on MatNullSpace 104 105 Input Parameters: 106 + sp - the null space context 107 . vec - the vector from which the null space is to be removed 108 - out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise 109 the removal is done in-place (in vec) 110 111 112 113 Level: advanced 114 115 .keywords: PC, null space, remove 116 117 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 118 @*/ 119 int MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 120 { 121 PetscScalar sum; 122 int j,n = sp->n,N,ierr; 123 Vec l = vec; 124 125 PetscFunctionBegin; 126 if (out) { 127 if (!sp->vec) { 128 ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 129 } 130 *out = sp->vec; 131 ierr = VecCopy(vec,*out);CHKERRQ(ierr); 132 l = *out; 133 } 134 135 if (sp->has_cnst) { 136 ierr = VecSum(l,&sum);CHKERRQ(ierr); 137 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 138 sum = sum/(-1.0*N); 139 ierr = VecShift(&sum,l);CHKERRQ(ierr); 140 } 141 142 for (j=0; j<n; j++) { 143 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 144 sum = -sum; 145 ierr = VecAXPY(&sum,sp->vecs[j],l);CHKERRQ(ierr); 146 } 147 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "MatNullSpaceTest" 153 /*@ 154 MatNullSpaceTest - Tests if the claimed null space is really a 155 null space of a matrix 156 157 Collective on MatNullSpace 158 159 Input Parameters: 160 + sp - the null space context 161 - mat - the matrix 162 163 Level: advanced 164 165 .keywords: PC, null space, remove 166 167 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 168 @*/ 169 int MatNullSpaceTest(MatNullSpace sp,Mat mat) 170 { 171 PetscScalar sum; 172 PetscReal nrm; 173 int j,n = sp->n,N,ierr,m; 174 Vec l,r; 175 MPI_Comm comm = sp->comm; 176 PetscTruth flg1,flg2; 177 178 PetscFunctionBegin; 179 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 180 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 181 182 if (!sp->vec) { 183 if (n) { 184 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 185 } else { 186 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 187 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 188 } 189 } 190 l = sp->vec; 191 192 if (sp->has_cnst) { 193 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 194 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 195 sum = 1.0/N; 196 ierr = VecSet(&sum,l);CHKERRQ(ierr); 197 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 198 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 199 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 200 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 201 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 202 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 203 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 204 ierr = VecDestroy(r);CHKERRQ(ierr); 205 } 206 207 for (j=0; j<n; j++) { 208 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 209 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 210 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %d is likely null vector",j);CHKERRQ(ierr);} 211 else {ierr = PetscPrintf(comm,"Null vector %d unlikely null vector ",j);CHKERRQ(ierr);} 212 ierr = PetscPrintf(comm,"|| A * v[%d] || = %g\n",j,nrm);CHKERRQ(ierr); 213 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 214 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 215 } 216 217 PetscFunctionReturn(0); 218 } 219 220