Actual source code: bvorthog.c

slepc-3.21.1 2024-04-26
Report Typos and Errors
```  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: }
```