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 orthogonalization routines
12 : */
13 :
14 : #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15 :
16 : /*
17 : BV_NormVecOrColumn - Compute the 2-norm of the working vector, irrespective of
18 : whether it is in a column or not
19 : */
20 11664 : static inline PetscErrorCode BV_NormVecOrColumn(BV bv,PetscInt j,Vec v,PetscReal *nrm)
21 : {
22 11664 : PetscFunctionBegin;
23 11664 : if (v) PetscCall(BVNormVec(bv,v,NORM_2,nrm));
24 10512 : else PetscCall(BVNormColumn(bv,j,NORM_2,nrm));
25 11664 : PetscFunctionReturn(PETSC_SUCCESS);
26 : }
27 :
28 : /*
29 : BVDotColumnInc - Same as BVDotColumn() but also including column j, which
30 : is multiplied by itself
31 : */
32 247348 : static inline PetscErrorCode BVDotColumnInc(BV X,PetscInt j,PetscScalar *q)
33 : {
34 247348 : PetscInt ksave;
35 247348 : Vec y;
36 :
37 247348 : PetscFunctionBegin;
38 247348 : PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
39 247348 : ksave = X->k;
40 247348 : X->k = j+1;
41 247348 : PetscCall(BVGetColumn(X,j,&y));
42 247348 : PetscUseTypeMethod(X,dotvec,y,q);
43 247348 : PetscCall(BVRestoreColumn(X,j,&y));
44 247348 : X->k = ksave;
45 247348 : PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
46 247348 : PetscFunctionReturn(PETSC_SUCCESS);
47 : }
48 :
49 : /*
50 : BVOrthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt
51 : */
52 19467 : static PetscErrorCode BVOrthogonalizeMGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onrm,PetscReal *nrm)
53 : {
54 19467 : PetscInt i;
55 19467 : PetscScalar dot;
56 19467 : PetscBool indef=bv->indef;
57 19467 : Vec vi,z,w=v;
58 19467 : const PetscScalar *omega;
59 :
60 19467 : PetscFunctionBegin;
61 19467 : if (!v) PetscCall(BVGetColumn(bv,j,&w));
62 19467 : if (onrm) PetscCall(BVNormVec(bv,w,NORM_2,onrm));
63 19467 : z = w;
64 19467 : if (indef) PetscCall(VecGetArrayRead(bv->omega,&omega));
65 207767 : for (i=-bv->nc;i<j;i++) {
66 188300 : if (which && i>=0 && !which[i]) continue;
67 82729 : PetscCall(BVGetColumn(bv,i,&vi));
68 : /* h_i = (v, v_i) */
69 82729 : if (bv->matrix) {
70 160 : PetscCall(BV_IPMatMult(bv,w));
71 160 : z = bv->Bx;
72 : }
73 82729 : PetscCall(VecDot(z,vi,&dot));
74 : /* v <- v - h_i v_i */
75 82729 : PetscCall(BV_SetValue(bv,i,0,c,dot));
76 82729 : if (indef) dot /= PetscRealPart(omega[bv->nc+i]);
77 82729 : PetscCall(VecAXPY(w,-dot,vi));
78 188300 : PetscCall(BVRestoreColumn(bv,i,&vi));
79 : }
80 19467 : if (nrm) PetscCall(BVNormVec(bv,w,NORM_2,nrm));
81 19467 : if (!v) PetscCall(BVRestoreColumn(bv,j,&w));
82 19467 : PetscCall(BV_AddCoefficients(bv,j,h,c));
83 19467 : if (indef) PetscCall(VecRestoreArrayRead(bv->omega,&omega));
84 19467 : PetscFunctionReturn(PETSC_SUCCESS);
85 : }
86 :
87 : /*
88 : BVOrthogonalizeCGS1 - Compute |v'| (estimated), |v| and one step of CGS with
89 : only one global synchronization
90 : */
91 302293 : static PetscErrorCode BVOrthogonalizeCGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
92 : {
93 302293 : PetscReal sum,beta;
94 :
95 302293 : PetscFunctionBegin;
96 : /* h = W^* v ; alpha = (v, v) */
97 302293 : bv->k = j;
98 302293 : if (onorm || norm) {
99 299929 : if (!v) {
100 247348 : PetscCall(BVDotColumnInc(bv,j,c));
101 247348 : PetscCall(BV_SquareRoot(bv,j,c,&beta));
102 : } else {
103 52581 : PetscCall(BVDotVec(bv,v,c));
104 52581 : PetscCall(BVNormVec(bv,v,NORM_2,&beta));
105 : }
106 : } else {
107 2364 : if (!v) PetscCall(BVDotColumn(bv,j,c));
108 0 : else PetscCall(BVDotVec(bv,v,c));
109 : }
110 :
111 : /* q = v - V h */
112 302293 : if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,j,c,PETSC_TRUE));
113 302293 : if (!v) PetscCall(BVMultColumn(bv,-1.0,1.0,j,c));
114 52581 : else PetscCall(BVMultVec(bv,-1.0,1.0,v,c));
115 302293 : if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,j,c,PETSC_FALSE));
116 :
117 : /* compute |v| */
118 302293 : if (onorm) *onorm = beta;
119 :
120 302293 : if (norm) {
121 299929 : if (PetscUnlikely(bv->indef)) PetscCall(BV_NormVecOrColumn(bv,j,v,norm));
122 : else {
123 : /* estimate |v'| from |v| */
124 296364 : PetscCall(BV_SquareSum(bv,j,c,&sum));
125 296364 : *norm = beta*beta-sum;
126 296364 : if (PetscUnlikely(*norm <= 0.0)) PetscCall(BV_NormVecOrColumn(bv,j,v,norm));
127 288299 : else *norm = PetscSqrtReal(*norm);
128 : }
129 : }
130 302293 : PetscCall(BV_AddCoefficients(bv,j,h,c));
131 302293 : PetscFunctionReturn(PETSC_SUCCESS);
132 : }
133 :
134 : #define BVOrthogonalizeGS1(a,b,c,d,e,f,g,h) (bv->ops->gramschmidt?(*bv->ops->gramschmidt):(mgs?BVOrthogonalizeMGS1:BVOrthogonalizeCGS1))(a,b,c,d,e,f,g,h)
135 :
136 : /*
137 : BVOrthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
138 :
139 : j - the index of the column to orthogonalize (cannot use both j and v)
140 : v - the vector to orthogonalize (cannot use both j and v)
141 : which - logical array indicating selected columns (only used in MGS)
142 : norm - (optional) norm of the vector after being orthogonalized
143 : lindep - (optional) flag indicating possible linear dependence
144 : */
145 228951 : static PetscErrorCode BVOrthogonalizeGS(BV bv,PetscInt j,Vec v,PetscBool *which,PetscReal *norm,PetscBool *lindep)
146 : {
147 228951 : PetscScalar *h,*c,*omega;
148 228951 : PetscReal onrm,nrm;
149 228951 : PetscInt k,l;
150 228951 : PetscBool mgs,dolindep,signature;
151 :
152 228951 : PetscFunctionBegin;
153 228951 : if (v) {
154 47369 : k = bv->k;
155 47369 : h = bv->h;
156 47369 : c = bv->c;
157 : } else {
158 : k = j;
159 : h = NULL;
160 : c = NULL;
161 : }
162 :
163 228951 : mgs = (bv->orthog_type==BV_ORTHOG_MGS)? PETSC_TRUE: PETSC_FALSE;
164 :
165 : /* if indefinite inner product, skip the computation of lindep */
166 228951 : if (bv->indef && lindep) *lindep = PETSC_FALSE;
167 228951 : dolindep = (!bv->indef && lindep)? PETSC_TRUE: PETSC_FALSE;
168 :
169 : /* if indefinite and we are orthogonalizing a column, the norm must always be computed */
170 228951 : signature = (bv->indef && !v)? PETSC_TRUE: PETSC_FALSE;
171 :
172 228951 : PetscCall(BV_CleanCoefficients(bv,k,h));
173 :
174 228951 : switch (bv->orthog_ref) {
175 :
176 222087 : case BV_ORTHOG_REFINE_IFNEEDED:
177 423276 : PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,&onrm,&nrm));
178 : /* repeat if ||q|| < eta ||h|| */
179 : l = 1;
180 331768 : while (l<3 && nrm && PetscAbsReal(nrm) < bv->orthog_eta*PetscAbsReal(onrm)) {
181 115005 : l++;
182 115005 : if (mgs||bv->indef) onrm = nrm;
183 442914 : PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,(mgs||bv->indef)?NULL:&onrm,&nrm));
184 : }
185 : /* linear dependence check: criterion not satisfied in the last iteration */
186 361646 : if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
187 : break;
188 :
189 78 : case BV_ORTHOG_REFINE_NEVER:
190 124 : PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL));
191 : /* compute ||v|| */
192 78 : if (norm || dolindep || signature) PetscCall(BV_NormVecOrColumn(bv,k,v,&nrm));
193 : /* linear dependence check: just test for exactly zero norm */
194 78 : if (dolindep) *lindep = PetscNot(nrm);
195 : break;
196 :
197 6786 : case BV_ORTHOG_REFINE_ALWAYS:
198 9104 : PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL));
199 15576 : PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,dolindep?&onrm:NULL,(norm||dolindep||signature)?&nrm:NULL));
200 : /* linear dependence check: criterion not satisfied in the second iteration */
201 7104 : if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
202 : break;
203 : }
204 228951 : if (signature) {
205 7593 : PetscCall(VecGetArray(bv->omega,&omega));
206 7593 : omega[bv->nc+k] = (nrm<0.0)? -1.0: 1.0;
207 7593 : PetscCall(VecRestoreArray(bv->omega,&omega));
208 : }
209 228951 : if (norm) {
210 213911 : *norm = nrm;
211 213911 : if (!v) { /* store norm value next to the orthogonalization coefficients */
212 181526 : if (dolindep && *lindep) PetscCall(BV_SetValue(bv,k,k,h,0.0));
213 180402 : else PetscCall(BV_SetValue(bv,k,k,h,nrm));
214 : }
215 : }
216 228951 : PetscFunctionReturn(PETSC_SUCCESS);
217 : }
218 :
219 : /*@
220 : BVOrthogonalizeVec - Orthogonalize a given vector with respect to all
221 : active columns.
222 :
223 : Collective
224 :
225 : Input Parameters:
226 : + bv - the basis vectors context
227 : - v - the vector
228 :
229 : Output Parameters:
230 : + H - (optional) coefficients computed during orthogonalization
231 : . norm - (optional) norm of the vector after being orthogonalized
232 : - lindep - (optional) flag indicating that refinement did not improve the quality
233 : of orthogonalization
234 :
235 : Notes:
236 : This function is equivalent to BVOrthogonalizeColumn() but orthogonalizes
237 : a vector as an argument rather than taking one of the BV columns. The
238 : vector is orthogonalized against all active columns (k) and the constraints.
239 : If H is given, it must have enough space to store k-l coefficients, where l
240 : is the number of leading columns.
241 :
242 : In the case of an indefinite inner product, the lindep parameter is not
243 : computed (set to false).
244 :
245 : Level: advanced
246 :
247 : .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVSetActiveColumns(), BVGetNumConstraints()
248 : @*/
249 47369 : PetscErrorCode BVOrthogonalizeVec(BV bv,Vec v,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
250 : {
251 47369 : PetscInt ksave,lsave;
252 :
253 47369 : PetscFunctionBegin;
254 47369 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
255 47369 : PetscValidHeaderSpecific(v,VEC_CLASSID,2);
256 47369 : PetscValidType(bv,1);
257 47369 : BVCheckSizes(bv,1);
258 47369 : PetscValidType(v,2);
259 47369 : PetscCheckSameComm(bv,1,v,2);
260 :
261 47369 : PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
262 47369 : ksave = bv->k;
263 47369 : lsave = bv->l;
264 47369 : bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
265 47369 : PetscCall(BV_AllocateCoeffs(bv));
266 47369 : PetscCall(BV_AllocateSignature(bv));
267 47369 : PetscCall(BVOrthogonalizeGS(bv,0,v,NULL,norm,lindep));
268 47369 : bv->k = ksave;
269 47369 : bv->l = lsave;
270 47369 : if (H) PetscCall(BV_StoreCoefficients(bv,bv->k,bv->h,H));
271 47369 : PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
272 47369 : PetscFunctionReturn(PETSC_SUCCESS);
273 : }
274 :
275 : /*@
276 : BVOrthogonalizeColumn - Orthogonalize one of the column vectors with respect to
277 : the previous ones.
278 :
279 : Collective
280 :
281 : Input Parameters:
282 : + bv - the basis vectors context
283 : - j - index of column to be orthogonalized
284 :
285 : Output Parameters:
286 : + H - (optional) coefficients computed during orthogonalization
287 : . norm - (optional) norm of the vector after being orthogonalized
288 : - lindep - (optional) flag indicating that refinement did not improve the quality
289 : of orthogonalization
290 :
291 : Notes:
292 : This function applies an orthogonal projector to project vector V[j] onto
293 : the orthogonal complement of the span of the columns V[0..j-1],
294 : where V[.] are the vectors of BV. The columns V[0..j-1] are assumed to be
295 : mutually orthonormal.
296 :
297 : Leading columns V[0..l-1] also participate in the orthogonalization, as well
298 : as the constraints. If H is given, it must have enough space to store
299 : j-l+1 coefficients (the last coefficient will contain the value norm, unless
300 : the norm argument is NULL).
301 :
302 : If a non-standard inner product has been specified with BVSetMatrix(),
303 : then the vector is B-orthogonalized, using the non-standard inner product
304 : defined by matrix B. The output vector satisfies V[j]'*B*V[0..j-1] = 0.
305 :
306 : This routine does not normalize the resulting vector, see BVOrthonormalizeColumn().
307 :
308 : In the case of an indefinite inner product, the lindep parameter is not
309 : computed (set to false).
310 :
311 : Level: advanced
312 :
313 : .seealso: BVSetOrthogonalization(), BVSetMatrix(), BVSetActiveColumns(), BVOrthogonalize(), BVOrthogonalizeVec(), BVGetNumConstraints(), BVOrthonormalizeColumn()
314 : @*/
315 72412 : PetscErrorCode BVOrthogonalizeColumn(BV bv,PetscInt j,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
316 : {
317 72412 : PetscInt ksave,lsave;
318 :
319 72412 : PetscFunctionBegin;
320 72412 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
321 217236 : PetscValidLogicalCollectiveInt(bv,j,2);
322 72412 : PetscValidType(bv,1);
323 72412 : BVCheckSizes(bv,1);
324 72412 : PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
325 72412 : PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,bv->m);
326 :
327 72412 : PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
328 72412 : ksave = bv->k;
329 72412 : lsave = bv->l;
330 72412 : bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
331 72412 : if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
332 72412 : PetscCall(BV_AllocateSignature(bv));
333 72412 : PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,norm,lindep));
334 72412 : bv->k = ksave;
335 72412 : bv->l = lsave;
336 72412 : if (H) PetscCall(BV_StoreCoefficients(bv,j,NULL,H));
337 72412 : PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
338 72412 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
339 72412 : PetscFunctionReturn(PETSC_SUCCESS);
340 : }
341 :
342 : /*@
343 : BVOrthonormalizeColumn - Orthonormalize one of the column vectors with respect to
344 : the previous ones.
345 :
346 : Collective
347 :
348 : Input Parameters:
349 : + bv - the basis vectors context
350 : . j - index of column to be orthonormalized
351 : - replace - whether it is allowed to set the vector randomly
352 :
353 : Output Parameters:
354 : + norm - (optional) norm of the vector after orthogonalization and before normalization
355 : - lindep - (optional) flag indicating that linear dependence was determined during
356 : orthogonalization
357 :
358 : Notes:
359 : This is equivalent to a call to BVOrthogonalizeColumn() followed by a
360 : call to BVScaleColumn() with the reciprocal of the norm.
361 :
362 : This function first orthogonalizes vector V[j] with respect to V[0..j-1],
363 : where V[.] are the vectors of BV. A byproduct of this computation is norm,
364 : the norm of the vector after orthogonalization. Secondly, it scales the
365 : vector with 1/norm, so that the resulting vector has unit norm.
366 :
367 : If after orthogonalization the vector V[j] is exactly zero, it cannot be normalized
368 : because norm=0. In that case, it could be left as zero or replaced by a random
369 : vector that is then orthonormalized. The latter is achieved by setting the
370 : argument replace to TRUE. The vector will be replaced by a random vector also
371 : if lindep was set to TRUE, even if the norm is not exactly zero.
372 :
373 : If the vector has been replaced by a random vector, the output arguments norm and
374 : lindep will be set according to the orthogonalization of this new vector.
375 :
376 : Level: advanced
377 :
378 : .seealso: BVOrthogonalizeColumn(), BVScaleColumn()
379 : @*/
380 100370 : PetscErrorCode BVOrthonormalizeColumn(BV bv,PetscInt j,PetscBool replace,PetscReal *norm,PetscBool *lindep)
381 : {
382 100370 : PetscScalar alpha;
383 100370 : PetscReal nrm;
384 100370 : PetscInt ksave,lsave;
385 100370 : PetscBool lndep;
386 :
387 100370 : PetscFunctionBegin;
388 100370 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
389 301110 : PetscValidLogicalCollectiveInt(bv,j,2);
390 100370 : PetscValidType(bv,1);
391 100370 : BVCheckSizes(bv,1);
392 100370 : PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
393 100370 : PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,bv->m);
394 :
395 : /* orthogonalize */
396 100370 : PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
397 100370 : ksave = bv->k;
398 100370 : lsave = bv->l;
399 100370 : bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
400 100370 : if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
401 100370 : PetscCall(BV_AllocateSignature(bv));
402 100370 : PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
403 100370 : if (replace && (nrm==0.0 || lndep)) {
404 1 : PetscCall(PetscInfo(bv,"Vector was linearly dependent, generating a new random vector\n"));
405 1 : PetscCall(BVSetRandomColumn(bv,j));
406 1 : PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
407 1 : if (nrm==0.0 || lndep) { /* yet another attempt */
408 0 : PetscCall(BVSetRandomColumn(bv,j));
409 0 : PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
410 : }
411 : }
412 100370 : bv->k = ksave;
413 100370 : bv->l = lsave;
414 100370 : PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
415 :
416 : /* scale */
417 100370 : if (nrm!=1.0 && nrm!=0.0) {
418 99743 : alpha = 1.0/nrm;
419 99743 : PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
420 99743 : PetscUseTypeMethod(bv,scale,j,alpha);
421 99743 : PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
422 : }
423 100370 : if (norm) *norm = nrm;
424 100370 : if (lindep) *lindep = lndep;
425 100370 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
426 100370 : PetscFunctionReturn(PETSC_SUCCESS);
427 : }
428 :
429 : /*@
430 : BVOrthogonalizeSomeColumn - Orthogonalize one of the column vectors with
431 : respect to some of the previous ones.
432 :
433 : Collective
434 :
435 : Input Parameters:
436 : + bv - the basis vectors context
437 : . j - index of column to be orthogonalized
438 : - which - logical array indicating selected columns
439 :
440 : Output Parameters:
441 : + H - (optional) coefficients computed during orthogonalization
442 : . norm - (optional) norm of the vector after being orthogonalized
443 : - lindep - (optional) flag indicating that refinement did not improve the quality
444 : of orthogonalization
445 :
446 : Notes:
447 : This function is similar to BVOrthogonalizeColumn(), but V[j] is
448 : orthogonalized only against columns V[i] having which[i]=PETSC_TRUE.
449 : The length of array which must be j at least.
450 :
451 : The use of this operation is restricted to MGS orthogonalization type.
452 :
453 : In the case of an indefinite inner product, the lindep parameter is not
454 : computed (set to false).
455 :
456 : Level: advanced
457 :
458 : .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization()
459 : @*/
460 8799 : PetscErrorCode BVOrthogonalizeSomeColumn(BV bv,PetscInt j,PetscBool *which,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
461 : {
462 8799 : PetscInt ksave,lsave;
463 :
464 8799 : PetscFunctionBegin;
465 8799 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
466 26397 : PetscValidLogicalCollectiveInt(bv,j,2);
467 8799 : PetscAssertPointer(which,3);
468 8799 : PetscValidType(bv,1);
469 8799 : BVCheckSizes(bv,1);
470 8799 : PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
471 8799 : PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,bv->m);
472 8799 : PetscCheck(bv->orthog_type==BV_ORTHOG_MGS,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation only available for MGS orthogonalization");
473 :
474 8799 : PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
475 8799 : ksave = bv->k;
476 8799 : lsave = bv->l;
477 8799 : bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
478 8799 : if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
479 8799 : PetscCall(BV_AllocateSignature(bv));
480 8799 : PetscCall(BVOrthogonalizeGS(bv,j,NULL,which,norm,lindep));
481 8799 : bv->k = ksave;
482 8799 : bv->l = lsave;
483 8799 : if (H) PetscCall(BV_StoreCoefficients(bv,j,NULL,H));
484 8799 : PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
485 8799 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
486 8799 : PetscFunctionReturn(PETSC_SUCCESS);
487 : }
488 :
489 : /*
490 : Block Gram-Schmidt: V2 = V2 - V1*R12, where R12 = V1'*V2
491 : */
492 124 : static PetscErrorCode BVOrthogonalize_BlockGS(BV V,Mat R)
493 : {
494 124 : BV V1;
495 :
496 124 : PetscFunctionBegin;
497 124 : PetscCall(BVGetSplit(V,&V1,NULL));
498 124 : PetscCall(BVDot(V,V1,R));
499 124 : PetscCall(BVMult(V,-1.0,1.0,V1,R));
500 124 : PetscCall(BVRestoreSplit(V,&V1,NULL));
501 124 : PetscFunctionReturn(PETSC_SUCCESS);
502 : }
503 :
504 : /*
505 : Orthogonalize a set of vectors with Gram-Schmidt, column by column.
506 : */
507 5142 : static PetscErrorCode BVOrthogonalize_GS(BV V,Mat R)
508 : {
509 5142 : PetscScalar *r=NULL;
510 5142 : PetscReal norm;
511 5142 : PetscInt j,ldr,lsave;
512 5142 : Vec v,w;
513 :
514 5142 : PetscFunctionBegin;
515 5142 : if (R) {
516 362 : PetscCall(MatDenseGetLDA(R,&ldr));
517 362 : PetscCall(MatDenseGetArray(R,&r));
518 : }
519 5142 : if (V->matrix) {
520 1514 : PetscCall(BVGetCachedBV(V,&V->cached));
521 1514 : PetscCall(BVSetActiveColumns(V->cached,V->l,V->k));
522 : }
523 41546 : for (j=V->l;j<V->k;j++) {
524 36404 : if (V->matrix && V->orthog_type==BV_ORTHOG_MGS) { /* fill cached BV */
525 0 : PetscCall(BVGetColumn(V->cached,j,&v));
526 0 : PetscCall(BVGetColumn(V,j,&w));
527 0 : PetscCall(MatMult(V->matrix,w,v));
528 0 : PetscCall(BVRestoreColumn(V,j,&w));
529 0 : PetscCall(BVRestoreColumn(V->cached,j,&v));
530 : }
531 36404 : if (R) {
532 6214 : PetscCall(BVOrthogonalizeColumn(V,j,NULL,&norm,NULL));
533 6214 : lsave = V->l;
534 6214 : V->l = -V->nc;
535 6214 : PetscCall(BV_StoreCoefficients(V,j,NULL,r+j*ldr));
536 6214 : V->l = lsave;
537 6214 : r[j+j*ldr] = norm;
538 30190 : } else PetscCall(BVOrthogonalizeColumn(V,j,NULL,&norm,NULL));
539 36404 : PetscCheck(norm,PetscObjectComm((PetscObject)V),PETSC_ERR_CONV_FAILED,"Breakdown in BVOrthogonalize due to a linearly dependent column");
540 36404 : if (V->matrix && V->orthog_type==BV_ORTHOG_CGS) { /* fill cached BV */
541 10129 : PetscCall(BVGetColumn(V->cached,j,&v));
542 10129 : PetscCall(VecCopy(V->Bx,v));
543 10129 : PetscCall(BVRestoreColumn(V->cached,j,&v));
544 : }
545 36404 : PetscCall(BVScaleColumn(V,j,1.0/norm));
546 : }
547 5142 : if (R) PetscCall(MatDenseRestoreArray(R,&r));
548 5142 : PetscFunctionReturn(PETSC_SUCCESS);
549 : }
550 :
551 : /*
552 : BV_GetBufferMat - Create auxiliary seqdense matrix that wraps the bv->buffer.
553 : */
554 470 : static inline PetscErrorCode BV_GetBufferMat(BV bv)
555 : {
556 470 : PetscInt ld;
557 470 : PetscScalar *array;
558 :
559 470 : PetscFunctionBegin;
560 470 : if (!bv->Abuffer) {
561 161 : if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
562 161 : ld = bv->m+bv->nc;
563 161 : PetscCall(VecGetArray(bv->buffer,&array));
564 161 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,bv->m,array,&bv->Abuffer));
565 161 : PetscCall(VecRestoreArray(bv->buffer,&array));
566 : }
567 470 : PetscFunctionReturn(PETSC_SUCCESS);
568 : }
569 :
570 : /*
571 : BV_StoreCoeffsBlock_Default - Copy the contents of the BV buffer to a dense Mat
572 : provided by the caller. Only columns l:k-1 are copied, restricting to the upper
573 : triangular part if tri=PETSC_TRUE.
574 : */
575 153 : static inline PetscErrorCode BV_StoreCoeffsBlock_Default(BV bv,Mat R,PetscBool tri)
576 : {
577 153 : const PetscScalar *bb;
578 153 : PetscScalar *rr;
579 153 : PetscInt j,ldr,ldb;
580 :
581 153 : PetscFunctionBegin;
582 153 : PetscCall(MatDenseGetLDA(R,&ldr));
583 153 : PetscCall(MatDenseGetArray(R,&rr));
584 153 : ldb = bv->m+bv->nc;
585 153 : PetscCall(VecGetArrayRead(bv->buffer,&bb));
586 1453 : for (j=bv->l;j<bv->k;j++) PetscCall(PetscArraycpy(rr+j*ldr,bb+j*ldb,(tri?(j+1):bv->k)+bv->nc));
587 153 : PetscCall(VecRestoreArrayRead(bv->buffer,&bb));
588 153 : PetscCall(MatDenseRestoreArray(R,&rr));
589 153 : PetscFunctionReturn(PETSC_SUCCESS);
590 : }
591 :
592 : /*
593 : Orthogonalize a set of vectors with Cholesky: R=chol(V'*V), Q=V*inv(R)
594 : */
595 105 : static PetscErrorCode BVOrthogonalize_Chol(BV V,Mat Rin)
596 : {
597 105 : Mat R,S;
598 :
599 105 : PetscFunctionBegin;
600 105 : PetscCall(BV_GetBufferMat(V));
601 105 : R = V->Abuffer;
602 105 : if (Rin) S = Rin; /* use Rin as a workspace for S */
603 73 : else S = R;
604 105 : if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
605 105 : PetscCall(BVDot(V,V,R));
606 105 : PetscCall(BVMatCholInv_LAPACK_Private(V,R,S));
607 105 : PetscCall(BVMultInPlace(V,S,V->l,V->k));
608 105 : if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
609 105 : PetscFunctionReturn(PETSC_SUCCESS);
610 : }
611 :
612 : /*
613 : Orthogonalize a set of vectors with the Tall-Skinny QR method
614 : */
615 187 : static PetscErrorCode BVOrthogonalize_TSQR(BV V,Mat Rin)
616 : {
617 187 : PetscScalar *pv,*r=NULL;
618 187 : PetscInt ldr;
619 187 : Mat R;
620 :
621 187 : PetscFunctionBegin;
622 187 : PetscCall(BV_GetBufferMat(V));
623 187 : R = V->Abuffer;
624 187 : if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
625 187 : PetscCall(MatDenseGetLDA(R,&ldr));
626 187 : PetscCall(MatDenseGetArray(R,&r));
627 187 : PetscCall(BVGetArray(V,&pv));
628 187 : PetscCall(BVOrthogonalize_LAPACK_TSQR(V,V->n,V->k-V->l,pv+(V->nc+V->l)*V->ld,V->ld,r+V->l*ldr+V->l,ldr));
629 187 : PetscCall(BVRestoreArray(V,&pv));
630 187 : PetscCall(MatDenseRestoreArray(R,&r));
631 187 : if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
632 187 : PetscFunctionReturn(PETSC_SUCCESS);
633 : }
634 :
635 : /*
636 : Orthogonalize a set of vectors with TSQR, but computing R only, then doing Q=V*inv(R)
637 : */
638 73 : static PetscErrorCode BVOrthogonalize_TSQRCHOL(BV V,Mat Rin)
639 : {
640 73 : PetscScalar *pv,*r=NULL;
641 73 : PetscInt ldr;
642 73 : Mat R,S;
643 :
644 73 : PetscFunctionBegin;
645 73 : PetscCall(BV_GetBufferMat(V));
646 73 : R = V->Abuffer;
647 73 : if (Rin) S = Rin; /* use Rin as a workspace for S */
648 57 : else S = R;
649 73 : if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
650 73 : PetscCall(MatDenseGetLDA(R,&ldr));
651 73 : PetscCall(MatDenseGetArray(R,&r));
652 73 : PetscCall(BVGetArray(V,&pv));
653 73 : PetscCall(BVOrthogonalize_LAPACK_TSQR_OnlyR(V,V->n,V->k-V->l,pv+(V->nc+V->l)*V->ld,V->ld,r+V->l*ldr+V->l,ldr));
654 73 : PetscCall(BVRestoreArray(V,&pv));
655 73 : PetscCall(MatDenseRestoreArray(R,&r));
656 73 : PetscCall(BVMatTriInv_LAPACK_Private(V,R,S));
657 73 : PetscCall(BVMultInPlace(V,S,V->l,V->k));
658 73 : if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
659 73 : PetscFunctionReturn(PETSC_SUCCESS);
660 : }
661 :
662 : /*
663 : Orthogonalize a set of vectors with SVQB
664 : */
665 105 : static PetscErrorCode BVOrthogonalize_SVQB(BV V,Mat Rin)
666 : {
667 105 : Mat R,S;
668 :
669 105 : PetscFunctionBegin;
670 105 : PetscCall(BV_GetBufferMat(V));
671 105 : R = V->Abuffer;
672 105 : if (Rin) S = Rin; /* use Rin as a workspace for S */
673 73 : else S = R;
674 105 : if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
675 105 : PetscCall(BVDot(V,V,R));
676 105 : PetscCall(BVMatSVQB_LAPACK_Private(V,R,S));
677 105 : PetscCall(BVMultInPlace(V,S,V->l,V->k));
678 105 : if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_FALSE));
679 105 : PetscFunctionReturn(PETSC_SUCCESS);
680 : }
681 :
682 : /*@
683 : BVOrthogonalize - Orthogonalize all columns (starting from the leading ones),
684 : that is, compute the QR decomposition.
685 :
686 : Collective
687 :
688 : Input Parameters:
689 : + V - basis vectors to be orthogonalized (or B-orthogonalized), modified on output
690 : - R - a sequential dense matrix (or NULL), on output the triangular factor of
691 : the QR decomposition
692 :
693 : Notes:
694 : On input, matrix R must be a square sequential dense Mat, with at least as many
695 : rows and columns as the number of active columns of V. The output satisfies
696 : V0 = V*R (where V0 represent the input V) and V'*V = I (or V'*B*V = I if an
697 : inner product matrix B has been specified with BVSetMatrix()).
698 :
699 : If V has leading columns, then they are not modified (are assumed to be already
700 : orthonormal) and the leading columns of R are not referenced. Let the
701 : decomposition be
702 : .vb
703 : [ V01 V02 ] = [ V1 V2 ] [ R11 R12 ]
704 : [ 0 R22 ]
705 : .ve
706 : then V1 is left unchanged (equal to V01) as well as R11 (it should satisfy
707 : V01 = V1*R11).
708 :
709 : Can pass NULL if R is not required.
710 :
711 : The method to be used for block orthogonalization can be set with
712 : BVSetOrthogonalization(). If set to GS, the computation is done column by
713 : column with successive calls to BVOrthogonalizeColumn(). Note that in the
714 : SVQB method the R factor is not upper triangular.
715 :
716 : If V is rank-deficient or very ill-conditioned, that is, one or more columns are
717 : (almost) linearly dependent with respect to the rest, then the algorithm may
718 : break down or result in larger numerical error. Linearly dependent columns are
719 : essentially replaced by random directions, and the corresponding diagonal entry
720 : in R is set to (nearly) zero.
721 :
722 : Level: intermediate
723 :
724 : .seealso: BVOrthogonalizeColumn(), BVOrthogonalizeVec(), BVSetMatrix(), BVSetActiveColumns(), BVSetOrthogonalization(), BVOrthogBlockType
725 : @*/
726 5612 : PetscErrorCode BVOrthogonalize(BV V,Mat R)
727 : {
728 5612 : PetscInt m,n;
729 :
730 5612 : PetscFunctionBegin;
731 5612 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
732 5612 : PetscValidType(V,1);
733 5612 : BVCheckSizes(V,1);
734 5612 : if (R) {
735 515 : PetscValidHeaderSpecific(R,MAT_CLASSID,2);
736 515 : PetscValidType(R,2);
737 515 : PetscCheckTypeName(R,MATSEQDENSE);
738 515 : PetscCall(MatGetSize(R,&m,&n));
739 515 : PetscCheck(m==n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument is not square, it has %" PetscInt_FMT " rows and %" PetscInt_FMT " columns",m,n);
740 515 : PetscCheck(n>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat size %" PetscInt_FMT " is smaller than the number of BV active columns %" PetscInt_FMT,n,V->k);
741 : }
742 5612 : PetscCheck(!V->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Not implemented for BV with constraints, use BVOrthogonalizeColumn() instead");
743 :
744 5612 : PetscCall(PetscLogEventBegin(BV_Orthogonalize,V,R,0,0));
745 5612 : switch (V->orthog_block) {
746 5142 : case BV_ORTHOG_BLOCK_GS: /* proceed column by column with Gram-Schmidt */
747 5142 : PetscCall(BVOrthogonalize_GS(V,R));
748 : break;
749 105 : case BV_ORTHOG_BLOCK_CHOL:
750 105 : PetscCall(BVOrthogonalize_Chol(V,R));
751 : break;
752 187 : case BV_ORTHOG_BLOCK_TSQR:
753 187 : PetscCheck(!V->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
754 187 : PetscCall(BVOrthogonalize_TSQR(V,R));
755 : break;
756 73 : case BV_ORTHOG_BLOCK_TSQRCHOL:
757 73 : PetscCheck(!V->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
758 73 : PetscCall(BVOrthogonalize_TSQRCHOL(V,R));
759 : break;
760 105 : case BV_ORTHOG_BLOCK_SVQB:
761 105 : PetscCall(BVOrthogonalize_SVQB(V,R));
762 : break;
763 : }
764 5612 : PetscCall(PetscLogEventEnd(BV_Orthogonalize,V,R,0,0));
765 5612 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
766 5612 : PetscFunctionReturn(PETSC_SUCCESS);
767 : }
|
