1 /*$Id: matnull.c,v 1.40 2001/09/07 20:09:09 bsmith Exp $*/ 2 /* 3 Routines to project vectors out of null spaces. 4 */ 5 6 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 7 #include "petscsys.h" 8 9 int MAT_NULLSPACE_COOKIE; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatNullSpaceCreate" 13 /*@C 14 MatNullSpaceCreate - Creates a data structure used to project vectors 15 out of null spaces. 16 17 Collective on MPI_Comm 18 19 Input Parameters: 20 + comm - the MPI communicator associated with the object 21 . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE 22 . n - number of vectors (excluding constant vector) in null space 23 - vecs - the vectors that span the null space (excluding the constant vector); 24 these vectors must be orthonormal 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(), PCNullSpaceAttach(), MatNullSpace 37 @*/ 38 int MatNullSpaceCreate(MPI_Comm comm,int has_cnst,int n,Vec *vecs,MatNullSpace *SP) 39 { 40 MatNullSpace sp; 41 42 PetscFunctionBegin; 43 PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_COOKIE,0,"MatNullSpace",comm,MatNullSpaceDestroy,0); 44 PetscLogObjectCreate(sp); 45 PetscLogObjectMemory(sp,sizeof(struct _p_MatNullSpace)); 46 47 sp->has_cnst = has_cnst; 48 sp->n = n; 49 sp->vecs = vecs; 50 sp->vec = PETSC_NULL; 51 52 *SP = sp; 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatNullSpaceDestroy" 58 /*@ 59 MatNullSpaceDestroy - Destroys a data structure used to project vectors 60 out of null spaces. 61 62 Collective on MatNullSpace 63 64 Input Parameter: 65 . sp - the null space context to be destroyed 66 67 Level: advanced 68 69 .keywords: PC, null space, destroy 70 71 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove() 72 @*/ 73 int MatNullSpaceDestroy(MatNullSpace sp) 74 { 75 int ierr; 76 77 PetscFunctionBegin; 78 if (--sp->refct > 0) PetscFunctionReturn(0); 79 80 if (sp->vec) {ierr = VecDestroy(sp->vec);CHKERRQ(ierr);} 81 82 PetscLogObjectDestroy(sp); 83 PetscHeaderDestroy(sp); 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatNullSpaceRemove" 89 /*@ 90 MatNullSpaceRemove - Removes all the components of a null space from a vector. 91 92 Collective on MatNullSpace 93 94 Input Parameters: 95 + sp - the null space context 96 - vec - the vector from which the null space is to be removed 97 98 Level: advanced 99 100 .keywords: PC, null space, remove 101 102 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 103 @*/ 104 int MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out) 105 { 106 PetscScalar sum; 107 int j,n = sp->n,N,ierr; 108 Vec l = vec; 109 110 PetscFunctionBegin; 111 if (out) { 112 if (!sp->vec) { 113 ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 114 } 115 *out = sp->vec; 116 ierr = VecCopy(vec,*out);CHKERRQ(ierr); 117 l = *out; 118 } 119 120 if (sp->has_cnst) { 121 ierr = VecSum(l,&sum);CHKERRQ(ierr); 122 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 123 sum = sum/(-1.0*N); 124 ierr = VecShift(&sum,l);CHKERRQ(ierr); 125 } 126 127 for (j=0; j<n; j++) { 128 ierr = VecDot(l,sp->vecs[j],&sum);CHKERRQ(ierr); 129 sum = -sum; 130 ierr = VecAXPY(&sum,sp->vecs[j],l);CHKERRQ(ierr); 131 } 132 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "MatNullSpaceTest" 138 /*@ 139 MatNullSpaceTest - Tests if the claimed null space is really a 140 null space of a matrix 141 142 Collective on MatNullSpace 143 144 Input Parameters: 145 + sp - the null space context 146 - mat - the matrix 147 148 Level: advanced 149 150 .keywords: PC, null space, remove 151 152 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy() 153 @*/ 154 int MatNullSpaceTest(MatNullSpace sp,Mat mat) 155 { 156 PetscScalar sum; 157 PetscReal nrm; 158 int j,n = sp->n,N,ierr,m; 159 Vec l,r; 160 MPI_Comm comm = sp->comm; 161 PetscTruth flg1,flg2; 162 163 PetscFunctionBegin; 164 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view",&flg1);CHKERRQ(ierr); 165 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_null_space_test_view_draw",&flg2);CHKERRQ(ierr); 166 167 if (!sp->vec) { 168 if (n) { 169 ierr = VecDuplicate(sp->vecs[0],&sp->vec);CHKERRQ(ierr); 170 } else { 171 ierr = MatGetLocalSize(mat,&m,PETSC_NULL);CHKERRQ(ierr); 172 ierr = VecCreateMPI(sp->comm,m,PETSC_DETERMINE,&sp->vec);CHKERRQ(ierr); 173 } 174 } 175 l = sp->vec; 176 177 if (sp->has_cnst) { 178 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 179 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 180 sum = 1.0/N; 181 ierr = VecSet(&sum,l);CHKERRQ(ierr); 182 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 183 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 184 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Constants are likely null vector");CHKERRQ(ierr);} 185 else {ierr = PetscPrintf(comm,"Constants are unlikely null vector ");CHKERRQ(ierr);} 186 ierr = PetscPrintf(comm,"|| A * 1 || = %g\n",nrm);CHKERRQ(ierr); 187 if (nrm > 1.e-7 && flg1) {ierr = VecView(r,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 188 if (nrm > 1.e-7 && flg2) {ierr = VecView(r,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 189 ierr = VecDestroy(r);CHKERRQ(ierr); 190 } 191 192 for (j=0; j<n; j++) { 193 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 194 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 195 if (nrm < 1.e-7) {ierr = PetscPrintf(comm,"Null vector %d is likely null vector",j);CHKERRQ(ierr);} 196 else {ierr = PetscPrintf(comm,"Null vector %d unlikely null vector ",j);CHKERRQ(ierr);} 197 ierr = PetscPrintf(comm,"|| A * v[%d] || = %g\n",j,nrm);CHKERRQ(ierr); 198 if (nrm > 1.e-7 && flg1) {ierr = VecView(l,PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);} 199 if (nrm > 1.e-7 && flg2) {ierr = VecView(l,PETSC_VIEWER_DRAW_(comm));CHKERRQ(ierr);} 200 } 201 202 PetscFunctionReturn(0); 203 } 204 205