Actual source code: bvbiorthog.c
slepc-3.21.0 2024-03-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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>
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: {
21: PetscInt i;
22: PetscScalar dot;
23: Vec vi,wi;
25: PetscFunctionBegin;
26: for (i=-V->nc;i<V->k;i++) {
27: PetscCall(BVGetColumn(W,i,&wi));
28: /* h_i = (v, w_i) */
29: PetscCall(VecDot(v,wi,&dot));
30: PetscCall(BVRestoreColumn(W,i,&wi));
31: /* v <- v - h_i v_i */
32: PetscCall(BV_SetValue(V,i,0,c,dot));
33: PetscCall(BVGetColumn(V,i,&vi));
34: PetscCall(VecAXPY(v,-dot,vi));
35: PetscCall(BVRestoreColumn(V,i,&vi));
36: }
37: PetscCall(BV_AddCoefficients(V,V->k,h,c));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: /*
42: BVBiorthogonalizeCGS1 - Compute one step of CGS bi-orthogonalization: v = (I-V*W')*v
43: */
44: static PetscErrorCode BVBiorthogonalizeCGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
45: {
46: PetscFunctionBegin;
47: /* h = W'*v */
48: PetscCall(BVDotVec(W,v,c));
50: /* v = v - V h */
51: PetscCall(BVMultVec(V,-1.0,1.0,v,c));
53: PetscCall(BV_AddCoefficients(V,V->k,h,c));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: #define BVBiorthogonalizeGS1(a,b,c,d,e) ((V->orthog_type==BV_ORTHOG_MGS)?BVBiorthogonalizeMGS1:BVBiorthogonalizeCGS1)(a,b,c,d,e)
59: /*
60: BVBiorthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
62: V, W - the two basis vectors objects
63: v - the vector to bi-orthogonalize
64: */
65: static PetscErrorCode BVBiorthogonalizeGS(BV V,BV W,Vec v)
66: {
67: PetscScalar *h,*c;
69: PetscFunctionBegin;
70: h = V->h;
71: c = V->c;
72: PetscCall(BV_CleanCoefficients(V,V->k,h));
73: PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
74: if (V->orthog_ref!=BV_ORTHOG_REFINE_NEVER) PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*@
79: BVBiorthogonalizeColumn - Bi-orthogonalize a column of two BV objects.
81: Collective
83: Input Parameters:
84: + V - first basis vectors context
85: . W - second basis vectors context
86: - j - index of column to be bi-orthonormalized
88: Notes:
89: This function bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
90: and V[0..j-1], respectively, so that W[0..j]'*V[0..j] = diagonal.
92: Level: advanced
94: .seealso: BVOrthogonalizeColumn(), BVBiorthonormalizeColumn()
95: @*/
96: PetscErrorCode BVBiorthogonalizeColumn(BV V,BV W,PetscInt j)
97: {
98: PetscInt ksavev,lsavev,ksavew,lsavew;
99: Vec y,z;
101: PetscFunctionBegin;
106: BVCheckSizes(V,1);
108: BVCheckSizes(W,2);
109: PetscCheckSameTypeAndComm(V,1,W,2);
110: PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
111: PetscCheck(j<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but V only has %" PetscInt_FMT " columns",j,V->m);
112: PetscCheck(j<W->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but W only has %" PetscInt_FMT " columns",j,W->m);
113: PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
114: PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
115: PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
116: PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
118: /* bi-orthogonalize */
119: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
120: ksavev = V->k;
121: lsavev = V->l;
122: ksavew = W->k;
123: lsavew = W->l;
124: V->k = j;
125: V->l = -V->nc; /* must also bi-orthogonalize against constraints and leading columns */
126: W->k = j;
127: W->l = -W->nc;
128: PetscCall(BV_AllocateCoeffs(V));
129: PetscCall(BV_AllocateCoeffs(W));
130: PetscCall(BVGetColumn(V,j,&y));
131: PetscCall(BVBiorthogonalizeGS(V,W,y));
132: PetscCall(BVRestoreColumn(V,j,&y));
133: PetscCall(BVGetColumn(W,j,&z));
134: PetscCall(BVBiorthogonalizeGS(W,V,z));
135: PetscCall(BVRestoreColumn(W,j,&z));
136: V->k = ksavev;
137: V->l = lsavev;
138: W->k = ksavew;
139: W->l = lsavew;
140: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
141: PetscCall(PetscObjectStateIncrease((PetscObject)V));
142: PetscCall(PetscObjectStateIncrease((PetscObject)W));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*@
147: BVBiorthonormalizeColumn - Bi-orthonormalize a column of two BV objects.
149: Collective
151: Input Parameters:
152: + V - first basis vectors context
153: . W - second basis vectors context
154: - j - index of column to be bi-orthonormalized
156: Output Parameters:
157: . delta - (optional) value used for normalization
159: Notes:
160: This function first bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
161: and V[0..j-1], respectively. Then, it scales the vectors with 1/delta, so
162: that the resulting vectors satisfy W[j]'*V[j] = 1.
164: Level: advanced
166: .seealso: BVOrthonormalizeColumn(), BVBiorthogonalizeColumn()
167: @*/
168: PetscErrorCode BVBiorthonormalizeColumn(BV V,BV W,PetscInt j,PetscReal *delta)
169: {
170: PetscScalar alpha;
171: PetscReal deltat;
172: PetscInt ksavev,lsavev,ksavew,lsavew;
173: Vec y,z;
175: PetscFunctionBegin;
180: BVCheckSizes(V,1);
182: BVCheckSizes(W,2);
183: PetscCheckSameTypeAndComm(V,1,W,2);
184: PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
185: PetscCheck(j<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but V only has %" PetscInt_FMT " columns",j,V->m);
186: PetscCheck(j<W->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but W only has %" PetscInt_FMT " columns",j,W->m);
187: PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
188: PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
189: PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
190: PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
192: /* bi-orthogonalize */
193: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
194: ksavev = V->k;
195: lsavev = V->l;
196: ksavew = W->k;
197: lsavew = W->l;
198: V->k = j;
199: V->l = -V->nc; /* must also bi-orthogonalize against constraints and leading columns */
200: W->k = j;
201: W->l = -W->nc;
202: PetscCall(BV_AllocateCoeffs(V));
203: PetscCall(BV_AllocateCoeffs(W));
204: PetscCall(BVGetColumn(V,j,&y));
205: PetscCall(BVBiorthogonalizeGS(V,W,y));
206: PetscCall(BVRestoreColumn(V,j,&y));
207: PetscCall(BVGetColumn(W,j,&z));
208: PetscCall(BVBiorthogonalizeGS(W,V,z));
209: PetscCall(BVRestoreColumn(W,j,&z));
210: V->k = ksavev;
211: V->l = lsavev;
212: W->k = ksavew;
213: W->l = lsavew;
214: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
216: /* scale */
217: PetscCall(PetscLogEventBegin(BV_Scale,V,0,0,0));
218: PetscCall(BVGetColumn(V,j,&y));
219: PetscCall(BVGetColumn(W,j,&z));
220: PetscCall(VecDot(z,y,&alpha));
221: PetscCall(BVRestoreColumn(V,j,&y));
222: PetscCall(BVRestoreColumn(W,j,&z));
223: deltat = PetscSqrtReal(PetscAbsScalar(alpha));
224: PetscUseTypeMethod(V,scale,j,1.0/PetscConj(alpha/deltat));
225: PetscUseTypeMethod(W,scale,j,1.0/deltat);
226: PetscCall(PetscLogEventEnd(BV_Scale,V,0,0,0));
227: if (delta) *delta = deltat;
228: PetscCall(PetscObjectStateIncrease((PetscObject)V));
229: PetscCall(PetscObjectStateIncrease((PetscObject)W));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }