Actual source code: bvorthog.c
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: [](sec:bv), `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: The columns $v_{0:j-1}$ are assumed to be mutually orthonormal.
296: Leading columns $v_{0:\ell-1}$ also participate in the orthogonalization, as well
297: as the constraints. If `H` is given, it must have enough space to store
298: $j-\ell+1$ coefficients (the last coefficient will contain the value `norm`, unless
299: the `norm` argument is `NULL`).
301: If a non-standard inner product has been specified with `BVSetMatrix()`,
302: then the vector is B-orthogonalized, using the non-standard inner product
303: defined by matrix $B$. The output vector satisfies $v_j^*Bv_{0:j-1} = 0$.
305: This routine does not normalize the resulting vector, see `BVOrthonormalizeColumn()`.
307: In the case of an indefinite inner product, the `lindep` parameter is not
308: computed (set to false).
310: Level: advanced
312: .seealso: [](sec:bv), `BVSetOrthogonalization()`, `BVSetMatrix()`, `BVSetActiveColumns()`, `BVOrthogonalize()`, `BVOrthogonalizeVec()`, `BVGetNumConstraints()`, `BVOrthonormalizeColumn()`
313: @*/
314: PetscErrorCode BVOrthogonalizeColumn(BV bv,PetscInt j,PetscScalar H[],PetscReal *norm,PetscBool *lindep)
315: {
316: PetscInt ksave,lsave;
318: PetscFunctionBegin;
322: BVCheckSizes(bv,1);
323: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
324: 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: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
327: ksave = bv->k;
328: lsave = bv->l;
329: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
330: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
331: PetscCall(BV_AllocateSignature(bv));
332: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,norm,lindep));
333: bv->k = ksave;
334: bv->l = lsave;
335: if (H) PetscCall(BV_StoreCoefficients(bv,j,NULL,H));
336: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
337: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: BVOrthonormalizeColumn - Orthonormalize one of the column vectors with respect to
343: the previous ones.
345: Collective
347: Input Parameters:
348: + bv - the basis vectors context
349: . j - index of column to be orthonormalized
350: - replace - whether it is allowed to set the vector randomly
352: Output Parameters:
353: + norm - (optional) norm of the vector after orthogonalization and before normalization
354: - lindep - (optional) flag indicating that linear dependence was determined during
355: orthogonalization
357: Notes:
358: This is equivalent to a call to `BVOrthogonalizeColumn()` followed by a
359: call to `BVScaleColumn()` with the reciprocal of the norm.
361: This function first orthogonalizes vector $v_j$ with respect to $v_{0:j-1}$.
362: A byproduct of this computation is `norm`, the norm of the vector after
363: orthogonalization. Secondly, it scales the vector with `1/norm`, so that the
364: resulting vector has unit norm.
366: If after orthogonalization the vector $v_j$ is exactly zero, it cannot be normalized
367: because `norm=0`. In that case, it could be left as zero or replaced by a random
368: vector that is then orthonormalized. The latter is achieved by setting the
369: argument `replace` to `PETSC_TRUE`. The vector will be replaced by a random vector also
370: if `lindep` was set to `PETSC_TRUE`, even if the norm is not exactly zero.
372: If the vector has been replaced by a random vector, the output arguments `norm` and
373: `lindep` will be set according to the orthogonalization of this new vector.
375: Level: advanced
377: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVScaleColumn()`
378: @*/
379: PetscErrorCode BVOrthonormalizeColumn(BV bv,PetscInt j,PetscBool replace,PetscReal *norm,PetscBool *lindep)
380: {
381: PetscScalar alpha;
382: PetscReal nrm;
383: PetscInt ksave,lsave;
384: PetscBool lndep;
386: PetscFunctionBegin;
390: BVCheckSizes(bv,1);
391: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
392: 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: /* orthogonalize */
395: PetscCall(PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0));
396: ksave = bv->k;
397: lsave = bv->l;
398: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
399: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
400: PetscCall(BV_AllocateSignature(bv));
401: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
402: if (replace && (nrm==0.0 || lndep)) {
403: PetscCall(PetscInfo(bv,"Vector was linearly dependent, generating a new random vector\n"));
404: PetscCall(BVSetRandomColumn(bv,j));
405: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
406: if (nrm==0.0 || lndep) { /* yet another attempt */
407: PetscCall(BVSetRandomColumn(bv,j));
408: PetscCall(BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep));
409: }
410: }
411: bv->k = ksave;
412: bv->l = lsave;
413: PetscCall(PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0));
415: /* scale */
416: if (nrm!=1.0 && nrm!=0.0) {
417: alpha = 1.0/nrm;
418: PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
419: PetscUseTypeMethod(bv,scale,j,alpha);
420: PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
421: }
422: if (norm) *norm = nrm;
423: if (lindep) *lindep = lndep;
424: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: /*@
429: BVOrthogonalizeSomeColumn - Orthogonalize one of the column vectors with
430: respect to some of the previous ones.
432: Collective
434: Input Parameters:
435: + bv - the basis vectors context
436: . j - index of column to be orthogonalized
437: - which - logical array indicating selected columns
439: Output Parameters:
440: + H - (optional) coefficients computed during orthogonalization
441: . norm - (optional) norm of the vector after being orthogonalized
442: - lindep - (optional) flag indicating that refinement did not improve the quality
443: of orthogonalization
445: Notes:
446: This function is similar to `BVOrthogonalizeColumn()`, but $v_j$ is
447: orthogonalized only against columns $v_i$ having `which[i]=PETSC_TRUE`.
448: The length of array `which` must be `j` at least.
450: The use of this operation is restricted to MGS orthogonalization type,
451: `BV_ORTHOG_MGS`.
453: In the case of an indefinite inner product, the `lindep` parameter is not
454: computed (set to false).
456: Level: advanced
458: .seealso: [](sec:bv), `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: /*MC
505: BV_ORTHOG_BLOCK_GS - Orthogonalize a set of vectors with Gram-Schmidt, column by column.
507: Note:
508: In this case, the block orthogonalization is done in the same way as the
509: orthogonalization of individual vectors, e.g., `BVOrthogonalizeColumn`, and
510: hence the settings provided via `BVSetOrthogonalization()` apply.
512: Level: advanced
514: .seealso: [](sec:bv), `BVOrthogBlockType`, `BVSetOrthogonalization()`, `BVOrthogonalize()`, `BV_ORTHOG_BLOCK_CHOL`, `BV_ORTHOG_BLOCK_TSQR`, `BV_ORTHOG_BLOCK_TSQRCHOL`, `BV_ORTHOG_BLOCK_SVQB`
515: M*/
516: static PetscErrorCode BVOrthogonalize_GS(BV V,Mat R)
517: {
518: PetscScalar *r=NULL;
519: PetscReal norm;
520: PetscInt j,ldr,lsave;
521: Vec v,w;
523: PetscFunctionBegin;
524: if (R) {
525: PetscCall(MatDenseGetLDA(R,&ldr));
526: PetscCall(MatDenseGetArray(R,&r));
527: }
528: if (V->matrix) {
529: PetscCall(BVGetCachedBV(V,&V->cached));
530: PetscCall(BVSetActiveColumns(V->cached,V->l,V->k));
531: }
532: for (j=V->l;j<V->k;j++) {
533: if (V->matrix && V->orthog_type==BV_ORTHOG_MGS) { /* fill cached BV */
534: PetscCall(BVGetColumn(V->cached,j,&v));
535: PetscCall(BVGetColumn(V,j,&w));
536: PetscCall(MatMult(V->matrix,w,v));
537: PetscCall(BVRestoreColumn(V,j,&w));
538: PetscCall(BVRestoreColumn(V->cached,j,&v));
539: }
540: if (R) {
541: PetscCall(BVOrthogonalizeColumn(V,j,NULL,&norm,NULL));
542: lsave = V->l;
543: V->l = -V->nc;
544: PetscCall(BV_StoreCoefficients(V,j,NULL,r+j*ldr));
545: V->l = lsave;
546: r[j+j*ldr] = norm;
547: } else PetscCall(BVOrthogonalizeColumn(V,j,NULL,&norm,NULL));
548: PetscCheck(norm,PetscObjectComm((PetscObject)V),PETSC_ERR_CONV_FAILED,"Breakdown in BVOrthogonalize due to a linearly dependent column");
549: if (V->matrix && V->orthog_type==BV_ORTHOG_CGS) { /* fill cached BV */
550: PetscCall(BVGetColumn(V->cached,j,&v));
551: PetscCall(VecCopy(V->Bx,v));
552: PetscCall(BVRestoreColumn(V->cached,j,&v));
553: }
554: PetscCall(BVScaleColumn(V,j,1.0/norm));
555: }
556: if (R) PetscCall(MatDenseRestoreArray(R,&r));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: /*
561: BV_GetBufferMat - Create auxiliary seqdense matrix that wraps the bv->buffer.
562: */
563: static inline PetscErrorCode BV_GetBufferMat(BV bv)
564: {
565: PetscInt ld;
566: PetscScalar *array;
568: PetscFunctionBegin;
569: if (!bv->Abuffer) {
570: if (!bv->buffer) PetscCall(BVGetBufferVec(bv,&bv->buffer));
571: ld = bv->m+bv->nc;
572: PetscCall(VecGetArray(bv->buffer,&array));
573: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,bv->m,array,&bv->Abuffer));
574: PetscCall(VecRestoreArray(bv->buffer,&array));
575: }
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: /*
580: BV_StoreCoeffsBlock_Default - Copy the contents of the BV buffer to a dense Mat
581: provided by the caller. Only columns l:k-1 are copied, restricting to the upper
582: triangular part if tri=PETSC_TRUE.
583: */
584: static inline PetscErrorCode BV_StoreCoeffsBlock_Default(BV bv,Mat R,PetscBool tri)
585: {
586: const PetscScalar *bb;
587: PetscScalar *rr;
588: PetscInt j,ldr,ldb;
590: PetscFunctionBegin;
591: PetscCall(MatDenseGetLDA(R,&ldr));
592: PetscCall(MatDenseGetArray(R,&rr));
593: ldb = bv->m+bv->nc;
594: PetscCall(VecGetArrayRead(bv->buffer,&bb));
595: for (j=bv->l;j<bv->k;j++) PetscCall(PetscArraycpy(rr+j*ldr,bb+j*ldb,(tri?(j+1):bv->k)+bv->nc));
596: PetscCall(VecRestoreArrayRead(bv->buffer,&bb));
597: PetscCall(MatDenseRestoreArray(R,&rr));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*MC
602: BV_ORTHOG_BLOCK_CHOL - Orthogonalize a set of vectors $V$ with the Cholesky QR
603: method, that is, compute $R=\mathrm{chol}(V^*V)$, then obtain the
604: orthogonalized vectors as $Q=V R^{-1}$.
606: Developer Note:
607: This is plain Cholesky QR, which may fail in case of ill-conditioned $V$.
608: An iterated version such as CholeskyQR2 should be preferred {cite:p}`Yam15`.
610: Level: advanced
612: .seealso: [](sec:bv), `BVOrthogBlockType`, `BVSetOrthogonalization()`, `BVOrthogonalize()`, `BV_ORTHOG_BLOCK_GS`, `BV_ORTHOG_BLOCK_TSQR`, `BV_ORTHOG_BLOCK_TSQRCHOL`, `BV_ORTHOG_BLOCK_SVQB`
613: M*/
614: static PetscErrorCode BVOrthogonalize_Chol(BV V,Mat Rin)
615: {
616: Mat R,S;
618: PetscFunctionBegin;
619: PetscCall(BV_GetBufferMat(V));
620: R = V->Abuffer;
621: if (Rin) S = Rin; /* use Rin as a workspace for S */
622: else S = R;
623: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
624: PetscCall(BVDot(V,V,R));
625: PetscCall(BVMatCholInv_LAPACK_Private(V,R,S));
626: PetscCall(BVMultInPlace(V,S,V->l,V->k));
627: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /*MC
632: BV_ORTHOG_BLOCK_TSQR - Orthogonalize a set of vectors with the Tall-Skinny QR method {cite:p}`Dem12`.
634: Level: advanced
636: .seealso: [](sec:bv), `BVOrthogBlockType`, `BVSetOrthogonalization()`, `BVOrthogonalize()`, `BV_ORTHOG_BLOCK_GS`, `BV_ORTHOG_BLOCK_CHOL`, `BV_ORTHOG_BLOCK_TSQRCHOL`, `BV_ORTHOG_BLOCK_SVQB`
637: M*/
638: static PetscErrorCode BVOrthogonalize_TSQR(BV V,Mat Rin)
639: {
640: PetscScalar *pv,*r=NULL;
641: PetscInt ldr;
642: Mat R;
644: PetscFunctionBegin;
645: PetscCall(BV_GetBufferMat(V));
646: R = V->Abuffer;
647: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
648: PetscCall(MatDenseGetLDA(R,&ldr));
649: PetscCall(MatDenseGetArray(R,&r));
650: PetscCall(BVGetArray(V,&pv));
651: 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));
652: PetscCall(BVRestoreArray(V,&pv));
653: PetscCall(MatDenseRestoreArray(R,&r));
654: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*MC
659: BV_ORTHOG_BLOCK_TSQRCHOL - Orthogonalize a set of vectors $V$ with the TSQR
660: method, but computing the triangular factor $R$ only, then retrieving the
661: orthogonalized vectors as $Q=V R^{-1}$.
663: Level: advanced
665: .seealso: [](sec:bv), `BVOrthogBlockType`, `BVSetOrthogonalization()`, `BVOrthogonalize()`, `BV_ORTHOG_BLOCK_GS`, `BV_ORTHOG_BLOCK_CHOL`, `BV_ORTHOG_BLOCK_TSQR`, `BV_ORTHOG_BLOCK_SVQB`
666: M*/
667: static PetscErrorCode BVOrthogonalize_TSQRCHOL(BV V,Mat Rin)
668: {
669: PetscScalar *pv,*r=NULL;
670: PetscInt ldr;
671: Mat R,S;
673: PetscFunctionBegin;
674: PetscCall(BV_GetBufferMat(V));
675: R = V->Abuffer;
676: if (Rin) S = Rin; /* use Rin as a workspace for S */
677: else S = R;
678: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
679: PetscCall(MatDenseGetLDA(R,&ldr));
680: PetscCall(MatDenseGetArray(R,&r));
681: PetscCall(BVGetArray(V,&pv));
682: 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));
683: PetscCall(BVRestoreArray(V,&pv));
684: PetscCall(MatDenseRestoreArray(R,&r));
685: PetscCall(BVMatTriInv_LAPACK_Private(V,R,S));
686: PetscCall(BVMultInPlace(V,S,V->l,V->k));
687: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: /*MC
692: BV_ORTHOG_BLOCK_SVQB - Orthogonalize a set of vectors with the SVQB method {cite:p}`Sta02`.
694: Note:
695: As opposed to other block orthogonalization methods, SVQB does not compute a
696: QR decomposition but a QB decomposition, i.e., the second factor is not
697: triangular.
699: Level: advanced
701: .seealso: [](sec:bv), `BVOrthogBlockType`, `BVSetOrthogonalization()`, `BVOrthogonalize()`, `BV_ORTHOG_BLOCK_GS`, `BV_ORTHOG_BLOCK_CHOL`, `BV_ORTHOG_BLOCK_TSQR`, `BV_ORTHOG_BLOCK_TSQRCHOL`
702: M*/
703: static PetscErrorCode BVOrthogonalize_SVQB(BV V,Mat Rin)
704: {
705: Mat R,S;
707: PetscFunctionBegin;
708: PetscCall(BV_GetBufferMat(V));
709: R = V->Abuffer;
710: if (Rin) S = Rin; /* use Rin as a workspace for S */
711: else S = R;
712: if (V->l) PetscCall(BVOrthogonalize_BlockGS(V,R));
713: PetscCall(BVDot(V,V,R));
714: PetscCall(BVMatSVQB_LAPACK_Private(V,R,S));
715: PetscCall(BVMultInPlace(V,S,V->l,V->k));
716: if (Rin) PetscCall(BV_StoreCoeffsBlock_Default(V,Rin,PETSC_FALSE));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: BVOrthogonalize - Orthogonalize all columns (starting from the leading ones),
722: that is, compute the QR decomposition.
724: Collective
726: Input Parameters:
727: + V - basis vectors to be orthogonalized (or $B$-orthogonalized), modified on output
728: - R - a sequential dense matrix (or `NULL`), on output the triangular factor of
729: the QR decomposition
731: Notes:
732: On input, matrix `R` must be a square sequential dense `Mat`, with at least as many
733: rows and columns as the number of active columns of `V`. The output satisfies
734: $V = \tilde V R$ (where $\tilde V$ represents `V` on exit) and $\tilde V^*\tilde V = I$
735: (or $\tilde V^*B\tilde V = I$ if an inner product matrix $B$ has been specified with
736: `BVSetMatrix()`).
738: If `V` has leading columns, then they are not modified (are assumed to be already
739: orthonormal) and the leading columns of `R` are not referenced. Let the
740: decomposition be
741: $$
742: \begin{bmatrix}
743: V_1 & V_2
744: \end{bmatrix}=
745: \begin{bmatrix}
746: \tilde V_1 & \tilde V_2
747: \end{bmatrix}
748: \begin{bmatrix}
749: R_{11} & R_{12} \\
750: 0 & R_{22}
751: \end{bmatrix}
752: $$
753: then $\tilde V_1$ is left unchanged (equal to $V_1$) as well as $R_{11}$ (they should
754: satisfy $V_1 = \tilde V_1 R_{11}$).
756: Can pass `NULL` if `R` is not required.
758: The method to be used for block orthogonalization can be set with
759: `BVSetOrthogonalization()`. If set to `BV_ORTHOG_BLOCK_GS`, the computation is done
760: column by column with successive calls to `BVOrthogonalizeColumn()`. Note that in the
761: `BV_ORTHOG_BLOCK_SVQB` method the `R` factor is not upper triangular.
763: If `V` is rank-deficient or very ill-conditioned, that is, one or more columns are
764: (almost) linearly dependent with respect to the rest, then the algorithm may
765: break down or result in larger numerical error. Linearly dependent columns are
766: essentially replaced by random directions, and the corresponding diagonal entry
767: in `R` is set to (nearly) zero.
769: Level: intermediate
771: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVOrthogonalizeVec()`, `BVSetMatrix()`, `BVSetActiveColumns()`, `BVSetOrthogonalization()`, `BVOrthogBlockType`
772: @*/
773: PetscErrorCode BVOrthogonalize(BV V,Mat R)
774: {
775: PetscInt m,n;
777: PetscFunctionBegin;
780: BVCheckSizes(V,1);
781: if (R) {
784: PetscCheckTypeName(R,MATSEQDENSE);
785: PetscCall(MatGetSize(R,&m,&n));
786: 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);
787: 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);
788: }
789: PetscCheck(!V->nc,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Not implemented for BV with constraints, use BVOrthogonalizeColumn() instead");
791: PetscCall(PetscLogEventBegin(BV_Orthogonalize,V,R,0,0));
792: switch (V->orthog_block) {
793: case BV_ORTHOG_BLOCK_GS: /* proceed column by column with Gram-Schmidt */
794: PetscCall(BVOrthogonalize_GS(V,R));
795: break;
796: case BV_ORTHOG_BLOCK_CHOL:
797: PetscCall(BVOrthogonalize_Chol(V,R));
798: break;
799: case BV_ORTHOG_BLOCK_TSQR:
800: PetscCheck(!V->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
801: PetscCall(BVOrthogonalize_TSQR(V,R));
802: break;
803: case BV_ORTHOG_BLOCK_TSQRCHOL:
804: PetscCheck(!V->matrix,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
805: PetscCall(BVOrthogonalize_TSQRCHOL(V,R));
806: break;
807: case BV_ORTHOG_BLOCK_SVQB:
808: PetscCall(BVOrthogonalize_SVQB(V,R));
809: break;
810: }
811: PetscCall(PetscLogEventEnd(BV_Orthogonalize,V,R,0,0));
812: PetscCall(PetscObjectStateIncrease((PetscObject)V));
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }