1 #ifndef lint 2 static char vcid[] = "$Id: pcnull.c,v 1.5 1996/04/25 23:59:27 balay Exp bsmith $"; 3 #endif 4 /* 5 Routines to project vectors out of null spaces. 6 */ 7 8 #include "petsc.h" 9 #include "src/pc/pcimpl.h" /*I "pc.h" I*/ 10 #include <stdio.h> 11 #include "src/sys/nreg.h" 12 #include "sys.h" 13 14 15 /*@C 16 PCNullSpaceCreate - Creates a data-structure used to project vectors 17 out of null spaces. 18 19 Input Parameters: 20 . comm - the MPI communicator associated with the object. 21 . has_cnst - if the null spaces contains the constant vector, PETSC_TRUE or 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 30 .keywords: PC, Null space 31 @*/ 32 int PCNullSpaceCreate(MPI_Comm comm, int has_cnst, int n, Vec *vecs,PCNullSpace *SP) 33 { 34 PCNullSpace sp; 35 36 PetscHeaderCreate(sp,_PCNullSpace,PCNULLSPACE_COOKIE,0,comm); 37 PLogObjectCreate(sp); 38 PLogObjectMemory(sp,sizeof(struct _PCNullSpace)); 39 40 sp->has_cnst = has_cnst; 41 sp->n = n; 42 sp->vecs = vecs; 43 44 *SP = sp; 45 return 0; 46 } 47 48 /*@ 49 PCNullSpaceDestroy - Destroys a data-structure used to project vectors 50 out of null spaces. 51 52 Input Parameter: 53 . SP - the null space context to be destroyed 54 55 .keywords: PC, Null space 56 @*/ 57 int PCNullSpaceDestroy(PCNullSpace sp) 58 { 59 PLogObjectDestroy(sp); 60 PetscHeaderDestroy(sp); 61 return 0; 62 } 63 64 /*@ 65 PCNullSpaceRemove - Removes all the components of a null space from a vector. 66 67 Input Parameters: 68 . sp - the null space context 69 . vec - the vector you want the null space removed from 70 71 72 .keywords: PC, Null space 73 @*/ 74 int PCNullSpaceRemove(PCNullSpace sp,Vec vec) 75 { 76 Scalar sum; 77 int j, n = sp->n, N,ierr; 78 79 if (sp->has_cnst) { 80 ierr = VecSum(vec,&sum); CHKERRQ(ierr); 81 ierr = VecGetSize(vec,&N); CHKERRQ(ierr); 82 sum = -sum/N; 83 ierr = VecShift(&sum,vec); CHKERRQ(ierr); 84 } 85 86 for ( j=0; j<n; j++ ) { 87 ierr = VecDot(vec,sp->vecs[j],&sum);CHKERRQ(ierr); 88 sum = -sum; 89 ierr = VecAYPX(&sum,sp->vecs[j],vec); CHKERRQ(ierr); 90 } 91 92 return 0; 93 } 94