Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15 :
16 : /*
17 : BVBiorthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt bi-orthogonalization
18 : */
19 112 : static PetscErrorCode BVBiorthogonalizeMGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
20 : {
21 112 : PetscInt i;
22 112 : PetscScalar dot;
23 112 : Vec vi,wi;
24 :
25 112 : PetscFunctionBegin;
26 448 : for (i=-V->nc;i<V->k;i++) {
27 336 : PetscCall(BVGetColumn(W,i,&wi));
28 : /* h_i = (v, w_i) */
29 336 : PetscCall(VecDot(v,wi,&dot));
30 336 : PetscCall(BVRestoreColumn(W,i,&wi));
31 : /* v <- v - h_i v_i */
32 336 : PetscCall(BV_SetValue(V,i,0,c,dot));
33 336 : PetscCall(BVGetColumn(V,i,&vi));
34 336 : PetscCall(VecAXPY(v,-dot,vi));
35 336 : PetscCall(BVRestoreColumn(V,i,&vi));
36 : }
37 112 : PetscCall(BV_AddCoefficients(V,V->k,h,c));
38 112 : PetscFunctionReturn(PETSC_SUCCESS);
39 : }
40 :
41 : /*
42 : BVBiorthogonalizeCGS1 - Compute one step of CGS bi-orthogonalization: v = (I-V*W')*v
43 : */
44 1756 : static PetscErrorCode BVBiorthogonalizeCGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
45 : {
46 1756 : PetscFunctionBegin;
47 : /* h = W'*v */
48 1756 : PetscCall(BVDotVec(W,v,c));
49 :
50 : /* v = v - V h */
51 1756 : PetscCall(BVMultVec(V,-1.0,1.0,v,c));
52 :
53 1756 : PetscCall(BV_AddCoefficients(V,V->k,h,c));
54 1756 : PetscFunctionReturn(PETSC_SUCCESS);
55 : }
56 :
57 : #define BVBiorthogonalizeGS1(a,b,c,d,e) ((V->orthog_type==BV_ORTHOG_MGS)?BVBiorthogonalizeMGS1:BVBiorthogonalizeCGS1)(a,b,c,d,e)
58 :
59 : /*
60 : BVBiorthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
61 :
62 : V, W - the two basis vectors objects
63 : v - the vector to bi-orthogonalize
64 : */
65 934 : static PetscErrorCode BVBiorthogonalizeGS(BV V,BV W,Vec v)
66 : {
67 934 : PetscScalar *h,*c;
68 :
69 934 : PetscFunctionBegin;
70 934 : h = V->h;
71 934 : c = V->c;
72 934 : PetscCall(BV_CleanCoefficients(V,V->k,h));
73 1812 : PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
74 1812 : if (V->orthog_ref!=BV_ORTHOG_REFINE_NEVER) PetscCall(BVBiorthogonalizeGS1(V,W,v,h,c));
75 934 : PetscFunctionReturn(PETSC_SUCCESS);
76 : }
77 :
78 : /*@
79 : BVBiorthogonalizeColumn - Bi-orthogonalize a column of two BV objects.
80 :
81 : Collective
82 :
83 : Input Parameters:
84 : + V - first basis vectors context
85 : . W - second basis vectors context
86 : - j - index of column to be bi-orthonormalized
87 :
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.
91 :
92 : Level: advanced
93 :
94 : .seealso: BVOrthogonalizeColumn(), BVBiorthonormalizeColumn()
95 : @*/
96 394 : PetscErrorCode BVBiorthogonalizeColumn(BV V,BV W,PetscInt j)
97 : {
98 394 : PetscInt ksavev,lsavev,ksavew,lsavew;
99 394 : Vec y,z;
100 :
101 394 : PetscFunctionBegin;
102 394 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
103 394 : PetscValidHeaderSpecific(W,BV_CLASSID,2);
104 1182 : PetscValidLogicalCollectiveInt(V,j,3);
105 394 : PetscValidType(V,1);
106 394 : BVCheckSizes(V,1);
107 394 : PetscValidType(W,2);
108 394 : BVCheckSizes(W,2);
109 394 : PetscCheckSameTypeAndComm(V,1,W,2);
110 394 : PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
111 394 : 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 394 : 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 394 : 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 394 : PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
115 394 : PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
116 394 : PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
117 :
118 : /* bi-orthogonalize */
119 394 : PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
120 394 : ksavev = V->k;
121 394 : lsavev = V->l;
122 394 : ksavew = W->k;
123 394 : lsavew = W->l;
124 394 : V->k = j;
125 394 : V->l = -V->nc; /* must also bi-orthogonalize against constraints and leading columns */
126 394 : W->k = j;
127 394 : W->l = -W->nc;
128 394 : PetscCall(BV_AllocateCoeffs(V));
129 394 : PetscCall(BV_AllocateCoeffs(W));
130 394 : PetscCall(BVGetColumn(V,j,&y));
131 394 : PetscCall(BVBiorthogonalizeGS(V,W,y));
132 394 : PetscCall(BVRestoreColumn(V,j,&y));
133 394 : PetscCall(BVGetColumn(W,j,&z));
134 394 : PetscCall(BVBiorthogonalizeGS(W,V,z));
135 394 : PetscCall(BVRestoreColumn(W,j,&z));
136 394 : V->k = ksavev;
137 394 : V->l = lsavev;
138 394 : W->k = ksavew;
139 394 : W->l = lsavew;
140 394 : PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
141 394 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
142 394 : PetscCall(PetscObjectStateIncrease((PetscObject)W));
143 394 : PetscFunctionReturn(PETSC_SUCCESS);
144 : }
145 :
146 : /*@
147 : BVBiorthonormalizeColumn - Bi-orthonormalize a column of two BV objects.
148 :
149 : Collective
150 :
151 : Input Parameters:
152 : + V - first basis vectors context
153 : . W - second basis vectors context
154 : - j - index of column to be bi-orthonormalized
155 :
156 : Output Parameters:
157 : . delta - (optional) value used for normalization
158 :
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.
163 :
164 : Level: advanced
165 :
166 : .seealso: BVOrthonormalizeColumn(), BVBiorthogonalizeColumn()
167 : @*/
168 73 : PetscErrorCode BVBiorthonormalizeColumn(BV V,BV W,PetscInt j,PetscReal *delta)
169 : {
170 73 : PetscScalar alpha;
171 73 : PetscReal deltat;
172 73 : PetscInt ksavev,lsavev,ksavew,lsavew;
173 73 : Vec y,z;
174 :
175 73 : PetscFunctionBegin;
176 73 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
177 73 : PetscValidHeaderSpecific(W,BV_CLASSID,2);
178 219 : PetscValidLogicalCollectiveInt(V,j,3);
179 73 : PetscValidType(V,1);
180 73 : BVCheckSizes(V,1);
181 73 : PetscValidType(W,2);
182 73 : BVCheckSizes(W,2);
183 73 : PetscCheckSameTypeAndComm(V,1,W,2);
184 73 : PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
185 73 : 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 73 : 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 73 : 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 73 : PetscCheck(!V->matrix && !W->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
189 73 : PetscCheck(!V->nc && !W->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
190 73 : PetscCheck(!V->ops->gramschmidt && !W->ops->gramschmidt,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");
191 :
192 : /* bi-orthogonalize */
193 73 : PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0));
194 73 : ksavev = V->k;
195 73 : lsavev = V->l;
196 73 : ksavew = W->k;
197 73 : lsavew = W->l;
198 73 : V->k = j;
199 73 : V->l = -V->nc; /* must also bi-orthogonalize against constraints and leading columns */
200 73 : W->k = j;
201 73 : W->l = -W->nc;
202 73 : PetscCall(BV_AllocateCoeffs(V));
203 73 : PetscCall(BV_AllocateCoeffs(W));
204 73 : PetscCall(BVGetColumn(V,j,&y));
205 73 : PetscCall(BVBiorthogonalizeGS(V,W,y));
206 73 : PetscCall(BVRestoreColumn(V,j,&y));
207 73 : PetscCall(BVGetColumn(W,j,&z));
208 73 : PetscCall(BVBiorthogonalizeGS(W,V,z));
209 73 : PetscCall(BVRestoreColumn(W,j,&z));
210 73 : V->k = ksavev;
211 73 : V->l = lsavev;
212 73 : W->k = ksavew;
213 73 : W->l = lsavew;
214 73 : PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0));
215 :
216 : /* scale */
217 73 : PetscCall(PetscLogEventBegin(BV_Scale,V,0,0,0));
218 73 : PetscCall(BVGetColumn(V,j,&y));
219 73 : PetscCall(BVGetColumn(W,j,&z));
220 73 : PetscCall(VecDot(z,y,&alpha));
221 73 : PetscCall(BVRestoreColumn(V,j,&y));
222 73 : PetscCall(BVRestoreColumn(W,j,&z));
223 73 : deltat = PetscSqrtReal(PetscAbsScalar(alpha));
224 73 : PetscUseTypeMethod(V,scale,j,1.0/PetscConj(alpha/deltat));
225 73 : PetscUseTypeMethod(W,scale,j,1.0/deltat);
226 73 : PetscCall(PetscLogEventEnd(BV_Scale,V,0,0,0));
227 73 : if (delta) *delta = deltat;
228 73 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
229 73 : PetscCall(PetscObjectStateIncrease((PetscObject)W));
230 73 : PetscFunctionReturn(PETSC_SUCCESS);
231 : }
|
