Actual source code: bvorthog.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 orthogonalization routines
12: */
14: #include <slepc/private/bvimpl.h>
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: static inline PetscErrorCode BV_NormVecOrColumn(BV bv,PetscInt j,Vec v,PetscReal *nrm)
21: {
22: PetscFunctionBegin;
23: if (v) PetscCall(BVNormVec(bv,v,NORM_2,nrm));
24: else PetscCall(BVNormColumn(bv,j,NORM_2,nrm));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*
29: BVDotColumnInc - Same as BVDotColumn() but also including column j, which
30: is multiplied by itself
31: */
32: static inline PetscErrorCode BVDotColumnInc(BV X,PetscInt j,PetscScalar *q)
33: {
34: PetscInt ksave;
35: Vec y;
37: PetscFunctionBegin;
38: PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
39: ksave = X->k;
40: X->k = j+1;
41: PetscCall(BVGetColumn(X,j,&y));
42: PetscUseTypeMethod(X,dotvec,y,q);
43: PetscCall(BVRestoreColumn(X,j,&y));
44: X->k = ksave;
45: PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: /*
50: BVOrthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt
51: */
52: static PetscErrorCode BVOrthogonalizeMGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onrm,PetscReal *nrm)
53: {
54: PetscInt i;
55: PetscScalar dot;
56: PetscBool indef=bv->indef;
57: Vec vi,z,w=v;
58: const PetscScalar *omega;
60: PetscFunctionBegin;
61: if (!v) PetscCall(BVGetColumn(bv,j,&w));
62: if (onrm) PetscCall(BVNormVec(bv,w,NORM_2,onrm));
63: z = w;
64: if (indef) PetscCall(VecGetArrayRead(bv->omega,&omega));
65: for (i=-bv->nc;i<j;i++) {
66: if (which && i>=0 && !which[i]) continue;
67: PetscCall(BVGetColumn(bv,i,&vi));
68: /* h_i = (v, v_i) */
69: if (bv->matrix) {
70: PetscCall(BV_IPMatMult(bv,w));
71: z = bv->Bx;
72: }
73: PetscCall(VecDot(z,vi,&dot));
74: /* v <- v - h_i v_i */
75: PetscCall(BV_SetValue(bv,i,0,c,dot));
76: if (indef) dot /= PetscRealPart(omega[bv->nc+i]);
77: PetscCall(VecAXPY(w,-dot,vi));
78: PetscCall(BVRestoreColumn(bv,i,&vi));
79: }
80: if (nrm) PetscCall(BVNormVec(bv,w,NORM_2,nrm));
81: if (!v) PetscCall(BVRestoreColumn(bv,j,&w));
82: PetscCall(BV_AddCoefficients(bv,j,h,c));
83: if (indef) PetscCall(VecRestoreArrayRead(bv->omega,&omega));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /*
88: BVOrthogonalizeCGS1 - Compute |v'| (estimated), |v| and one step of CGS with
89: only one global synchronization
90: */
91: static PetscErrorCode BVOrthogonalizeCGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
92: {
93: PetscReal sum,beta;
95: PetscFunctionBegin;
96: /* h = W^* v ; alpha = (v, v) */
97: bv->k = j;
98: if (onorm || norm) {
99: if (!v) {
100: PetscCall(BVDotColumnInc(bv,j,c));
101: PetscCall(BV_SquareRoot(bv,j,c,&beta));
102: } else {
103: PetscCall(BVDotVec(bv,v,c));
104: PetscCall(BVNormVec(bv,v,NORM_2,&beta));
105: }
106: } else {
107: if (!v) PetscCall(BVDotColumn(bv,j,c));
108: else PetscCall(BVDotVec(bv,v,c));
109: }
111: /* q = v - V h */
112: if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,j,c,PETSC_TRUE));
113: if (!v) PetscCall(BVMultColumn(bv,-1.0,1.0,j,c));
114: else PetscCall(BVMultVec(bv,-1.0,1.0,v,c));
115: if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,j,c,PETSC_FALSE));
117: /* compute |v| */
118: if (onorm) *onorm = beta;
120: if (norm) {
121: if (PetscUnlikely(bv->indef)) PetscCall(BV_NormVecOrColumn(bv,j,v,norm));
122: else {
123: /* estimate |v'| from |v| */
124: PetscCall(BV_SquareSum(bv,j,c,&sum));
125: *norm = beta*beta-sum;
126: if (PetscUnlikely(*norm <= 0.0)) PetscCall(BV_NormVecOrColumn(bv,j,v,norm));
127: else *norm = PetscSqrtReal(*norm);
128: }
129: }
130: PetscCall(BV_AddCoefficients(bv,j,h,c));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
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)
136: /*
137: BVOrthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
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: static PetscErrorCode BVOrthogonalizeGS(BV bv,PetscInt j,Vec v,PetscBool *which,PetscReal *norm,PetscBool *lindep)
146: {
147: PetscScalar *h,*c,*omega;
148: PetscReal onrm,nrm;
149: PetscInt k,l;
150: PetscBool mgs,dolindep,signature;
152: PetscFunctionBegin;
153: if (v) {
154: k = bv->k;
155: h = bv->h;
156: c = bv->c;
157: } else {
158: k = j;
159: h = NULL;
160: c = NULL;
161: }
163: mgs = (bv->orthog_type==BV_ORTHOG_MGS)? PETSC_TRUE: PETSC_FALSE;
165: /* if indefinite inner product, skip the computation of lindep */
166: if (bv->indef && lindep) *lindep = PETSC_FALSE;
167: dolindep = (!bv->indef && lindep)? PETSC_TRUE: PETSC_FALSE;
169: /* if indefinite and we are orthogonalizing a column, the norm must always be computed */
170: signature = (bv->indef && !v)? PETSC_TRUE: PETSC_FALSE;
172: PetscCall(BV_CleanCoefficients(bv,k,h));
174: switch (bv->orthog_ref) {
176: case BV_ORTHOG_REFINE_IFNEEDED:
177: PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,&onrm,&nrm));
178: /* repeat if ||q|| < eta ||h|| */
179: l = 1;
180: while (l<3 && nrm && PetscAbsReal(nrm) < bv->orthog_eta*PetscAbsReal(onrm)) {
181: l++;
182: if (mgs||bv->indef) onrm = nrm;
183: 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: if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
187: break;
189: case BV_ORTHOG_REFINE_NEVER:
190: PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL));
191: /* compute ||v|| */
192: if (norm || dolindep || signature) PetscCall(BV_NormVecOrColumn(bv,k,v,&nrm));
193: /* linear dependence check: just test for exactly zero norm */
194: if (dolindep) *lindep = PetscNot(nrm);
195: break;
197: case BV_ORTHOG_REFINE_ALWAYS:
198: PetscCall(BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL));
199: 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: if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
202: break;
203: }
204: if (signature) {
205: PetscCall(VecGetArray(bv->omega,&omega));
206: omega[bv->nc+k] = (nrm<0.0)? -1.0: 1.0;
207: PetscCall(VecRestoreArray(bv->omega,&omega));
208: }
209: if (norm) {
210: *norm = nrm;
211: if (!v) { /* store norm value next to the orthogonalization coefficients */
212: if (dolindep && *lindep) PetscCall(BV_SetValue(bv,k,k,h,0.0));
213: else PetscCall(BV_SetValue(bv,k,k,h,nrm));
214: }
215: }
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: BVOrthogonalizeVec - Orthogonalize a given vector with respect to all
221: active columns.
223: Collective
225: Input Parameters:
226: + bv - the basis vectors context
227: - v - the vector
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
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.
242: In the case of an indefinite inner product, the lindep parameter is not
243: computed (set to false).
245: Level: advanced
247: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVSetActiveColumns(), BVGetNumConstraints()
248: @*/
249: PetscErrorCode BVOrthogonalizeVec(BV bv,Vec v,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
250: {
251: PetscInt ksave,lsave;
253: PetscFunctionBegin;
257: BVCheckSizes(bv,1);
259: PetscCheckSameComm(bv,1,v,2);
261: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
262: ksave = bv->k;
263: lsave = bv->l;
264: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
265: PetscCall(BV_AllocateCoeffs(bv));
266: PetscCall(BV_AllocateSignature(bv));
267: PetscCall(BVOrthogonalizeGS(bv,0,v,NULL,norm,lindep));
268: bv->k = ksave;
269: bv->l = lsave;
270: if (H) PetscCall(BV_StoreCoefficients(bv,bv->k,bv->h,H));
271: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: BVOrthogonalizeColumn - Orthogonalize one of the column vectors with respect to
277: the previous ones.
279: Collective
281: Input Parameters:
282: + bv - the basis vectors context
283: - j - index of column to be orthogonalized
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
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.
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).
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.
306: This routine does not normalize the resulting vector, see BVOrthonormalizeColumn().
308: In the case of an indefinite inner product, the lindep parameter is not
309: computed (set to false).
311: Level: advanced
313: .seealso: BVSetOrthogonalization(), BVSetMatrix(), BVSetActiveColumns(), BVOrthogonalize(), BVOrthogonalizeVec(), BVGetNumConstraints(), BVOrthonormalizeColumn()
314: @*/
315: PetscErrorCode BVOrthogonalizeColumn(BV bv,PetscInt j,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
316: {
317: PetscInt ksave,lsave;
319: PetscFunctionBegin;
323: BVCheckSizes(bv,1);
324: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
325: 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);
327: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
328: ksave = bv->k;
329: lsave = bv->l;
330: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
331: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
332: PetscCall(BV_AllocateSignature(bv));
333: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,norm,lindep));
334: bv->k = ksave;
335: bv->l = lsave;
336: if (H) PetscCall(BV_StoreCoefficients(bv,j,NULL,H));
337: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
338: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@
343: BVOrthonormalizeColumn - Orthonormalize one of the column vectors with respect to
344: the previous ones.
346: Collective
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
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
358: Notes:
359: This is equivalent to a call to BVOrthogonalizeColumn() followed by a
360: call to BVScaleColumn() with the reciprocal of the norm.
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.
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.
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.
376: Level: advanced
378: .seealso: BVOrthogonalizeColumn(), BVScaleColumn()
379: @*/
380: PetscErrorCode BVOrthonormalizeColumn(BV bv,PetscInt j,PetscBool replace,PetscReal *norm,PetscBool *lindep)
381: {
382: PetscScalar alpha;
383: PetscReal nrm;
384: PetscInt ksave,lsave;
385: PetscBool lndep;
387: PetscFunctionBegin;
391: BVCheckSizes(bv,1);
392: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
393: 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);
395: /* orthogonalize */
396: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
397: ksave = bv->k;
398: lsave = bv->l;
399: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
400: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
401: PetscCall(BV_AllocateSignature(bv));
402: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
403: if (replace && (nrm==0.0 || lndep)) {
404: PetscCall(PetscInfo(bv,"Vector was linearly dependent, generating a new random vector\n"));
405: PetscCall(BVSetRandomColumn(bv,j));
406: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
407: if (nrm==0.0 || lndep) { /* yet another attempt */
408: PetscCall(BVSetRandomColumn(bv,j));
409: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
410: }
411: }
412: bv->k = ksave;
413: bv->l = lsave;
414: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
416: /* scale */
417: if (nrm!=1.0 && nrm!=0.0) {
418: alpha = 1.0/nrm;
419: PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
420: PetscUseTypeMethod(bv,scale,j,alpha);
421: PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
422: }
423: if (norm) *norm = nrm;
424: if (lindep) *lindep = lndep;
425: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: /*@
430: BVOrthogonalizeSomeColumn - Orthogonalize one of the column vectors with
431: respect to some of the previous ones.
433: Collective
435: Input Parameters:
436: + bv - the basis vectors context
437: . j - index of column to be orthogonalized
438: - which - logical array indicating selected columns
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
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.
451: The use of this operation is restricted to MGS orthogonalization type.
453: In the case of an indefinite inner product, the lindep parameter is not
454: computed (set to false).
456: Level: advanced
458: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization()
459: @*/
460: PetscErrorCode BVOrthogonalizeSomeColumn(BV bv,PetscInt j,PetscBool *which,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
461: {
462: PetscInt ksave,lsave;
464: PetscFunctionBegin;
467: PetscAssertPointer(which,3);
469: BVCheckSizes(bv,1);
470: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
471: 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: PetscCheck(bv->orthog_type==BV_ORTHOG_MGS,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation only available for MGS orthogonalization");
474: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
475: ksave = bv->k;
476: lsave = bv->l;
477: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
478: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
479: PetscCall(BV_AllocateSignature(bv));
480: PetscCall(BVOrthogonalizeGS(bv,j,NULL,which,norm,lindep));
481: bv->k = ksave;
482: bv->l = lsave;
483: if (H) PetscCall(BV_StoreCoefficients(bv,j,NULL,H));
484: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
485: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /*
490: Block Gram-Schmidt: V2 = V2 - V1*R12, where R12 = V1'*V2
491: */
492: static PetscErrorCode BVOrthogonalize_BlockGS(BV V,Mat R)
493: {
494: BV V1;
496: PetscFunctionBegin;
497: PetscCall(BVGetSplit(V,&V1,NULL));
498: PetscCall(BVDot(V,V1,R));
499: PetscCall(BVMult(V,-1.0,1.0,V1,R));
500: PetscCall(BVRestoreSplit(V,&V1,NULL));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: /*
505: Orthogonalize a set of vectors with Gram-Schmidt, column by column.
506: */
507: static PetscErrorCode BVOrthogonalize_GS(BV V,Mat R)
508: {
509: PetscScalar *r=NULL;
510: PetscReal norm;
511: PetscInt j,ldr,lsave;
512: Vec v,w;
514: PetscFunctionBegin;
515: if (R) {
516: PetscCall(MatDenseGetLDA(R,&ldr));
517: PetscCall(MatDenseGetArray(R,&r));
518: }
519: if (V->matrix) {
520: PetscCall(BVGetCachedBV(V,&V->cached));
521: PetscCall(BVSetActiveColumns(V->cached,V->l,V->k));
522: }
523: for (j=V->l;j<V->k;j++) {
524: if (V->matrix && V->orthog_type==BV_ORTHOG_MGS) { /* fill cached BV */
525: PetscCall(BVGetColumn(V->cached,j,&v));
526: PetscCall(BVGetColumn(V,j,&w));
527: PetscCall(MatMult(V->matrix,w,v));
528: PetscCall(BVRestoreColumn(V,j,&w));
529: PetscCall(BVRestoreColumn(V->cached,j,&v));
530: }
531: if (R) {
532: PetscCall(BVOrthogonalizeColumn(V,j,NULL,&norm,NULL));
533: lsave = V->l;
534: V->l = -V->nc;
535: PetscCall(BV_StoreCoefficients(V,j,NULL,r+j*ldr));
536: V->l = lsave;
537: r[j+j*ldr] = norm;
538: } else PetscCall(BVOrthogonalizeColumn(V,j,NULL,&norm,NULL));
539: PetscCheck(norm,PetscObjectComm((PetscObject)V),PETSC_ERR_CONV_FAILED,"Breakdown in BVOrthogonalize due to a linearly dependent column");
540: if (V->matrix && V->orthog_type==BV_ORTHOG_CGS) { /* fill cached BV */
541: PetscCall(BVGetColumn(V->cached,j,&v));
542: PetscCall(VecCopy(V->Bx,v));
543: PetscCall(BVRestoreColumn(V->cached,j,&v));
544: }
545: PetscCall(BVScaleColumn(V,j,1.0/norm));
546: }
547: if (R) PetscCall(MatDenseRestoreArray(R,&r));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: /*
552: BV_GetBufferMat - Create auxiliary seqdense matrix that wraps the bv->buffer.
553: */
554: static inline PetscErrorCode BV_GetBufferMat(BV bv)
555: {
556: PetscInt ld;
557: PetscScalar *array;
559: PetscFunctionBegin;
560: if (!bv->Abuffer) {
561: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
562: ld = bv->m+bv->nc;
563: PetscCall(VecGetArray(bv->buffer,&array));
564: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,bv->m,array,&bv->Abuffer));
565: PetscCall(VecRestoreArray(bv->buffer,&array));
566: }
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
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: static inline PetscErrorCode BV_StoreCoeffsBlock_Default(BV bv,Mat R,PetscBool tri)
576: {
577: const PetscScalar *bb;
578: PetscScalar *rr;
579: PetscInt j,ldr,ldb;
581: PetscFunctionBegin;
582: PetscCall(MatDenseGetLDA(R,&ldr));
583: PetscCall(MatDenseGetArray(R,&rr));
584: ldb = bv->m+bv->nc;
585: PetscCall(VecGetArrayRead(bv->buffer,&bb));
586: for (j=bv->l;j<bv->k;j++) PetscCall(PetscArraycpy(rr+j*ldr,bb+j*ldb,(tri?(j+1):bv->k)+bv->nc));
587: PetscCall(VecRestoreArrayRead(bv->buffer,&bb));
588: PetscCall(MatDenseRestoreArray(R,&rr));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*
593: Orthogonalize a set of vectors with Cholesky: R=chol(V'*V), Q=V*inv(R)
594: */
595: static PetscErrorCode BVOrthogonalize_Chol(BV V,Mat Rin)
596: {
597: Mat R,S;
599: PetscFunctionBegin;
600: PetscCall(BV_GetBufferMat(V));
601: R = V->Abuffer;
602: if (Rin) S = Rin; /* use Rin as a workspace for S */
603: else S = R;
604: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
605: PetscCall(BVDot(V,V,R));
606: PetscCall(BVMatCholInv_LAPACK_Private(V,R,S));
607: PetscCall(BVMultInPlace(V,S,V->l,V->k));
608: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: /*
613: Orthogonalize a set of vectors with the Tall-Skinny QR method
614: */
615: static PetscErrorCode BVOrthogonalize_TSQR(BV V,Mat Rin)
616: {
617: PetscScalar *pv,*r=NULL;
618: PetscInt ldr;
619: Mat R;
621: PetscFunctionBegin;
622: PetscCall(BV_GetBufferMat(V));
623: R = V->Abuffer;
624: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
625: PetscCall(MatDenseGetLDA(R,&ldr));
626: PetscCall(MatDenseGetArray(R,&r));
627: PetscCall(BVGetArray(V,&pv));
628: 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: PetscCall(BVRestoreArray(V,&pv));
630: PetscCall(MatDenseRestoreArray(R,&r));
631: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /*
636: Orthogonalize a set of vectors with TSQR, but computing R only, then doing Q=V*inv(R)
637: */
638: static PetscErrorCode BVOrthogonalize_TSQRCHOL(BV V,Mat Rin)
639: {
640: PetscScalar *pv,*r=NULL;
641: PetscInt ldr;
642: Mat R,S;
644: PetscFunctionBegin;
645: PetscCall(BV_GetBufferMat(V));
646: R = V->Abuffer;
647: if (Rin) S = Rin; /* use Rin as a workspace for S */
648: else S = R;
649: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
650: PetscCall(MatDenseGetLDA(R,&ldr));
651: PetscCall(MatDenseGetArray(R,&r));
652: PetscCall(BVGetArray(V,&pv));
653: 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: PetscCall(BVRestoreArray(V,&pv));
655: PetscCall(MatDenseRestoreArray(R,&r));
656: PetscCall(BVMatTriInv_LAPACK_Private(V,R,S));
657: PetscCall(BVMultInPlace(V,S,V->l,V->k));
658: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*
663: Orthogonalize a set of vectors with SVQB
664: */
665: static PetscErrorCode BVOrthogonalize_SVQB(BV V,Mat Rin)
666: {
667: Mat R,S;
669: PetscFunctionBegin;
670: PetscCall(BV_GetBufferMat(V));
671: R = V->Abuffer;
672: if (Rin) S = Rin; /* use Rin as a workspace for S */
673: else S = R;
674: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
675: PetscCall(BVDot(V,V,R));
676: PetscCall(BVMatSVQB_LAPACK_Private(V,R,S));
677: PetscCall(BVMultInPlace(V,S,V->l,V->k));
678: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_FALSE));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: /*@
683: BVOrthogonalize - Orthogonalize all columns (starting from the leading ones),
684: that is, compute the QR decomposition.
686: Collective
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
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()).
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).
709: Can pass NULL if R is not required.
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.
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.
722: Level: intermediate
724: .seealso: BVOrthogonalizeColumn(), BVOrthogonalizeVec(), BVSetMatrix(), BVSetActiveColumns(), BVSetOrthogonalization(), BVOrthogBlockType
725: @*/
726: PetscErrorCode BVOrthogonalize(BV V,Mat R)
727: {
728: PetscInt m,n;
730: PetscFunctionBegin;
733: BVCheckSizes(V,1);
734: if (R) {
737: PetscCheckTypeName(R,MATSEQDENSE);
738: PetscCall(MatGetSize(R,&m,&n));
739: 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: 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: PetscCheck(!V->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Not implemented for BV with constraints, use BVOrthogonalizeColumn() instead");
744: PetscCall(PetscLogEventBegin(BV_Orthogonalize,V,R,0,0));
745: switch (V->orthog_block) {
746: case BV_ORTHOG_BLOCK_GS: /* proceed column by column with Gram-Schmidt */
747: PetscCall(BVOrthogonalize_GS(V,R));
748: break;
749: case BV_ORTHOG_BLOCK_CHOL:
750: PetscCall(BVOrthogonalize_Chol(V,R));
751: break;
752: case BV_ORTHOG_BLOCK_TSQR:
753: PetscCheck(!V->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
754: PetscCall(BVOrthogonalize_TSQR(V,R));
755: break;
756: case BV_ORTHOG_BLOCK_TSQRCHOL:
757: PetscCheck(!V->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
758: PetscCall(BVOrthogonalize_TSQRCHOL(V,R));
759: break;
760: case BV_ORTHOG_BLOCK_SVQB:
761: PetscCall(BVOrthogonalize_SVQB(V,R));
762: break;
763: }
764: PetscCall(PetscLogEventEnd(BV_Orthogonalize,V,R,0,0));
765: PetscCall(PetscObjectStateIncrease((PetscObject)V));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }