Actual source code: bvbiorthog.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: BV bi-orthogonalization routines
12: */
14: #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
16: /*
17: BVBiorthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt bi-orthogonalization
18: */
19: static PetscErrorCode BVBiorthogonalizeMGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
20: {
22: PetscInt i;
23: PetscScalar dot;
24: Vec vi,wi;
27: for (i=-V->nc;i<V->k;i++) {
28: BVGetColumn(W,i,&wi);
29: /* h_i = ( v, w_i ) */
30: VecDot(v,wi,&dot);
31: BVRestoreColumn(W,i,&wi);
32: /* v <- v - h_i v_i */
33: BV_SetValue(V,i,0,c,dot);
34: BVGetColumn(V,i,&vi);
35: VecAXPY(v,-dot,vi);
36: BVRestoreColumn(V,i,&vi);
37: }
38: BV_AddCoefficients(V,V->k,h,c);
39: return(0);
40: }
42: /*
43: BVBiorthogonalizeCGS1 - Compute one step of CGS bi-orthogonalization: v = (I-V*W')*v
44: */
45: static PetscErrorCode BVBiorthogonalizeCGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
46: {
50: /* h = W'*v */
51: BVDotVec(W,v,c);
53: /* v = v - V h */
54: BVMultVec(V,-1.0,1.0,v,c);
56: BV_AddCoefficients(V,V->k,h,c);
57: return(0);
58: }
60: #define BVBiorthogonalizeGS1(a,b,c,d,e) ((V->orthog_type==BV_ORTHOG_MGS)?BVBiorthogonalizeMGS1:BVBiorthogonalizeCGS1)(a,b,c,d,e)
62: /*
63: BVBiorthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
65: V, W - the two basis vectors objects
66: v - the vector to bi-orthogonalize
67: */
68: static PetscErrorCode BVBiorthogonalizeGS(BV V,BV W,Vec v)
69: {
71: PetscScalar *h,*c;
74: h = V->h;
75: c = V->c;
76: BV_CleanCoefficients(V,V->k,h);
77: BVBiorthogonalizeGS1(V,W,v,h,c);
78: if (V->orthog_ref!=BV_ORTHOG_REFINE_NEVER) {
79: BVBiorthogonalizeGS1(V,W,v,h,c);
80: }
81: return(0);
82: }
84: /*@
85: BVBiorthogonalizeColumn - Bi-orthogonalize a column of two BV objects.
87: Collective on BV
89: Input Parameters:
90: + V,W - two basis vectors contexts
91: - j - index of column to be bi-orthonormalized
93: Notes:
94: This function bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
95: and V[0..j-1], respectively, so that W[0..j]'*V[0..j] = diagonal.
97: Level: advanced
99: .seealso: BVOrthogonalizeColumn(), BVBiorthonormalizeColumn()
100: @*/
101: PetscErrorCode BVBiorthogonalizeColumn(BV V,BV W,PetscInt j)
102: {
104: PetscInt ksavev,lsavev,ksavew,lsavew;
105: Vec y,z;
112: BVCheckSizes(V,1);
114: BVCheckSizes(W,2);
116: if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
117: if (j>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but V only has %D columns",j,V->m);
118: if (j>=W->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but W only has %D columns",j,W->m);
119: if (V->n!=W->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, W %D",V->n,W->n);
120: if (V->matrix || W->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
121: if (V->nc || W->nc) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
122: if (V->ops->gramschmidt || W->ops->gramschmidt) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
124: /* bi-orthogonalize */
125: PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0);
126: ksavev = V->k;
127: lsavev = V->l;
128: ksavew = W->k;
129: lsavew = W->l;
130: V->k = j;
131: V->l = -V->nc; /* must also bi-orthogonalize against constraints and leading columns */
132: W->k = j;
133: W->l = -W->nc;
134: BV_AllocateCoeffs(V);
135: BV_AllocateCoeffs(W);
136: BVGetColumn(V,j,&y);
137: BVBiorthogonalizeGS(V,W,y);
138: BVRestoreColumn(V,j,&y);
139: BVGetColumn(W,j,&z);
140: BVBiorthogonalizeGS(W,V,z);
141: BVRestoreColumn(W,j,&z);
142: V->k = ksavev;
143: V->l = lsavev;
144: W->k = ksavew;
145: W->l = lsavew;
146: PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0);
147: PetscObjectStateIncrease((PetscObject)V);
148: PetscObjectStateIncrease((PetscObject)W);
149: return(0);
150: }
151: /*@
152: BVBiorthonormalizeColumn - Bi-orthonormalize a column of two BV objects.
154: Collective on BV
156: Input Parameters:
157: + V,W - two basis vectors contexts
158: - j - index of column to be bi-orthonormalized
160: Output Parameters:
161: . delta - (optional) value used for normalization
163: Notes:
164: This function first bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
165: and V[0..j-1], respectively. Then, it scales the vectors with 1/delta, so
166: that the resulting vectors satisfy W[j]'*V[j] = 1.
168: Level: advanced
170: .seealso: BVOrthonormalizeColumn(), BVBiorthogonalizeColumn()
171: @*/
172: PetscErrorCode BVBiorthonormalizeColumn(BV V,BV W,PetscInt j,PetscReal *delta)
173: {
175: PetscScalar alpha;
176: PetscReal deltat;
177: PetscInt ksavev,lsavev,ksavew,lsavew;
178: Vec y,z;
185: BVCheckSizes(V,1);
187: BVCheckSizes(W,2);
189: if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
190: if (j>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but V only has %D columns",j,V->m);
191: if (j>=W->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but W only has %D columns",j,W->m);
192: if (V->n!=W->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, W %D",V->n,W->n);
193: if (V->matrix || W->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
194: if (V->nc || W->nc) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
195: if (V->ops->gramschmidt || W->ops->gramschmidt) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
197: /* bi-orthogonalize */
198: PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0);
199: ksavev = V->k;
200: lsavev = V->l;
201: ksavew = W->k;
202: lsavew = W->l;
203: V->k = j;
204: V->l = -V->nc; /* must also bi-orthogonalize against constraints and leading columns */
205: W->k = j;
206: W->l = -W->nc;
207: BV_AllocateCoeffs(V);
208: BV_AllocateCoeffs(W);
209: BVGetColumn(V,j,&y);
210: BVBiorthogonalizeGS(V,W,y);
211: BVRestoreColumn(V,j,&y);
212: BVGetColumn(W,j,&z);
213: BVBiorthogonalizeGS(W,V,z);
214: BVRestoreColumn(W,j,&z);
215: V->k = ksavev;
216: V->l = lsavev;
217: W->k = ksavew;
218: W->l = lsavew;
219: PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0);
221: /* scale */
222: PetscLogEventBegin(BV_Scale,V,0,0,0);
223: BVGetColumn(V,j,&y);
224: BVGetColumn(W,j,&z);
225: VecDot(y,z,&alpha);
226: BVRestoreColumn(V,j,&y);
227: BVRestoreColumn(W,j,&z);
228: deltat = PetscSqrtReal(PetscAbsScalar(alpha));
229: if (V->n) { (*V->ops->scale)(V,j,1.0/PetscConj(alpha/deltat)); }
230: if (W->n) { (*W->ops->scale)(W,j,1.0/deltat); }
231: PetscLogEventEnd(BV_Scale,V,0,0,0);
232: if (delta) *delta = deltat;
233: PetscObjectStateIncrease((PetscObject)V);
234: PetscObjectStateIncrease((PetscObject)W);
235: return(0);
236: }