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: }