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