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