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