Actual source code: bvorthog.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>          /*I   "slepcbv.h"   I*/

 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: PETSC_STATIC_INLINE PetscErrorCode BV_NormVecOrColumn(BV bv,PetscInt j,Vec v,PetscReal *nrm)
 21: {

 25:   if (v) { BVNormVec(bv,v,NORM_2,nrm); }
 26:   else { BVNormColumn(bv,j,NORM_2,nrm); }
 27:   return(0);
 28: }

 30: /*
 31:    BVDotColumnInc - Same as BVDotColumn() but also including column j, which
 32:    is multiplied by itself
 33: */
 34: PETSC_STATIC_INLINE PetscErrorCode BVDotColumnInc(BV X,PetscInt j,PetscScalar *q)
 35: {
 37:   PetscInt       ksave;
 38:   Vec            y;

 41:   PetscLogEventBegin(BV_DotVec,X,0,0,0);
 42:   ksave = X->k;
 43:   X->k = j+1;
 44:   BVGetColumn(X,j,&y);
 45:   (*X->ops->dotvec)(X,y,q);
 46:   BVRestoreColumn(X,j,&y);
 47:   X->k = ksave;
 48:   PetscLogEventEnd(BV_DotVec,X,0,0,0);
 49:   return(0);
 50: }

 52: /*
 53:    BVOrthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt
 54: */
 55: static PetscErrorCode BVOrthogonalizeMGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onrm,PetscReal *nrm)
 56: {
 58:   PetscInt          i;
 59:   PetscScalar       dot;
 60:   Vec               vi,z,w=v;
 61:   const PetscScalar *omega;

 64:   if (!v) { BVGetColumn(bv,j,&w); }
 65:   if (onrm) { BVNormVec(bv,w,NORM_2,onrm); }
 66:   z = w;
 67:   if (bv->indef) {
 68:     VecGetArrayRead(bv->omega,&omega);
 69:   }
 70:   for (i=-bv->nc;i<j;i++) {
 71:     if (which && i>=0 && !which[i]) continue;
 72:     BVGetColumn(bv,i,&vi);
 73:     /* h_i = ( v, v_i ) */
 74:     if (bv->matrix) {
 75:       BV_IPMatMult(bv,w);
 76:       z = bv->Bx;
 77:     }
 78:     VecDot(z,vi,&dot);
 79:     /* v <- v - h_i v_i */
 80:     BV_SetValue(bv,i,0,c,dot);
 81:     if (bv->indef) dot /= PetscRealPart(omega[bv->nc+i]);
 82:     VecAXPY(w,-dot,vi);
 83:     BVRestoreColumn(bv,i,&vi);
 84:   }
 85:   if (nrm) { BVNormVec(bv,w,NORM_2,nrm); }
 86:   if (!v) { BVRestoreColumn(bv,j,&w); }
 87:   BV_AddCoefficients(bv,j,h,c);
 88:   if (bv->indef) {
 89:     VecRestoreArrayRead(bv->omega,&omega);
 90:   }
 91:   return(0);
 92: }

 94: /*
 95:    BVOrthogonalizeCGS1 - Compute |v'| (estimated), |v| and one step of CGS with
 96:    only one global synchronization
 97: */
 98: static PetscErrorCode BVOrthogonalizeCGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
 99: {
101:   PetscReal      sum,beta;

104:   /* h = W^* v ; alpha = (v, v) */
105:   bv->k = j;
106:   if (onorm || norm) {
107:     if (!v) {
108:       BVDotColumnInc(bv,j,c);
109:       BV_SquareRoot(bv,j,c,&beta);
110:     } else {
111:       BVDotVec(bv,v,c);
112:       BVNormVec(bv,v,NORM_2,&beta);
113:     }
114:   } else {
115:     if (!v) { BVDotColumn(bv,j,c); }
116:     else { BVDotVec(bv,v,c); }
117:   }

119:   /* q = v - V h */
120:   if (bv->indef) { BV_ApplySignature(bv,j,c,PETSC_TRUE); }
121:   if (!v) { BVMultColumn(bv,-1.0,1.0,j,c); }
122:   else { BVMultVec(bv,-1.0,1.0,v,c); }
123:   if (bv->indef) { BV_ApplySignature(bv,j,c,PETSC_FALSE); }

125:   /* compute |v| */
126:   if (onorm) *onorm = beta;

128:   if (norm) {
129:     if (bv->indef) {
130:       BV_NormVecOrColumn(bv,j,v,norm);
131:     } else {
132:       /* estimate |v'| from |v| */
133:       BV_SquareSum(bv,j,c,&sum);
134:       *norm = beta*beta-sum;
135:       if (*norm <= 0.0) {
136:         BV_NormVecOrColumn(bv,j,v,norm);
137:       } else *norm = PetscSqrtReal(*norm);
138:     }
139:   }
140:   BV_AddCoefficients(bv,j,h,c);
141:   return(0);
142: }

144: #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)

146: /*
147:    BVOrthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt

149:    j      - the index of the column to orthogonalize (cannot use both j and v)
150:    v      - the vector to orthogonalize (cannot use both j and v)
151:    which  - logical array indicating selected columns (only used in MGS)
152:    norm   - (optional) norm of the vector after being orthogonalized
153:    lindep - (optional) flag indicating possible linear dependence
154: */
155: static PetscErrorCode BVOrthogonalizeGS(BV bv,PetscInt j,Vec v,PetscBool *which,PetscReal *norm,PetscBool *lindep)
156: {
158:   PetscScalar    *h,*c,*omega;
159:   PetscReal      onrm,nrm;
160:   PetscInt       k,l;
161:   PetscBool      mgs,dolindep,signature;

164:   if (v) {
165:     k = bv->k;
166:     h = bv->h;
167:     c = bv->c;
168:   } else {
169:     k = j;
170:     h = NULL;
171:     c = NULL;
172:   }

174:   mgs = (bv->orthog_type==BV_ORTHOG_MGS)? PETSC_TRUE: PETSC_FALSE;

176:   /* if indefinite inner product, skip the computation of lindep */
177:   if (bv->indef && lindep) *lindep = PETSC_FALSE;
178:   dolindep = (!bv->indef && lindep)? PETSC_TRUE: PETSC_FALSE;

180:   /* if indefinite and we are orthogonalizing a column, the norm must always be computed */
181:   signature = (bv->indef && !v)? PETSC_TRUE: PETSC_FALSE;

183:   BV_CleanCoefficients(bv,k,h);

185:   switch (bv->orthog_ref) {

187:   case BV_ORTHOG_REFINE_IFNEEDED:
188:     BVOrthogonalizeGS1(bv,k,v,which,h,c,&onrm,&nrm);
189:     /* repeat if ||q|| < eta ||h|| */
190:     l = 1;
191:     while (l<3 && nrm && PetscAbsReal(nrm) < bv->orthog_eta*PetscAbsReal(onrm)) {
192:       l++;
193:       if (mgs||bv->indef) onrm = nrm;
194:       BVOrthogonalizeGS1(bv,k,v,which,h,c,(mgs||bv->indef)?NULL:&onrm,&nrm);
195:     }
196:     /* linear dependence check: criterion not satisfied in the last iteration */
197:     if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
198:     break;

200:   case BV_ORTHOG_REFINE_NEVER:
201:     BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL);
202:     /* compute ||v|| */
203:     if (norm || dolindep || signature) {
204:       BV_NormVecOrColumn(bv,k,v,&nrm);
205:     }
206:     /* linear dependence check: just test for exactly zero norm */
207:     if (dolindep) *lindep = PetscNot(nrm);
208:     break;

210:   case BV_ORTHOG_REFINE_ALWAYS:
211:     BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL);
212:     BVOrthogonalizeGS1(bv,k,v,which,h,c,dolindep?&onrm:NULL,(norm||dolindep||signature)?&nrm:NULL);
213:     /* linear dependence check: criterion not satisfied in the second iteration */
214:     if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
215:     break;
216:   }
217:   if (signature) {
218:     VecGetArray(bv->omega,&omega);
219:     omega[bv->nc+k] = (nrm<0.0)? -1.0: 1.0;
220:     VecRestoreArray(bv->omega,&omega);
221:   }
222:   if (norm) {
223:     *norm = nrm;
224:     if (!v) { /* store norm value next to the orthogonalization coefficients */
225:       if (dolindep && *lindep) { BV_SetValue(bv,k,k,h,0.0); }
226:       else { BV_SetValue(bv,k,k,h,nrm); }
227:     }
228:   }
229:   return(0);
230: }

232: /*@
233:    BVOrthogonalizeVec - Orthogonalize a given vector with respect to all
234:    active columns.

236:    Collective on BV

238:    Input Parameters:
239: +  bv     - the basis vectors context
240: -  v      - the vector

242:    Output Parameters:
243: +  H      - (optional) coefficients computed during orthogonalization
244: .  norm   - (optional) norm of the vector after being orthogonalized
245: -  lindep - (optional) flag indicating that refinement did not improve the quality
246:             of orthogonalization

248:    Notes:
249:    This function is equivalent to BVOrthogonalizeColumn() but orthogonalizes
250:    a vector as an argument rather than taking one of the BV columns. The
251:    vector is orthogonalized against all active columns (k) and the constraints.
252:    If H is given, it must have enough space to store k-l coefficients, where l
253:    is the number of leading columns.

255:    In the case of an indefinite inner product, the lindep parameter is not
256:    computed (set to false).

258:    Level: advanced

260: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVSetActiveColumns(), BVGetNumConstraints()
261: @*/
262: PetscErrorCode BVOrthogonalizeVec(BV bv,Vec v,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
263: {
265:   PetscInt       ksave,lsave;

271:   BVCheckSizes(bv,1);

275:   PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
276:   ksave = bv->k;
277:   lsave = bv->l;
278:   bv->l = -bv->nc;  /* must also orthogonalize against constraints and leading columns */
279:   BV_AllocateCoeffs(bv);
280:   BV_AllocateSignature(bv);
281:   BVOrthogonalizeGS(bv,0,v,NULL,norm,lindep);
282:   bv->k = ksave;
283:   bv->l = lsave;
284:   if (H) { BV_StoreCoefficients(bv,bv->k,bv->h,H); }
285:   PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
286:   return(0);
287: }

289: /*@
290:    BVOrthogonalizeColumn - Orthogonalize one of the column vectors with respect to
291:    the previous ones.

293:    Collective on BV

295:    Input Parameters:
296: +  bv     - the basis vectors context
297: -  j      - index of column to be orthogonalized

299:    Output Parameters:
300: +  H      - (optional) coefficients computed during orthogonalization
301: .  norm   - (optional) norm of the vector after being orthogonalized
302: -  lindep - (optional) flag indicating that refinement did not improve the quality
303:             of orthogonalization

305:    Notes:
306:    This function applies an orthogonal projector to project vector V[j] onto
307:    the orthogonal complement of the span of the columns of V[0..j-1],
308:    where V[.] are the vectors of BV. The columns V[0..j-1] are assumed to be
309:    mutually orthonormal.

311:    Leading columns V[0..l-1] also participate in the orthogonalization, as well
312:    as the constraints. If H is given, it must have enough space to store
313:    j-l+1 coefficients (the last coefficient will contain the value norm, unless
314:    the norm argument is NULL).

316:    If a non-standard inner product has been specified with BVSetMatrix(),
317:    then the vector is B-orthogonalized, using the non-standard inner product
318:    defined by matrix B. The output vector satisfies V[j]'*B*V[0..j-1] = 0.

320:    This routine does not normalize the resulting vector, see BVOrthonormalizeColumn().

322:    In the case of an indefinite inner product, the lindep parameter is not
323:    computed (set to false).

325:    Level: advanced

327: .seealso: BVSetOrthogonalization(), BVSetMatrix(), BVSetActiveColumns(), BVOrthogonalize(), BVOrthogonalizeVec(), BVGetNumConstraints(), BVOrthonormalizeColumn()
328: @*/
329: PetscErrorCode BVOrthogonalizeColumn(BV bv,PetscInt j,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
330: {
332:   PetscInt       ksave,lsave;

338:   BVCheckSizes(bv,1);
339:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
340:   if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,bv->m);

342:   PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
343:   ksave = bv->k;
344:   lsave = bv->l;
345:   bv->l = -bv->nc;  /* must also orthogonalize against constraints and leading columns */
346:   if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
347:   BV_AllocateSignature(bv);
348:   BVOrthogonalizeGS(bv,j,NULL,NULL,norm,lindep);
349:   bv->k = ksave;
350:   bv->l = lsave;
351:   if (H) { BV_StoreCoefficients(bv,j,NULL,H); }
352:   PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
353:   PetscObjectStateIncrease((PetscObject)bv);
354:   return(0);
355: }

357: /*@
358:    BVOrthonormalizeColumn - Orthonormalize one of the column vectors with respect to
359:    the previous ones. This is equivalent to a call to BVOrthogonalizeColumn()
360:    followed by a call to BVScaleColumn() with the reciprocal of the norm.

362:    Collective on BV

364:    Input Parameters:
365: +  bv      - the basis vectors context
366: .  j       - index of column to be orthonormalized
367: -  replace - whether it is allowed to set the vector randomly

369:    Output Parameters:
370: +  norm    - (optional) norm of the vector after orthogonalization and before normalization
371: -  lindep  - (optional) flag indicating that linear dependence was determined during
372:              orthogonalization

374:    Notes:
375:    This function first orthogonalizes vector V[j] with respect to V[0..j-1],
376:    where V[.] are the vectors of BV. A byproduct of this computation is norm,
377:    the norm of the vector after orthogonalization. Secondly, it scales the
378:    vector with 1/norm, so that the resulting vector has unit norm.

380:    If after orthogonalization the vector V[j] is exactly zero, it cannot be normalized
381:    because norm=0. In that case, it could be left as zero or replaced by a random
382:    vector that is then orthonormalized. The latter is achieved by setting the
383:    argument replace to TRUE. The vector will be replaced by a random vector also
384:    if lindep was set to TRUE, even if the norm is not exaclty zero.

386:    If the vector has been replaced by a random vector, the output arguments norm and
387:    lindep will be set according to the orthogonalization of this new vector.

389:    Level: advanced

391: .seealso: BVOrthogonalizeColumn(), BVScaleColumn()
392: @*/
393: PetscErrorCode BVOrthonormalizeColumn(BV bv,PetscInt j,PetscBool replace,PetscReal *norm,PetscBool *lindep)
394: {
396:   PetscScalar    alpha;
397:   PetscReal      nrm;
398:   PetscInt       ksave,lsave;
399:   PetscBool      lndep;

405:   BVCheckSizes(bv,1);
406:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
407:   if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,bv->m);

409:   /* orthogonalize */
410:   PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
411:   ksave = bv->k;
412:   lsave = bv->l;
413:   bv->l = -bv->nc;  /* must also orthogonalize against constraints and leading columns */
414:   if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
415:   BV_AllocateSignature(bv);
416:   BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep);
417:   if (replace && (nrm==0.0 || lndep)) {
418:     PetscInfo(bv,"Vector was linearly dependent, generating a new random vector\n");
419:     BVSetRandomColumn(bv,j);
420:     BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep);
421:     if (nrm==0.0 || lndep) {  /* yet another attempt */
422:       BVSetRandomColumn(bv,j);
423:       BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep);
424:     }
425:   }
426:   bv->k = ksave;
427:   bv->l = lsave;
428:   PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);

430:   /* scale */
431:   if (nrm!=1.0 && nrm!=0.0) {
432:     alpha = 1.0/nrm;
433:     PetscLogEventBegin(BV_Scale,bv,0,0,0);
434:     if (bv->n) {
435:       (*bv->ops->scale)(bv,j,alpha);
436:     }
437:     PetscLogEventEnd(BV_Scale,bv,0,0,0);
438:   }
439:   if (norm) *norm = nrm;
440:   if (lindep) *lindep = lndep;
441:   PetscObjectStateIncrease((PetscObject)bv);
442:   return(0);
443: }

445: /*@
446:    BVOrthogonalizeSomeColumn - Orthogonalize one of the column vectors with
447:    respect to some of the previous ones.

449:    Collective on BV

451:    Input Parameters:
452: +  bv     - the basis vectors context
453: .  j      - index of column to be orthogonalized
454: -  which  - logical array indicating selected columns

456:    Output Parameters:
457: +  H      - (optional) coefficients computed during orthogonalization
458: .  norm   - (optional) norm of the vector after being orthogonalized
459: -  lindep - (optional) flag indicating that refinement did not improve the quality
460:             of orthogonalization

462:    Notes:
463:    This function is similar to BVOrthogonalizeColumn(), but V[j] is
464:    orthogonalized only against columns V[i] having which[i]=PETSC_TRUE.
465:    The length of array which must be j at least.

467:    The use of this operation is restricted to MGS orthogonalization type.

469:    In the case of an indefinite inner product, the lindep parameter is not
470:    computed (set to false).

472:    Level: advanced

474: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization()
475: @*/
476: PetscErrorCode BVOrthogonalizeSomeColumn(BV bv,PetscInt j,PetscBool *which,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
477: {
479:   PetscInt       ksave,lsave;

486:   BVCheckSizes(bv,1);
487:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
488:   if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,bv->m);
489:   if (bv->orthog_type!=BV_ORTHOG_MGS) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation only available for MGS orthogonalization");

491:   PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
492:   ksave = bv->k;
493:   lsave = bv->l;
494:   bv->l = -bv->nc;  /* must also orthogonalize against constraints and leading columns */
495:   if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
496:   BV_AllocateSignature(bv);
497:   BVOrthogonalizeGS(bv,j,NULL,which,norm,lindep);
498:   bv->k = ksave;
499:   bv->l = lsave;
500:   if (H) { BV_StoreCoefficients(bv,j,NULL,H); }
501:   PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
502:   PetscObjectStateIncrease((PetscObject)bv);
503:   return(0);
504: }

506: /*
507:    Block Gram-Schmidt: V2 = V2 - V1*R12, where R12 = V1'*V2
508:  */
509: static PetscErrorCode BVOrthogonalize_BlockGS(BV V,Mat R)
510: {
512:   BV             V1;

515:   BVGetSplit(V,&V1,NULL);
516:   BVDot(V,V1,R);
517:   BVMult(V,-1.0,1.0,V1,R);
518:   BVRestoreSplit(V,&V1,NULL);
519:   return(0);
520: }

522: /*
523:    Orthogonalize a set of vectors with Gram-Schmidt, column by column.
524:  */
525: static PetscErrorCode BVOrthogonalize_GS(BV V,Mat R)
526: {
528:   PetscScalar    *r=NULL;
529:   PetscReal      norm;
530:   PetscInt       j,ldr,lsave;
531:   Vec            v,w;

534:   if (R) {
535:     MatGetSize(R,&ldr,NULL);
536:     MatDenseGetArray(R,&r);
537:   }
538:   if (V->matrix) {
539:     BVGetCachedBV(V,&V->cached);
540:     BVSetActiveColumns(V->cached,V->l,V->k);
541:   }
542:   for (j=V->l;j<V->k;j++) {
543:     if (V->matrix && V->orthog_type==BV_ORTHOG_MGS) {  /* fill cached BV */
544:       BVGetColumn(V->cached,j,&v);
545:       BVGetColumn(V,j,&w);
546:       MatMult(V->matrix,w,v);
547:       BVRestoreColumn(V,j,&w);
548:       BVRestoreColumn(V->cached,j,&v);
549:     }
550:     if (R) {
551:       BVOrthogonalizeColumn(V,j,NULL,&norm,NULL);
552:       lsave = V->l;
553:       V->l = -V->nc;
554:       BV_StoreCoefficients(V,j,NULL,r+j*ldr);
555:       V->l = lsave;
556:       r[j+j*ldr] = norm;
557:     } else {
558:       BVOrthogonalizeColumn(V,j,NULL,&norm,NULL);
559:     }
560:     if (!norm) SETERRQ(PetscObjectComm((PetscObject)V),1,"Breakdown in BVOrthogonalize due to a linearly dependent column");
561:     if (V->matrix && V->orthog_type==BV_ORTHOG_CGS) {  /* fill cached BV */
562:       BVGetColumn(V->cached,j,&v);
563:       VecCopy(V->Bx,v);
564:       BVRestoreColumn(V->cached,j,&v);
565:     }
566:     BVScaleColumn(V,j,1.0/norm);
567:   }
568:   if (R) { MatDenseRestoreArray(R,&r); }
569:   return(0);
570: }

572: /*
573:   BV_GetBufferMat - Create auxiliary seqdense matrix that wraps the bv->buffer.
574: */
575: PETSC_STATIC_INLINE PetscErrorCode BV_GetBufferMat(BV bv)
576: {
578:   PetscInt       ld;
579:   PetscScalar    *array;

582:   if (!bv->Abuffer) {
583:     if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
584:     ld = bv->m+bv->nc;
585:     VecGetArray(bv->buffer,&array);
586:     MatCreateSeqDense(PETSC_COMM_SELF,ld,bv->m,array,&bv->Abuffer);
587:     VecRestoreArray(bv->buffer,&array);
588:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Abuffer);
589:   }
590:   return(0);
591: }

593: /*
594:    BV_StoreCoeffsBlock_Default - Copy the contents of the BV buffer to a dense Mat
595:    provided by the caller. Only columns l:k-1 are copied, restricting to the upper
596:    triangular part if tri=PETSC_TRUE.
597: */
598: PETSC_STATIC_INLINE PetscErrorCode BV_StoreCoeffsBlock_Default(BV bv,Mat R,PetscBool tri)
599: {
600:   PetscErrorCode    ierr;
601:   const PetscScalar *bb;
602:   PetscScalar       *rr;
603:   PetscInt          j,ldr,ldb;

606:   MatGetSize(R,&ldr,NULL);
607:   MatDenseGetArray(R,&rr);
608:   ldb  = bv->m+bv->nc;
609:   VecGetArrayRead(bv->buffer,&bb);
610:   for (j=bv->l;j<bv->k;j++) {
611:     PetscMemcpy(rr+j*ldr,bb+j*ldb,((tri?(j+1):bv->k)+bv->nc)*sizeof(PetscScalar));
612:   }
613:   VecRestoreArrayRead(bv->buffer,&bb);
614:   MatDenseRestoreArray(R,&rr);
615:   return(0);
616: }

618: /*
619:    Orthogonalize a set of vectors with Cholesky: R=chol(V'*V), Q=V*inv(R)
620:  */
621: static PetscErrorCode BVOrthogonalize_Chol(BV V,Mat Rin)
622: {
624:   Mat            R,S;

627:   BV_GetBufferMat(V);
628:   R = V->Abuffer;
629:   if (Rin) S = Rin;   /* use Rin as a workspace for S */
630:   else S = R;
631:   if (V->l) { BVOrthogonalize_BlockGS(V,R); }
632:   BVDot(V,V,R);
633:   BVMatCholInv_LAPACK_Private(V,R,S);
634:   BVMultInPlace(V,S,V->l,V->k);
635:   if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE); }
636:   return(0);
637: }

639: /*
640:    Orthogonalize a set of vectors with the Tall-Skinny QR method
641:  */
642: static PetscErrorCode BVOrthogonalize_TSQR(BV V,Mat Rin)
643: {
645:   PetscScalar    *pv,*r=NULL;
646:   PetscInt       ldr;
647:   Mat            R;

650:   BV_GetBufferMat(V);
651:   R = V->Abuffer;
652:   if (V->l) { BVOrthogonalize_BlockGS(V,R); }
653:   MatGetSize(R,&ldr,NULL);
654:   MatDenseGetArray(R,&r);
655:   BVGetArray(V,&pv);
656:   BVOrthogonalize_LAPACK_TSQR(V,V->n,V->k-V->l,pv+(V->nc+V->l)*V->n,r+V->l*ldr+V->l,ldr);
657:   BVRestoreArray(V,&pv);
658:   MatDenseRestoreArray(R,&r);
659:   if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE); }
660:   return(0);
661: }

663: /*
664:    Orthogonalize a set of vectors with TSQR, but computing R only, then doing Q=V*inv(R)
665:  */
666: static PetscErrorCode BVOrthogonalize_TSQRCHOL(BV V,Mat Rin)
667: {
669:   PetscScalar    *pv,*r=NULL;
670:   PetscInt       ldr;
671:   Mat            R,S;

674:   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) { BVOrthogonalize_BlockGS(V,R); }
679:   MatGetSize(R,&ldr,NULL);
680:   MatDenseGetArray(R,&r);
681:   BVGetArray(V,&pv);
682:   BVOrthogonalize_LAPACK_TSQR_OnlyR(V,V->n,V->k-V->l,pv+(V->nc+V->l)*V->n,r+V->l*ldr+V->l,ldr);
683:   BVRestoreArray(V,&pv);
684:   MatDenseRestoreArray(R,&r);
685:   BVMatTriInv_LAPACK_Private(V,R,S);
686:   BVMultInPlace(V,S,V->l,V->k);
687:   if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE); }
688:   return(0);
689: }

691: /*
692:    Orthogonalize a set of vectors with SVQB
693:  */
694: static PetscErrorCode BVOrthogonalize_SVQB(BV V,Mat Rin)
695: {
697:   Mat            R,S;

700:   BV_GetBufferMat(V);
701:   R = V->Abuffer;
702:   if (Rin) S = Rin;   /* use Rin as a workspace for S */
703:   else S = R;
704:   if (V->l) { BVOrthogonalize_BlockGS(V,R); }
705:   BVDot(V,V,R);
706:   BVMatSVQB_LAPACK_Private(V,R,S);
707:   BVMultInPlace(V,S,V->l,V->k);
708:   if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_FALSE); }
709:   return(0);
710: }

712: /*@
713:    BVOrthogonalize - Orthogonalize all columns (starting from the leading ones),
714:    that is, compute the QR decomposition.

716:    Collective on BV

718:    Input Parameters:
719: +  V - basis vectors to be orthogonalized (or B-orthogonalized)
720: -  R - a sequential dense matrix (or NULL)

722:    Output Parameters:
723: +  V - the modified basis vectors
724: -  R - (if not NULL) the triangular factor of the QR decomposition

726:    Notes:
727:    On input, matrix R must be a square sequential dense Mat, with at least as many
728:    rows and columns as the number of active columns of V. The output satisfies
729:    V0 = V*R (where V0 represent the input V) and V'*V = I (or V'*B*V = I if an
730:    inner product matrix B has been specified with BVSetMatrix()).

732:    If V has leading columns, then they are not modified (are assumed to be already
733:    orthonormal) and the leading columns of R are not referenced. Let the
734:    decomposition be
735: .vb
736:    [ V01 V02 ] = [ V1 V2 ] [ R11 R12 ]
737:                            [  0  R22 ]
738: .ve
739:    then V1 is left unchanged (equal to V01) as well as R11 (it should satisfy
740:    V01 = V1*R11).

742:    Can pass NULL if R is not required.

744:    The method to be used for block orthogonalization can be set with
745:    BVSetOrthogonalization(). If set to GS, the computation is done column by
746:    column with successive calls to BVOrthogonalizeColumn(). Note that in the
747:    SVQB method the R factor is not upper triangular.

749:    If V is rank-deficient or very ill-conditioned, that is, one or more columns are
750:    (almost) linearly dependent with respect to the rest, then the algorithm may
751:    break down or result in larger numerical error. Linearly dependent columns are
752:    essentially replaced by random directions, and the corresponding diagonal entry
753:    in R is set to (nearly) zero.

755:    Level: intermediate

757: .seealso: BVOrthogonalizeColumn(), BVOrthogonalizeVec(), BVSetMatrix(), BVSetActiveColumns(), BVSetOrthogonalization(), BVOrthogBlockType
758: @*/
759: PetscErrorCode BVOrthogonalize(BV V,Mat R)
760: {
762:   PetscBool      match;
763:   PetscInt       m,n;

768:   BVCheckSizes(V,1);
769:   if (R) {
772:     PetscObjectTypeCompare((PetscObject)R,MATSEQDENSE,&match);
773:     if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Mat argument must be of type seqdense");
774:     MatGetSize(R,&m,&n);
775:     if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument is not square, it has %D rows and %D columns",m,n);
776:     if (n<V->k) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat size %D is smaller than the number of BV active columns %D",n,V->k);
777:   }
778:   if (V->nc) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Not implemented for BV with constraints, use BVOrthogonalizeColumn() instead");

780:   PetscLogEventBegin(BV_Orthogonalize,V,R,0,0);
781:   switch (V->orthog_block) {
782:   case BV_ORTHOG_BLOCK_GS: /* proceed column by column with Gram-Schmidt */
783:     BVOrthogonalize_GS(V,R);
784:     break;
785:   case BV_ORTHOG_BLOCK_CHOL:
786:     BVOrthogonalize_Chol(V,R);
787:     break;
788:   case BV_ORTHOG_BLOCK_TSQR:
789:     if (V->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
790:     BVOrthogonalize_TSQR(V,R);
791:     break;
792:   case BV_ORTHOG_BLOCK_TSQRCHOL:
793:     if (V->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
794:     BVOrthogonalize_TSQRCHOL(V,R);
795:     break;
796:   case BV_ORTHOG_BLOCK_SVQB:
797:     BVOrthogonalize_SVQB(V,R);
798:     break;
799:   }
800:   PetscLogEventEnd(BV_Orthogonalize,V,R,0,0);
801:   PetscObjectStateIncrease((PetscObject)V);
802:   return(0);
803: }