1 #ifndef lint 2 static char vcid[] = "$Id: pcnull.c,v 1.7 1996/12/16 20:10:56 balay Exp balay $"; 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 #undef __FUNC__ 16 #define __FUNC__ "PCNullSpaceCreate" 17 /*@C 18 PCNullSpaceCreate - Creates a data-structure used to project vectors 19 out of null spaces. 20 21 Input Parameters: 22 . comm - the MPI communicator associated with the object. 23 . has_cnst - if the null spaces contains the constant vector, PETSC_TRUE or PETSC_FALSE 24 . n - number of vectors (excluding constant vector) in null space 25 . vecs - the vectors that span the null space (excluding the constant vector) 26 . these vectors must be orthonormal 27 28 Output Parameter: 29 . SP - the null space context 30 31 32 .keywords: PC, Null space 33 @*/ 34 int PCNullSpaceCreate(MPI_Comm comm, int has_cnst, int n, Vec *vecs,PCNullSpace *SP) 35 { 36 PCNullSpace sp; 37 38 PetscHeaderCreate(sp,_PCNullSpace,PCNULLSPACE_COOKIE,0,comm); 39 PLogObjectCreate(sp); 40 PLogObjectMemory(sp,sizeof(struct _PCNullSpace)); 41 42 sp->has_cnst = has_cnst; 43 sp->n = n; 44 sp->vecs = vecs; 45 46 *SP = sp; 47 return 0; 48 } 49 50 #undef __FUNC__ 51 #define __FUNC__ "PCNullSpaceDestroy" 52 /*@ 53 PCNullSpaceDestroy - Destroys a data-structure used to project vectors 54 out of null spaces. 55 56 Input Parameter: 57 . SP - the null space context to be destroyed 58 59 .keywords: PC, Null space 60 @*/ 61 int PCNullSpaceDestroy(PCNullSpace sp) 62 { 63 PLogObjectDestroy(sp); 64 PetscHeaderDestroy(sp); 65 return 0; 66 } 67 68 #undef __FUNC__ 69 #define __FUNC__ "PCNullSpaceRemove" 70 /*@ 71 PCNullSpaceRemove - Removes all the components of a null space from a vector. 72 73 Input Parameters: 74 . sp - the null space context 75 . vec - the vector you want the null space removed from 76 77 78 .keywords: PC, Null space 79 @*/ 80 int PCNullSpaceRemove(PCNullSpace sp,Vec vec) 81 { 82 Scalar sum; 83 int j, n = sp->n, N,ierr; 84 85 if (sp->has_cnst) { 86 ierr = VecSum(vec,&sum); CHKERRQ(ierr); 87 ierr = VecGetSize(vec,&N); CHKERRQ(ierr); 88 sum = -sum/N; 89 ierr = VecShift(&sum,vec); CHKERRQ(ierr); 90 } 91 92 for ( j=0; j<n; j++ ) { 93 ierr = VecDot(vec,sp->vecs[j],&sum);CHKERRQ(ierr); 94 sum = -sum; 95 ierr = VecAYPX(&sum,sp->vecs[j],vec); CHKERRQ(ierr); 96 } 97 98 return 0; 99 } 100