Actual source code: bvbasic.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:    Basic BV routines
     12: */

     14: #include <slepc/private/bvimpl.h>

     16: PetscBool         BVRegisterAllCalled = PETSC_FALSE;
     17: PetscFunctionList BVList = NULL;

     19: /*@
     20:    BVSetType - Selects the type for the `BV` object.

     22:    Logically Collective

     24:    Input Parameters:
     25: +  bv   - the basis vectors context
     26: -  type - a known type

     28:    Options Database Key:
     29: .  -bv_type type - sets `BV` type

     31:    Level: intermediate

     33: .seealso: [](sec:bv), `BVGetType()`
     34: @*/
     35: PetscErrorCode BVSetType(BV bv,BVType type)
     36: {
     37:   PetscErrorCode (*r)(BV);
     38:   PetscBool      match;

     40:   PetscFunctionBegin;
     42:   PetscAssertPointer(type,2);

     44:   PetscCall(PetscObjectTypeCompare((PetscObject)bv,type,&match));
     45:   if (match) PetscFunctionReturn(PETSC_SUCCESS);
     46:   PetscCall(PetscStrcmp(type,BVTENSOR,&match));
     47:   PetscCheck(!match,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Use BVCreateTensor() to create a BV of type tensor");

     49:   PetscCall(PetscFunctionListFind(BVList,type,&r));
     50:   PetscCheck(r,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested BV type %s",type);

     52:   PetscTryTypeMethod(bv,destroy);
     53:   PetscCall(PetscMemzero(bv->ops,sizeof(struct _BVOps)));

     55:   PetscCall(PetscObjectChangeTypeName((PetscObject)bv,type));
     56:   if (bv->n < 0 && bv->N < 0) {
     57:     bv->ops->create = r;
     58:   } else {
     59:     PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
     60:     PetscCall((*r)(bv));
     61:     PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
     62:   }
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: /*@
     67:    BVGetType - Gets the `BV` type name (as a string) from the `BV` context.

     69:    Not Collective

     71:    Input Parameter:
     72: .  bv - the basis vectors context

     74:    Output Parameter:
     75: .  type - name of the type of basis vectors

     77:    Level: intermediate

     79: .seealso: [](sec:bv), `BVSetType()`
     80: @*/
     81: PetscErrorCode BVGetType(BV bv,BVType *type)
     82: {
     83:   PetscFunctionBegin;
     85:   PetscAssertPointer(type,2);
     86:   *type = ((PetscObject)bv)->type_name;
     87:   PetscFunctionReturn(PETSC_SUCCESS);
     88: }

     90: /*@
     91:    BVSetSizes - Sets the local and global sizes, and the number of columns.

     93:    Collective

     95:    Input Parameters:
     96: +  bv - the basis vectors
     97: .  n  - the local size (or `PETSC_DECIDE` to have it set)
     98: .  N  - the global size (or `PETSC_DECIDE`)
     99: -  m  - the number of columns

    101:    Notes:
    102:    `n` and `N` cannot be both `PETSC_DECIDE`.
    103:    If one process calls this with `N` equal to `PETSC_DECIDE` then all processes must,
    104:    otherwise the program will hang.

    106:    Level: beginner

    108: .seealso: [](sec:bv), `BVSetSizesFromVec()`, `BVGetSizes()`, `BVResize()`
    109: @*/
    110: PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)
    111: {
    112:   PetscInt       ma;
    113:   PetscMPIInt    size;

    115:   PetscFunctionBegin;
    119:   PetscCheck(N<0 || n<=N,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT,n,N);
    120:   PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
    121:   PetscCheck((bv->n<0 && bv->N<0) || (bv->n==n && bv->N==N),PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change/reset vector sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global",n,N,bv->n,bv->N);
    122:   PetscCheck(bv->m<=0 || bv->m==m,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change the number of columns to %" PetscInt_FMT " after previously setting it to %" PetscInt_FMT "; use BVResize()",m,bv->m);
    123:   PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
    124:   bv->n = n;
    125:   bv->N = N;
    126:   bv->m = m;
    127:   bv->k = m;
    128:   /* create layout and get actual dimensions */
    129:   PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)bv),&bv->map));
    130:   PetscCall(PetscLayoutSetSize(bv->map,bv->N));
    131:   PetscCall(PetscLayoutSetLocalSize(bv->map,bv->n));
    132:   PetscCall(PetscLayoutSetUp(bv->map));
    133:   PetscCall(PetscLayoutGetSize(bv->map,&bv->N));
    134:   PetscCall(PetscLayoutGetLocalSize(bv->map,&bv->n));
    135:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
    136:   PetscCall(BVSetVecType(bv,(size==1)?VECSEQ:VECMPI));
    137:   if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
    138:     PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
    139:     PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
    140:   }
    141:   if (bv->ops->create) {
    142:     PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
    143:     PetscUseTypeMethod(bv,create);
    144:     PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
    145:     bv->ops->create = NULL;
    146:     bv->defersfo = PETSC_FALSE;
    147:   }
    148:   PetscFunctionReturn(PETSC_SUCCESS);
    149: }

    151: /*@
    152:    BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
    153:    Local and global sizes are specified indirectly by passing a template vector.

    155:    Collective

    157:    Input Parameters:
    158: +  bv - the basis vectors
    159: .  t  - the template vector
    160: -  m  - the number of columns

    162:    Notes:
    163:    Apart from the local and global sizes, the template vector `t` is also used
    164:    to get the vector type, which will be set as the `BV` vector type with
    165:    `BVSetVecType()`. This is relevant in cases where computation must be done
    166:    on the GPU.

    168:    If `t` is a nested vector (type `VECNEST`), then the vector type will be
    169:    obtained from the first subvector. The rationale is that `BV` will operate
    170:    internally with regular vectors, even though the user problem is defined via
    171:    nested matrices and vectors.

    173:    Level: beginner

    175: .seealso: [](sec:bv), `BVSetSizes()`, `BVGetSizes()`, `BVResize()`, `BVSetVecType()`
    176: @*/
    177: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
    178: {
    179:   PetscInt       ma;
    180:   PetscLayout    map;
    181:   VecType        vtype;
    182:   PetscBool      nest;
    183:   Vec            v;

    185:   PetscFunctionBegin;
    188:   PetscCheckSameComm(bv,1,t,2);
    190:   PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
    191:   PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
    192:   PetscCall(PetscObjectTypeCompare((PetscObject)t,VECNEST,&nest));
    193:   if (nest) {
    194:     PetscCall(VecNestGetSubVec(t,0,&v));
    195:     PetscCall(VecGetType(v,&vtype));
    196:   } else PetscCall(VecGetType(t,&vtype));
    197:   PetscCall(BVSetVecType(bv,vtype));
    198:   PetscCall(VecGetLayout(t,&map));
    199:   PetscCall(PetscLayoutReference(map,&bv->map));
    200:   PetscCall(VecGetSize(t,&bv->N));
    201:   PetscCall(VecGetLocalSize(t,&bv->n));
    202:   if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
    203:     PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
    204:     PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
    205:   }
    206:   bv->m = m;
    207:   bv->k = m;
    208:   if (bv->ops->create) {
    209:     PetscUseTypeMethod(bv,create);
    210:     bv->ops->create = NULL;
    211:     bv->defersfo = PETSC_FALSE;
    212:   }
    213:   PetscFunctionReturn(PETSC_SUCCESS);
    214: }

    216: /*@
    217:    BVGetSizes - Returns the local and global sizes, and the number of columns.

    219:    Not Collective

    221:    Input Parameter:
    222: .  bv - the basis vectors

    224:    Output Parameters:
    225: +  n  - the local size
    226: .  N  - the global size
    227: -  m  - the number of columns

    229:    Note:
    230:    Normal usage requires that `bv` has already been given its sizes, otherwise
    231:    the call fails. However, this function can also be used to determine if
    232:    a `BV` object has been initialized completely (sizes and type). For this,
    233:    call with `n=NULL` and `N=NULL`, then a return value of `m=0` indicates that
    234:    the `BV` object is not ready for use yet.

    236:    Level: beginner

    238: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
    239: @*/
    240: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
    241: {
    242:   PetscFunctionBegin;
    243:   if (!bv) {
    244:     if (m && !n && !N) *m = 0;
    245:     PetscFunctionReturn(PETSC_SUCCESS);
    246:   }
    248:   if (n || N) BVCheckSizes(bv,1);
    249:   if (n) *n = bv->n;
    250:   if (N) *N = bv->N;
    251:   if (m) *m = bv->m;
    252:   if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
    253:   PetscFunctionReturn(PETSC_SUCCESS);
    254: }

    256: /*@
    257:    BVSetNumConstraints - Set the number of constraints.

    259:    Logically Collective

    261:    Input Parameters:
    262: +  V  - basis vectors
    263: -  nc - number of constraints

    265:    Notes:
    266:    This function sets the number of constraints to `nc` and marks all remaining
    267:    columns as regular. Normal usage would be to call `BVInsertConstraints()` instead.

    269:    If `nc` is smaller than the previously set value, then some of the constraints
    270:    are discarded. In particular, using `nc=0` removes all constraints preserving
    271:    the content of regular columns.

    273:    Level: developer

    275: .seealso: [](sec:bv), `BVInsertConstraints()`
    276: @*/
    277: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
    278: {
    279:   PetscInt       total,diff,i;
    280:   Vec            x,y;

    282:   PetscFunctionBegin;
    285:   PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
    287:   BVCheckSizes(V,1);
    288:   PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");

    290:   diff = nc-V->nc;
    291:   if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
    292:   total = V->nc+V->m;
    293:   PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
    294:   if (diff<0) {  /* lessen constraints, shift contents of BV */
    295:     for (i=0;i<V->m;i++) {
    296:       PetscCall(BVGetColumn(V,i,&x));
    297:       PetscCall(BVGetColumn(V,i+diff,&y));
    298:       PetscCall(VecCopy(x,y));
    299:       PetscCall(BVRestoreColumn(V,i,&x));
    300:       PetscCall(BVRestoreColumn(V,i+diff,&y));
    301:     }
    302:   }
    303:   V->nc = nc;
    304:   V->ci[0] = -V->nc-1;
    305:   V->ci[1] = -V->nc-1;
    306:   V->m = total-nc;
    307:   V->l = PetscMin(V->l,V->m);
    308:   V->k = PetscMin(V->k,V->m);
    309:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: /*@
    314:    BVGetNumConstraints - Returns the number of constraints.

    316:    Not Collective

    318:    Input Parameter:
    319: .  bv - the basis vectors

    321:    Output Parameter:
    322: .  nc - the number of constraints

    324:    Level: advanced

    326: .seealso: [](sec:bv), `BVGetSizes()`, `BVInsertConstraints()`
    327: @*/
    328: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
    329: {
    330:   PetscFunctionBegin;
    332:   PetscAssertPointer(nc,2);
    333:   *nc = bv->nc;
    334:   PetscFunctionReturn(PETSC_SUCCESS);
    335: }

    337: /*@
    338:    BVResize - Change the number of columns.

    340:    Collective

    342:    Input Parameters:
    343: +  bv   - the basis vectors
    344: .  m    - the new number of columns
    345: -  copy - a flag indicating whether current values should be kept

    347:    Note:
    348:    Internal storage is reallocated. If the `copy` flag is set to true, then
    349:    the contents are copied to the leading part of the new space.

    351:    Level: advanced

    353: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
    354: @*/
    355: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
    356: {
    357:   PetscScalar       *array;
    358:   const PetscScalar *omega;
    359:   Vec               v;

    361:   PetscFunctionBegin;
    366:   PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
    367:   PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
    368:   if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
    369:   BVCheckOp(bv,1,resize);

    371:   PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
    372:   PetscUseTypeMethod(bv,resize,m,copy);
    373:   PetscCall(VecDestroy(&bv->buffer));
    374:   PetscCall(BVDestroy(&bv->cached));
    375:   PetscCall(PetscFree2(bv->h,bv->c));
    376:   if (bv->omega) {
    377:     if (bv->cuda) {
    378: #if defined(PETSC_HAVE_CUDA)
    379:       PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
    380: #else
    381:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
    382: #endif
    383:     } else if (bv->hip) {
    384: #if defined(PETSC_HAVE_HIP)
    385:       PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
    386: #else
    387:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
    388: #endif
    389:     } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
    390:     if (copy) {
    391:       PetscCall(VecGetArray(v,&array));
    392:       PetscCall(VecGetArrayRead(bv->omega,&omega));
    393:       PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
    394:       PetscCall(VecRestoreArrayRead(bv->omega,&omega));
    395:       PetscCall(VecRestoreArray(v,&array));
    396:     } else PetscCall(VecSet(v,1.0));
    397:     PetscCall(VecDestroy(&bv->omega));
    398:     bv->omega = v;
    399:   }
    400:   bv->m = m;
    401:   bv->k = PetscMin(bv->k,m);
    402:   bv->l = PetscMin(bv->l,m);
    403:   PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
    404:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
    405:   PetscFunctionReturn(PETSC_SUCCESS);
    406: }

    408: /*@
    409:    BVSetActiveColumns - Specify the columns that will be involved in operations.

    411:    Logically Collective

    413:    Input Parameters:
    414: +  bv - the basis vectors context
    415: .  l  - number of leading columns
    416: -  k  - number of active columns

    418:    Notes:
    419:    In operations such as `BVMult()` or `BVDot()`, only the first `k` columns are
    420:    considered. This is useful when the `BV` is filled from left to right, so
    421:    the last `m-k` columns do not have relevant information.

    423:    Also in operations such as `BVMult()` or `BVDot()`, the first `l` columns are
    424:    normally not included in the computation. See the manual page of each
    425:    operation.

    427:    In orthogonalization operations, the first `l` columns are treated
    428:    differently, they participate in the orthogonalization but the computed
    429:    coefficients are not stored.

    431:    Use `PETSC_CURRENT` to leave any of the values unchanged. Use `PETSC_DETERMINE`
    432:    to set `l` to the minimum value (`0`) and `k` to the maximum (`m`).

    434:    Level: intermediate

    436: .seealso: [](sec:bv), `BVGetActiveColumns()`, `BVSetSizes()`
    437: @*/
    438: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
    439: {
    440:   PetscFunctionBegin;
    444:   BVCheckSizes(bv,1);
    445:   if (PetscUnlikely(k == PETSC_DETERMINE)) {
    446:     bv->k = bv->m;
    447:   } else if (k != PETSC_CURRENT) {
    448:     PetscCheck(k>=0 && k<=bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k (%" PetscInt_FMT "). Must be between 0 and m (%" PetscInt_FMT ")",k,bv->m);
    449:     bv->k = k;
    450:   }
    451:   if (PetscUnlikely(l == PETSC_DETERMINE)) {
    452:     bv->l = 0;
    453:   } else if (l != PETSC_CURRENT) {
    454:     PetscCheck(l>=0 && l<=bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l (%" PetscInt_FMT "). Must be between 0 and k (%" PetscInt_FMT ")",l,bv->k);
    455:     bv->l = l;
    456:   }
    457:   PetscFunctionReturn(PETSC_SUCCESS);
    458: }

    460: /*@
    461:    BVGetActiveColumns - Returns the current active dimensions.

    463:    Not Collective

    465:    Input Parameter:
    466: .  bv - the basis vectors context

    468:    Output Parameters:
    469: +  l  - number of leading columns
    470: -  k  - number of active columns

    472:    Level: intermediate

    474: .seealso: [](sec:bv), `BVSetActiveColumns()`
    475: @*/
    476: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
    477: {
    478:   PetscFunctionBegin;
    480:   if (l) *l = bv->l;
    481:   if (k) *k = bv->k;
    482:   PetscFunctionReturn(PETSC_SUCCESS);
    483: }

    485: /*@
    486:    BVSetMatrix - Specifies the inner product to be used in orthogonalization.

    488:    Collective

    490:    Input Parameters:
    491: +  bv    - the basis vectors context
    492: .  B     - a Hermitian matrix (may be `NULL`)
    493: -  indef - a flag indicating if the matrix is indefinite

    495:    Notes:
    496:    This is used to specify a non-standard inner product, whose matrix
    497:    representation is given by `B`. Then, all inner products required during
    498:    orthogonalization are computed as $(x,y)_B=y^*Bx$ rather than the
    499:    standard form $(x,y)=y^*x$.

    501:    Matrix `B` must be real symmetric (or complex Hermitian). A genuine inner
    502:    product requires that `B` is also positive (semi-)definite. However, we
    503:    also allow for an indefinite `B` (setting `indef=PETSC_TRUE`), in which
    504:    case the orthogonalization uses an indefinite inner product.

    506:    This affects operations `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, and variants.

    508:    Setting `B=NULL` has the same effect as if the identity matrix was passed.

    510:    Level: advanced

    512: .seealso: [](sec:bv), `BVGetMatrix()`, `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, `BVSetDefiniteTolerance()`
    513: @*/
    514: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
    515: {
    516:   PetscInt       m,n;

    518:   PetscFunctionBegin;
    521:   if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
    522:     if (B) {
    524:       PetscCall(MatGetLocalSize(B,&m,&n));
    525:       PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
    526:       PetscCheck(!bv->m || bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension BV %" PetscInt_FMT ", Mat %" PetscInt_FMT,bv->n,n);
    527:     }
    528:     PetscCall(PetscObjectReference((PetscObject)B));
    529:     PetscCall(MatDestroy(&bv->matrix));
    530:     bv->matrix = B;
    531:     bv->indef  = indef;
    532:     PetscCall(PetscObjectStateIncrease((PetscObject)bv));
    533:     if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
    534:     if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
    535:   }
    536:   PetscFunctionReturn(PETSC_SUCCESS);
    537: }

    539: /*@
    540:    BVGetMatrix - Retrieves the matrix representation of the inner product.

    542:    Not Collective

    544:    Input Parameter:
    545: .  bv    - the basis vectors context

    547:    Output Parameters:
    548: +  B     - the matrix of the inner product (may be `NULL`)
    549: -  indef - the flag indicating if the matrix is indefinite

    551:    Level: advanced

    553: .seealso: [](sec:bv), `BVSetMatrix()`
    554: @*/
    555: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
    556: {
    557:   PetscFunctionBegin;
    559:   if (B)     *B     = bv->matrix;
    560:   if (indef) *indef = bv->indef;
    561:   PetscFunctionReturn(PETSC_SUCCESS);
    562: }

    564: /*@
    565:    BVApplyMatrix - Multiplies a vector by the matrix representation of the
    566:    inner product.

    568:    Neighbor-wise Collective

    570:    Input Parameters:
    571: +  bv - the basis vectors context
    572: -  x  - the vector

    574:    Output Parameter:
    575: .  y  - the result

    577:    Note:
    578:    If no matrix was specified this function copies the vector.

    580:    Level: advanced

    582: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrixBV()`
    583: @*/
    584: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
    585: {
    586:   PetscFunctionBegin;
    590:   if (bv->matrix) {
    591:     PetscCall(BV_IPMatMult(bv,x));
    592:     PetscCall(VecCopy(bv->Bx,y));
    593:   } else PetscCall(VecCopy(x,y));
    594:   PetscFunctionReturn(PETSC_SUCCESS);
    595: }

    597: /*@
    598:    BVApplyMatrixBV - Multiplies the `BV` vectors by the matrix representation
    599:    of the inner product.

    601:    Neighbor-wise Collective

    603:    Input Parameter:
    604: .  X - the basis vectors context

    606:    Output Parameter:
    607: .  Y - the basis vectors to store the result (optional)

    609:    Note:
    610:    This function computes $Y = B X$, where $B$ is the matrix given with
    611:    `BVSetMatrix()`. This operation is computed as in `BVMatMult()`.
    612:    If no matrix was specified, then it just copies $Y = X$.

    614:    If no `Y` is given, the result is stored internally in the cached `BV`.

    616:    Level: developer

    618: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrix()`, `BVMatMult()`, `BVGetCachedBV()`
    619: @*/
    620: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
    621: {
    622:   PetscFunctionBegin;
    624:   if (Y) {
    626:     if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
    627:     else PetscCall(BVCopy(X,Y));
    628:   } else PetscCall(BV_IPMatMultBV(X));
    629:   PetscFunctionReturn(PETSC_SUCCESS);
    630: }

    632: /*@
    633:    BVSetSignature - Sets the signature matrix to be used in orthogonalization.

    635:    Logically Collective

    637:    Input Parameters:
    638: +  bv    - the basis vectors context
    639: -  omega - a vector representing the diagonal of the signature matrix

    641:    Note:
    642:    The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.

    644:    Level: developer

    646: .seealso: [](sec:bv), `BVSetMatrix()`, `BVGetSignature()`
    647: @*/
    648: PetscErrorCode BVSetSignature(BV bv,Vec omega)
    649: {
    650:   PetscInt          i,n;
    651:   const PetscScalar *pomega;
    652:   PetscScalar       *intern;

    654:   PetscFunctionBegin;
    656:   BVCheckSizes(bv,1);

    660:   PetscCall(VecGetSize(omega,&n));
    661:   PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
    662:   PetscCall(BV_AllocateSignature(bv));
    663:   if (bv->indef) {
    664:     PetscCall(VecGetArrayRead(omega,&pomega));
    665:     PetscCall(VecGetArray(bv->omega,&intern));
    666:     for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
    667:     PetscCall(VecRestoreArray(bv->omega,&intern));
    668:     PetscCall(VecRestoreArrayRead(omega,&pomega));
    669:   } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
    670:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
    671:   PetscFunctionReturn(PETSC_SUCCESS);
    672: }

    674: /*@
    675:    BVGetSignature - Retrieves the signature matrix from last orthogonalization.

    677:    Not Collective

    679:    Input Parameter:
    680: .  bv - the basis vectors context

    682:    Output Parameter:
    683: .  omega - a vector representing the diagonal of the signature matrix

    685:    Note:
    686:    The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.

    688:    Level: developer

    690: .seealso: [](sec:bv), `BVSetMatrix()`, `BVSetSignature()`
    691: @*/
    692: PetscErrorCode BVGetSignature(BV bv,Vec omega)
    693: {
    694:   PetscInt          i,n;
    695:   PetscScalar       *pomega;
    696:   const PetscScalar *intern;

    698:   PetscFunctionBegin;
    700:   BVCheckSizes(bv,1);

    704:   PetscCall(VecGetSize(omega,&n));
    705:   PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
    706:   if (bv->indef && bv->omega) {
    707:     PetscCall(VecGetArray(omega,&pomega));
    708:     PetscCall(VecGetArrayRead(bv->omega,&intern));
    709:     for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
    710:     PetscCall(VecRestoreArrayRead(bv->omega,&intern));
    711:     PetscCall(VecRestoreArray(omega,&pomega));
    712:   } else PetscCall(VecSet(omega,1.0));
    713:   PetscFunctionReturn(PETSC_SUCCESS);
    714: }

    716: /*@
    717:    BVSetBufferVec - Attach a vector object to be used as buffer space for
    718:    several operations.

    720:    Collective

    722:    Input Parameters:
    723: +  bv     - the basis vectors context)
    724: -  buffer - the vector

    726:    Notes:
    727:    Use `BVGetBufferVec()` to retrieve the vector (for example, to free it
    728:    at the end of the computations).

    730:    The vector must be sequential of length $(n_c+m)m$, where $m$ is the number
    731:    of columns of `bv` and $n_c$ is the number of constraints.

    733:    Level: developer

    735: .seealso: [](sec:bv), `BVGetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`
    736: @*/
    737: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
    738: {
    739:   PetscInt       ld,n;
    740:   PetscMPIInt    size;

    742:   PetscFunctionBegin;
    745:   BVCheckSizes(bv,1);
    746:   PetscCall(VecGetSize(buffer,&n));
    747:   ld = bv->m+bv->nc;
    748:   PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
    749:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
    750:   PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");

    752:   PetscCall(PetscObjectReference((PetscObject)buffer));
    753:   PetscCall(VecDestroy(&bv->buffer));
    754:   bv->buffer = buffer;
    755:   PetscFunctionReturn(PETSC_SUCCESS);
    756: }

    758: /*@
    759:    BVGetBufferVec - Obtain the buffer vector associated with the `BV` object.

    761:    Collective

    763:    Input Parameter:
    764: .  bv - the basis vectors context

    766:    Output Parameter:
    767: .  buffer - vector

    769:    Notes:
    770:    The vector is created if not available previously. It is a sequential vector
    771:    of length $(n_c+m) m$, where $m$ is the number of columns of `bv` and $n_c$ is the number
    772:    of constraints.

    774:    Developer Notes:
    775:    The buffer vector is viewed as a column-major matrix with leading dimension equal
    776:    to $n_c+m$, and $m$ columns at most. In the most common usage, it has the structure
    777:    $$
    778:      \left[\begin{array}{c|c}
    779:        s & \begin{array}{c} C \\\hline H \end{array}
    780:      \end{array}\right]
    781:    $$
    782:    where $H$ is an upper Hessenberg matrix of order $m (m-1)$, $C$ contains coefficients
    783:    related to orthogonalization against constraints (first $n_c$ rows), and $s$ is the
    784:    first column that contains scratch values computed during Gram-Schmidt
    785:    orthogonalization. In particular, `BVDotColumn()` and `BVMultColumn()` use $s$ to
    786:    store the coefficients.

    788:    Level: developer

    790: .seealso: [](sec:bv), `BVSetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`, `BVDotColumn()`, `BVMultColumn()`
    791: @*/
    792: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
    793: {
    794:   PetscInt       ld;

    796:   PetscFunctionBegin;
    798:   PetscAssertPointer(buffer,2);
    799:   BVCheckSizes(bv,1);
    800:   if (!bv->buffer) {
    801:     ld = bv->m+bv->nc;
    802:     PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
    803:     PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
    804:     PetscCall(VecSetType(bv->buffer,bv->vtype));
    805:   }
    806:   *buffer = bv->buffer;
    807:   PetscFunctionReturn(PETSC_SUCCESS);
    808: }

    810: /*@
    811:    BVSetRandomContext - Sets the `PetscRandom` object associated with the `BV`,
    812:    to be used in operations that need random numbers.

    814:    Collective

    816:    Input Parameters:
    817: +  bv   - the basis vectors context
    818: -  rand - the random number generator context

    820:    Level: advanced

    822: .seealso: [](sec:bv), `BVGetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`, `PetscRandomCreate()`
    823: @*/
    824: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
    825: {
    826:   PetscFunctionBegin;
    829:   PetscCheckSameComm(bv,1,rand,2);
    830:   PetscCall(PetscObjectReference((PetscObject)rand));
    831:   PetscCall(PetscRandomDestroy(&bv->rand));
    832:   bv->rand = rand;
    833:   PetscFunctionReturn(PETSC_SUCCESS);
    834: }

    836: /*@
    837:    BVGetRandomContext - Gets the `PetscRandom` object associated with the `BV`.

    839:    Collective

    841:    Input Parameter:
    842: .  bv - the basis vectors context

    844:    Output Parameter:
    845: .  rand - the random number generator context

    847:    Level: advanced

    849: .seealso: [](sec:bv), `BVSetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`
    850: @*/
    851: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
    852: {
    853:   PetscFunctionBegin;
    855:   PetscAssertPointer(rand,2);
    856:   if (!bv->rand) {
    857:     PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
    858:     if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
    859:     if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
    860:     if (bv->rrandom) {
    861:       PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
    862:       PetscCall(PetscRandomSeed(bv->rand));
    863:     }
    864:   }
    865:   *rand = bv->rand;
    866:   PetscFunctionReturn(PETSC_SUCCESS);
    867: }

    869: /*@
    870:    BVSetFromOptions - Sets `BV` options from the options database.

    872:    Collective

    874:    Input Parameter:
    875: .  bv - the basis vectors context

    877:    Note:
    878:    To see all options, run your program with the `-help` option.

    880:    Level: beginner

    882: .seealso: [](sec:bv), `BVSetOptionsPrefix()`
    883: @*/
    884: PetscErrorCode BVSetFromOptions(BV bv)
    885: {
    886:   char               type[256];
    887:   PetscBool          flg1,flg2,flg3,flg4;
    888:   PetscReal          r;
    889:   BVOrthogType       otype;
    890:   BVOrthogRefineType orefine;
    891:   BVOrthogBlockType  oblock;

    893:   PetscFunctionBegin;
    895:   PetscCall(BVRegisterAll());
    896:   PetscObjectOptionsBegin((PetscObject)bv);
    897:     PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
    898:     if (flg1) PetscCall(BVSetType(bv,type));
    899:     else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));

    901:     otype = bv->orthog_type;
    902:     PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
    903:     orefine = bv->orthog_ref;
    904:     PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
    905:     r = bv->orthog_eta;
    906:     PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
    907:     oblock = bv->orthog_block;
    908:     PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
    909:     if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));

    911:     PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));

    913:     PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
    914:     if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));

    916:     /* undocumented option to generate random vectors that are independent of the number of processes */
    917:     PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));

    919:     if (bv->ops->create) bv->defersfo = PETSC_TRUE;   /* defer call to setfromoptions */
    920:     else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
    921:     PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
    922:   PetscOptionsEnd();
    923:   bv->sfocalled = PETSC_TRUE;
    924:   PetscFunctionReturn(PETSC_SUCCESS);
    925: }

    927: /*@
    928:    BVSetOrthogonalization - Specifies the method used for the orthogonalization of
    929:    vectors (classical or modified Gram-Schmidt with or without refinement), and
    930:    for the block-orthogonalization (simultaneous orthogonalization of a set of
    931:    vectors).

    933:    Logically Collective

    935:    Input Parameters:
    936: +  bv     - the basis vectors context
    937: .  type   - the method of vector orthogonalization
    938: .  refine - type of refinement
    939: .  eta    - parameter for selective refinement
    940: -  block  - the method of block orthogonalization

    942:    Options Database Keys:
    943: +  -bv_orthog_type (cgs|mgs)                     - the vector orthogonalization
    944: .  -bv_orthog_refine (never|ifneeded|always)     - the refinement type, default is `ifneeded`
    945: .  -bv_orthog_eta eta                            - the value of `eta`
    946: -  -bv_orthog_block (gs|chol|tsqr|tsqrchol|svqb) - the block-orthogonalization method

    948:    Notes:
    949:    The default settings work well for most problems.

    951:    The parameter `eta` should be a real value between 0 and 1, that is used only when
    952:    the refinement type is `ifneeded`. Use `PETSC_DETERMINE` to set a reasonable
    953:    default value. Use `PETSC_CURRENT` to leave the current value unchanged.

    955:    When using several processes, `BV_ORTHOG_MGS` is likely to result in bad scalability.

    957:    If the method set for block orthogonalization is `BV_ORTHOG_BLOCK_GS`, then the computation
    958:    is done column by column with the vector orthogonalization.

    960:    Level: advanced

    962: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVGetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
    963: @*/
    964: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
    965: {
    966:   PetscFunctionBegin;
    972:   switch (type) {
    973:     case BV_ORTHOG_CGS:
    974:     case BV_ORTHOG_MGS:
    975:       bv->orthog_type = type;
    976:       break;
    977:     default:
    978:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
    979:   }
    980:   switch (refine) {
    981:     case BV_ORTHOG_REFINE_NEVER:
    982:     case BV_ORTHOG_REFINE_IFNEEDED:
    983:     case BV_ORTHOG_REFINE_ALWAYS:
    984:       bv->orthog_ref = refine;
    985:       break;
    986:     default:
    987:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
    988:   }
    989:   if (eta == (PetscReal)PETSC_DETERMINE) {
    990:     bv->orthog_eta = 0.7071;
    991:   } else if (eta != (PetscReal)PETSC_CURRENT) {
    992:     PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
    993:     bv->orthog_eta = eta;
    994:   }
    995:   switch (block) {
    996:     case BV_ORTHOG_BLOCK_GS:
    997:     case BV_ORTHOG_BLOCK_CHOL:
    998:     case BV_ORTHOG_BLOCK_TSQR:
    999:     case BV_ORTHOG_BLOCK_TSQRCHOL:
    1000:     case BV_ORTHOG_BLOCK_SVQB:
    1001:       bv->orthog_block = block;
    1002:       break;
    1003:     default:
    1004:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
    1005:   }
    1006:   PetscFunctionReturn(PETSC_SUCCESS);
    1007: }

    1009: /*@
    1010:    BVGetOrthogonalization - Gets the orthogonalization settings from the `BV` object.

    1012:    Not Collective

    1014:    Input Parameter:
    1015: .  bv - basis vectors context

    1017:    Output Parameters:
    1018: +  type   - the method of vector orthogonalization
    1019: .  refine - type of refinement
    1020: .  eta    - parameter for selective refinement
    1021: -  block  - the method of block orthogonalization

    1023:    Level: advanced

    1025: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVSetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
    1026: @*/
    1027: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
    1028: {
    1029:   PetscFunctionBegin;
    1031:   if (type)   *type   = bv->orthog_type;
    1032:   if (refine) *refine = bv->orthog_ref;
    1033:   if (eta)    *eta    = bv->orthog_eta;
    1034:   if (block)  *block  = bv->orthog_block;
    1035:   PetscFunctionReturn(PETSC_SUCCESS);
    1036: }

    1038: /*@
    1039:    BVSetMatMultMethod - Specifies the method used for the `BVMatMult()` operation.

    1041:    Logically Collective

    1043:    Input Parameters:
    1044: +  bv     - the basis vectors context
    1045: -  method - the method for the `BVMatMult()` operation

    1047:    Options Database Key:
    1048: .  -bv_matmult (vecs|mat) - choose the method

    1050:    Note:
    1051:    The default is `BV_MATMULT_MAT` except in the case of `BVVECS`.

    1053:    Level: advanced

    1055: .seealso: [](sec:bv), `BVMatMult()`, `BVGetMatMultMethod()`, `BVMatMultType`
    1056: @*/
    1057: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
    1058: {
    1059:   PetscFunctionBegin;
    1062:   switch (method) {
    1063:     case BV_MATMULT_VECS:
    1064:     case BV_MATMULT_MAT:
    1065:       bv->vmm = method;
    1066:       break;
    1067:     case BV_MATMULT_MAT_SAVE:
    1068:       PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
    1069:       bv->vmm = BV_MATMULT_MAT;
    1070:       break;
    1071:     default:
    1072:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
    1073:   }
    1074:   PetscFunctionReturn(PETSC_SUCCESS);
    1075: }

    1077: /*@
    1078:    BVGetMatMultMethod - Gets the method used for the `BVMatMult()` operation.

    1080:    Not Collective

    1082:    Input Parameter:
    1083: .  bv - basis vectors context

    1085:    Output Parameter:
    1086: .  method - the method for the `BVMatMult()` operation

    1088:    Level: advanced

    1090: .seealso: [](sec:bv), `BVMatMult()`, `BVSetMatMultMethod()`, `BVMatMultType`
    1091: @*/
    1092: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
    1093: {
    1094:   PetscFunctionBegin;
    1096:   PetscAssertPointer(method,2);
    1097:   *method = bv->vmm;
    1098:   PetscFunctionReturn(PETSC_SUCCESS);
    1099: }

    1101: /*@
    1102:    BVGetColumn - Returns a `Vec` object that contains the entries of the
    1103:    requested column of the basis vectors object.

    1105:    Logically Collective

    1107:    Input Parameters:
    1108: +  bv - the basis vectors context
    1109: -  j  - the index of the requested column

    1111:    Output Parameter:
    1112: .  v  - vector containing the jth column

    1114:    Notes:
    1115:    The returned `Vec` must be seen as a reference (not a copy) of the `BV`
    1116:    column, that is, modifying the `Vec` will change the `BV` entries as well.

    1118:    The returned `Vec` must not be destroyed. `BVRestoreColumn()` must be
    1119:    called when it is no longer needed. At most, two columns can be fetched,
    1120:    that is, this function can only be called twice before the corresponding
    1121:    `BVRestoreColumn()` is invoked.

    1123:    A negative index `j` selects the `i`-th constraint, where `i=-j`. Constraints
    1124:    should not be modified.

    1126:    Level: beginner

    1128: .seealso: [](sec:bv), `BVRestoreColumn()`, `BVInsertConstraints()`
    1129: @*/
    1130: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
    1131: {
    1132:   PetscInt       l;

    1134:   PetscFunctionBegin;
    1137:   BVCheckSizes(bv,1);
    1138:   BVCheckOp(bv,1,getcolumn);
    1140:   PetscCheck(j>=0 || -j<=bv->nc,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %" PetscInt_FMT " but only %" PetscInt_FMT " are available",-j,bv->nc);
    1141:   PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %" PetscInt_FMT " but only %" PetscInt_FMT " are available",j,bv->m);
    1142:   PetscCheck(j!=bv->ci[0] && j!=bv->ci[1],PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Column %" PetscInt_FMT " already fetched in a previous call to BVGetColumn",j);
    1143:   l = BVAvailableVec;
    1144:   PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
    1145:   PetscUseTypeMethod(bv,getcolumn,j,v);
    1146:   bv->ci[l] = j;
    1147:   PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
    1148:   PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
    1149:   *v = bv->cv[l];
    1150:   PetscFunctionReturn(PETSC_SUCCESS);
    1151: }

    1153: /*@
    1154:    BVRestoreColumn - Restore a column obtained with `BVGetColumn()`.

    1156:    Logically Collective

    1158:    Input Parameters:
    1159: +  bv - the basis vectors context
    1160: .  j  - the index of the column
    1161: -  v  - vector obtained with `BVGetColumn()`

    1163:    Note:
    1164:    The arguments must match the corresponding call to `BVGetColumn()`.

    1166:    Level: beginner

    1168: .seealso: [](sec:bv), `BVGetColumn()`
    1169: @*/
    1170: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
    1171: {
    1172:   PetscObjectId    id;
    1173:   PetscObjectState st;
    1174:   PetscInt         l;

    1176:   PetscFunctionBegin;
    1179:   BVCheckSizes(bv,1);
    1181:   PetscAssertPointer(v,3);
    1183:   PetscCheck(j>=0 || -j<=bv->nc,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %" PetscInt_FMT " but only %" PetscInt_FMT " are available",-j,bv->nc);
    1184:   PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %" PetscInt_FMT " but only %" PetscInt_FMT " are available",j,bv->m);
    1185:   PetscCheck(j==bv->ci[0] || j==bv->ci[1],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Column %" PetscInt_FMT " has not been fetched with a call to BVGetColumn",j);
    1186:   l = (j==bv->ci[0])? 0: 1;
    1187:   PetscCall(PetscObjectGetId((PetscObject)*v,&id));
    1188:   PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
    1189:   PetscCall(VecGetState(*v,&st));
    1190:   if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
    1191:   PetscUseTypeMethod(bv,restorecolumn,j,v);
    1192:   bv->ci[l] = -bv->nc-1;
    1193:   bv->st[l] = -1;
    1194:   bv->id[l] = 0;
    1195:   *v = NULL;
    1196:   PetscFunctionReturn(PETSC_SUCCESS);
    1197: }

    1199: /*@C
    1200:    BVGetArray - Returns a pointer to a contiguous array that contains this
    1201:    process' portion of the `BV` data.

    1203:    Logically Collective

    1205:    Input Parameter:
    1206: .  bv - the basis vectors context

    1208:    Output Parameter:
    1209: .  a  - location to put pointer to the array

    1211:    Notes:
    1212:    `BVRestoreArray()` must be called when access to the array is no longer needed.
    1213:    This operation may imply a data copy, for `BV` types that do not store
    1214:    data contiguously in memory.

    1216:    The pointer will normally point to the first entry of the first column,
    1217:    but if the `BV` has constraints then these go before the regular columns.

    1219:    Note that for manipulating the pointer to the `BV` array, one must take into
    1220:    account the leading dimension, which might be different from the local
    1221:    number of rows, see `BVGetLeadingDimension()`.

    1223:    Use `BVGetArrayRead()` for read-only access.

    1225:    Level: advanced

    1227: .seealso: [](sec:bv), `BVRestoreArray()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArrayRead()`, `BVType`
    1228: @*/
    1229: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
    1230: {
    1231:   PetscFunctionBegin;
    1234:   BVCheckSizes(bv,1);
    1235:   BVCheckOp(bv,1,getarray);
    1236:   PetscUseTypeMethod(bv,getarray,a);
    1237:   PetscFunctionReturn(PETSC_SUCCESS);
    1238: }

    1240: /*@C
    1241:    BVRestoreArray - Restore the `BV` object after `BVGetArray()` has been called.

    1243:    Logically Collective

    1245:    Input Parameters:
    1246: +  bv - the basis vectors context
    1247: -  a  - location of pointer to array obtained from `BVGetArray()`

    1249:    Note:
    1250:    This operation may imply a data copy, for `BV` types that do not store
    1251:    data contiguously in memory.

    1253:    Level: advanced

    1255: .seealso: [](sec:bv), `BVGetArray()`
    1256: @*/
    1257: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
    1258: {
    1259:   PetscFunctionBegin;
    1262:   BVCheckSizes(bv,1);
    1263:   PetscTryTypeMethod(bv,restorearray,a);
    1264:   if (a) *a = NULL;
    1265:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
    1266:   PetscFunctionReturn(PETSC_SUCCESS);
    1267: }

    1269: /*@C
    1270:    BVGetArrayRead - Returns a read-only pointer to a contiguous array that
    1271:    contains this processor's portion of the `BV` data.

    1273:    Not Collective

    1275:    Input Parameter:
    1276: .  bv - the basis vectors context

    1278:    Output Parameter:
    1279: .  a  - location to put pointer to the array

    1281:    Notes:
    1282:    `BVRestoreArrayRead()` must be called when access to the array is no
    1283:    longer needed. This operation may imply a data copy, for `BV` types that
    1284:    do not store data contiguously in memory.

    1286:    The pointer will normally point to the first entry of the first column,
    1287:    but if the `BV` has constraints then these go before the regular columns.

    1289:    Level: advanced

    1291: .seealso: [](sec:bv), `BVRestoreArrayRead()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArray()`, `BVType`
    1292: @*/
    1293: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
    1294: {
    1295:   PetscFunctionBegin;
    1298:   BVCheckSizes(bv,1);
    1299:   BVCheckOp(bv,1,getarrayread);
    1300:   PetscUseTypeMethod(bv,getarrayread,a);
    1301:   PetscFunctionReturn(PETSC_SUCCESS);
    1302: }

    1304: /*@C
    1305:    BVRestoreArrayRead - Restore the `BV` object after `BVGetArrayRead()` has
    1306:    been called.

    1308:    Not Collective

    1310:    Input Parameters:
    1311: +  bv - the basis vectors context
    1312: -  a  - location of pointer to array obtained from `BVGetArrayRead()`

    1314:    Level: advanced

    1316: .seealso: [](sec:bv), `BVGetArrayRead()`
    1317: @*/
    1318: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
    1319: {
    1320:   PetscFunctionBegin;
    1323:   BVCheckSizes(bv,1);
    1324:   PetscTryTypeMethod(bv,restorearrayread,a);
    1325:   if (a) *a = NULL;
    1326:   PetscFunctionReturn(PETSC_SUCCESS);
    1327: }

    1329: /*@
    1330:    BVCreateVec - Creates a new `Vec` object with the same type and dimensions
    1331:    as the columns of the basis vectors object.

    1333:    Collective

    1335:    Input Parameter:
    1336: .  bv - the basis vectors context

    1338:    Output Parameter:
    1339: .  v  - the new vector

    1341:    Note:
    1342:    The user is responsible for destroying the returned vector.

    1344:    Level: beginner

    1346: .seealso: [](sec:bv), `BVCreateMat()`, `BVCreateVecEmpty()`
    1347: @*/
    1348: PetscErrorCode BVCreateVec(BV bv,Vec *v)
    1349: {
    1350:   PetscFunctionBegin;
    1352:   BVCheckSizes(bv,1);
    1353:   PetscAssertPointer(v,2);
    1354:   PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
    1355:   PetscCall(VecSetLayout(*v,bv->map));
    1356:   PetscCall(VecSetType(*v,bv->vtype));
    1357:   PetscCall(VecSetUp(*v));
    1358:   PetscFunctionReturn(PETSC_SUCCESS);
    1359: }

    1361: /*@
    1362:    BVCreateVecEmpty - Creates a new `Vec` object with the same type and dimensions
    1363:    as the columns of the basis vectors object, but without internal array.

    1365:    Collective

    1367:    Input Parameter:
    1368: .  bv - the basis vectors context

    1370:    Output Parameter:
    1371: .  v  - the new vector

    1373:    Note:
    1374:    This works as `BVCreateVec()`, but the new vector does not have the array allocated,
    1375:    so the intended usage is with `VecPlaceArray()`.

    1377:    Level: developer

    1379: .seealso: [](sec:bv), `BVCreateVec()`
    1380: @*/
    1381: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
    1382: {
    1383:   PetscBool  standard,cuda,hip,mpi;
    1384:   PetscInt   N,nloc,bs;

    1386:   PetscFunctionBegin;
    1388:   BVCheckSizes(bv,1);
    1389:   PetscAssertPointer(v,2);

    1391:   PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
    1392:   PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
    1393:   PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
    1394:   if (standard || cuda || hip) {
    1395:     PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
    1396:     PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
    1397:     PetscCall(PetscLayoutGetSize(bv->map,&N));
    1398:     PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
    1399:     if (cuda) {
    1400: #if defined(PETSC_HAVE_CUDA)
    1401:       if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
    1402:       else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
    1403: #endif
    1404:     } else if (hip) {
    1405: #if defined(PETSC_HAVE_HIP)
    1406:       if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
    1407:       else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
    1408: #endif
    1409:     } else {
    1410:       if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
    1411:       else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
    1412:     }
    1413:   } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
    1414:   PetscFunctionReturn(PETSC_SUCCESS);
    1415: }

    1417: /*@
    1418:    BVSetVecType - Set the vector type to be used when creating vectors via `BVCreateVec()`.

    1420:    Collective

    1422:    Input Parameters:
    1423: +  bv - the basis vectors context
    1424: -  vtype - the vector type

    1426:    Level: advanced

    1428:    Note:
    1429:    This is not needed if the `BV` object is set up with `BVSetSizesFromVec()`, but may be
    1430:    required in the case of `BVSetSizes()` if one wants to work with non-standard vectors.

    1432: .seealso: [](sec:bv), `BVCreateVec()`, `BVGetVecType()`, `BVSetSizesFromVec()`, `BVSetSizes()`
    1433: @*/
    1434: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
    1435: {
    1436:   PetscBool   std;
    1437:   PetscMPIInt size;

    1439:   PetscFunctionBegin;
    1441:   PetscCall(PetscFree(bv->vtype));
    1442:   PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
    1443:   if (std) {
    1444:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
    1445:     PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
    1446:   } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
    1447:   PetscFunctionReturn(PETSC_SUCCESS);
    1448: }

    1450: /*@
    1451:    BVGetVecType - Get the vector type to be used when creating vectors via `BVCreateVec()`.

    1453:    Not Collective

    1455:    Input Parameter:
    1456: .  bv - the basis vectors context

    1458:    Output Parameter:
    1459: .  vtype - the vector type

    1461:    Level: advanced

    1463: .seealso: [](sec:bv), `BVSetVecType()`, `BVCreateVec()`
    1464: @*/
    1465: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
    1466: {
    1467:   PetscFunctionBegin;
    1469:   PetscAssertPointer(vtype,2);
    1470:   *vtype = bv->vtype;
    1471:   PetscFunctionReturn(PETSC_SUCCESS);
    1472: }

    1474: /*@
    1475:    BVCreateMat - Creates a new `Mat` object of dense type and copies the contents
    1476:    of the `BV` object.

    1478:    Collective

    1480:    Input Parameter:
    1481: .  bv - the basis vectors context

    1483:    Output Parameter:
    1484: .  A  - the new matrix

    1486:    Notes:
    1487:    The user is responsible for destroying the returned matrix.

    1489:    The matrix contains all columns of the `BV`, not just the active columns.

    1491:    Level: intermediate

    1493: .seealso: [](sec:bv), `BVCreateFromMat()`, `BVCreateVec()`, `BVGetMat()`
    1494: @*/
    1495: PetscErrorCode BVCreateMat(BV bv,Mat *A)
    1496: {
    1497:   PetscInt ksave,lsave;
    1498:   Mat      B;

    1500:   PetscFunctionBegin;
    1502:   BVCheckSizes(bv,1);
    1503:   PetscAssertPointer(A,2);

    1505:   PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
    1506:   lsave = bv->l;
    1507:   ksave = bv->k;
    1508:   bv->l = 0;
    1509:   bv->k = bv->m;
    1510:   PetscCall(BVGetMat(bv,&B));
    1511:   PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
    1512:   PetscCall(BVRestoreMat(bv,&B));
    1513:   bv->l = lsave;
    1514:   bv->k = ksave;
    1515:   PetscFunctionReturn(PETSC_SUCCESS);
    1516: }

    1518: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
    1519: {
    1520:   PetscScalar *vv,*aa;
    1521:   PetscBool   create=PETSC_FALSE;
    1522:   PetscInt    m,cols;

    1524:   PetscFunctionBegin;
    1525:   m = bv->k-bv->l;
    1526:   if (!bv->Aget) create=PETSC_TRUE;
    1527:   else {
    1528:     PetscCall(MatDenseGetArray(bv->Aget,&aa));
    1529:     PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
    1530:     PetscCall(MatGetSize(bv->Aget,NULL,&cols));
    1531:     if (cols!=m) {
    1532:       PetscCall(MatDestroy(&bv->Aget));
    1533:       create=PETSC_TRUE;
    1534:     }
    1535:   }
    1536:   PetscCall(BVGetArray(bv,&vv));
    1537:   if (create) {
    1538:     PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,m,bv->ld,vv,&bv->Aget)); /* pass a pointer to avoid allocation of storage */
    1539:     PetscCall(MatDenseReplaceArray(bv->Aget,NULL));  /* replace with a null pointer, the value after BVRestoreMat */
    1540:   }
    1541:   PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld));  /* set the actual pointer */
    1542:   *A = bv->Aget;
    1543:   PetscFunctionReturn(PETSC_SUCCESS);
    1544: }

    1546: /*@
    1547:    BVGetMat - Returns a `Mat` object of dense type that shares the memory of
    1548:    the `BV` object.

    1550:    Collective

    1552:    Input Parameter:
    1553: .  bv - the basis vectors context

    1555:    Output Parameter:
    1556: .  A  - the matrix

    1558:    Notes:
    1559:    The returned matrix contains only the active columns. If the content of
    1560:    the `Mat` is modified, these changes are also done in the `BV` object. The
    1561:    user must call `BVRestoreMat()` when no longer needed.

    1563:    This operation implies a call to `BVGetArray()`, which may result in data
    1564:    copies.

    1566:    Level: advanced

    1568: .seealso: [](sec:bv), `BVRestoreMat()`, `BVCreateMat()`, `BVGetArray()`
    1569: @*/
    1570: PetscErrorCode BVGetMat(BV bv,Mat *A)
    1571: {
    1572:   char name[64];

    1574:   PetscFunctionBegin;
    1576:   BVCheckSizes(bv,1);
    1577:   PetscAssertPointer(A,2);
    1578:   PetscUseTypeMethod(bv,getmat,A);
    1579:   if (((PetscObject)bv)->name) {   /* set A's name based on BV name */
    1580:     PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
    1581:     PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
    1582:     PetscCall(PetscObjectSetName((PetscObject)*A,name));
    1583:   }
    1584:   PetscFunctionReturn(PETSC_SUCCESS);
    1585: }

    1587: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
    1588: {
    1589:   PetscScalar *vv,*aa;

    1591:   PetscFunctionBegin;
    1592:   PetscCall(MatDenseGetArray(bv->Aget,&aa));
    1593:   vv = aa-(bv->nc+bv->l)*bv->ld;
    1594:   PetscCall(MatDenseResetArray(bv->Aget));
    1595:   PetscCall(BVRestoreArray(bv,&vv));
    1596:   *A = NULL;
    1597:   PetscFunctionReturn(PETSC_SUCCESS);
    1598: }

    1600: /*@
    1601:    BVRestoreMat - Restores the Mat obtained with `BVGetMat()`.

    1603:    Logically Collective

    1605:    Input Parameters:
    1606: +  bv - the basis vectors context
    1607: -  A  - the fetched matrix

    1609:    Note:
    1610:    A call to this function must match a previous call of `BVGetMat()`.
    1611:    The effect is that the contents of the `Mat` are copied back to the
    1612:    `BV` internal data structures.

    1614:    Level: advanced

    1616: .seealso: [](sec:bv), `BVGetMat()`, `BVRestoreArray()`
    1617: @*/
    1618: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
    1619: {
    1620:   PetscFunctionBegin;
    1622:   BVCheckSizes(bv,1);
    1623:   PetscAssertPointer(A,2);
    1624:   PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
    1625:   PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
    1626:   PetscUseTypeMethod(bv,restoremat,A);
    1627:   PetscFunctionReturn(PETSC_SUCCESS);
    1628: }

    1630: /*
    1631:    Copy all user-provided attributes of V to another BV object W
    1632:  */
    1633: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
    1634: {
    1635:   PetscFunctionBegin;
    1636:   PetscCall(PetscLayoutReference(V->map,&W->map));
    1637:   PetscCall(BVSetVecType(W,V->vtype));
    1638:   W->ld           = V->ld;
    1639:   PetscCall(BVSetType(W,((PetscObject)V)->type_name));
    1640:   W->orthog_type  = V->orthog_type;
    1641:   W->orthog_ref   = V->orthog_ref;
    1642:   W->orthog_eta   = V->orthog_eta;
    1643:   W->orthog_block = V->orthog_block;
    1644:   PetscCall(PetscObjectReference((PetscObject)V->matrix));
    1645:   W->matrix       = V->matrix;
    1646:   W->indef        = V->indef;
    1647:   W->vmm          = V->vmm;
    1648:   W->rrandom      = V->rrandom;
    1649:   W->deftol       = V->deftol;
    1650:   PetscCall(PetscObjectReference((PetscObject)V->rand));
    1651:   W->rand         = V->rand;
    1652:   W->sfocalled    = V->sfocalled;
    1653:   PetscTryTypeMethod(V,duplicate,W);
    1654:   PetscCall(PetscObjectStateIncrease((PetscObject)W));
    1655:   PetscFunctionReturn(PETSC_SUCCESS);
    1656: }

    1658: /*@
    1659:    BVDuplicate - Creates a new basis vectors object of the same type and
    1660:    dimensions as an existing one.

    1662:    Collective

    1664:    Input Parameter:
    1665: .  V - basis vectors context

    1667:    Output Parameter:
    1668: .  W - location to put the new `BV`

    1670:    Notes:
    1671:    The new `BV` has the same type and dimensions as `V`. Also, the inner
    1672:    product matrix and orthogonalization options are copied.

    1674:    `BVDuplicate()` DOES NOT COPY the entries, but rather allocates storage
    1675:    for the new basis vectors. Use `BVCopy()` to copy the content.

    1677:    Level: intermediate

    1679: .seealso: [](sec:bv), `BVDuplicateResize()`, `BVCreate()`, `BVSetSizesFromVec()`, `BVCopy()`
    1680: @*/
    1681: PetscErrorCode BVDuplicate(BV V,BV *W)
    1682: {
    1683:   PetscFunctionBegin;
    1686:   BVCheckSizes(V,1);
    1687:   PetscAssertPointer(W,2);
    1688:   PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
    1689:   (*W)->N = V->N;
    1690:   (*W)->n = V->n;
    1691:   (*W)->m = V->m;
    1692:   (*W)->k = V->m;
    1693:   PetscCall(BVDuplicate_Private(V,*W));
    1694:   PetscFunctionReturn(PETSC_SUCCESS);
    1695: }

    1697: /*@
    1698:    BVDuplicateResize - Creates a new basis vectors object of the same type and
    1699:    dimensions as an existing one, but with possibly different number of columns.

    1701:    Collective

    1703:    Input Parameters:
    1704: +  V - basis vectors context
    1705: -  m - the new number of columns

    1707:    Output Parameter:
    1708: .  W - location to put the new BV

    1710:    Note:
    1711:    This is equivalent to a call to `BVDuplicate()` followed by `BVResize()`. The
    1712:    contents of `V` are not copied to `W`.

    1714:    Level: intermediate

    1716: .seealso: [](sec:bv), `BVDuplicate()`, `BVResize()`
    1717: @*/
    1718: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
    1719: {
    1720:   PetscFunctionBegin;
    1723:   BVCheckSizes(V,1);
    1725:   PetscAssertPointer(W,3);
    1726:   PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
    1727:   (*W)->N = V->N;
    1728:   (*W)->n = V->n;
    1729:   (*W)->m = m;
    1730:   (*W)->k = m;
    1731:   PetscCall(BVDuplicate_Private(V,*W));
    1732:   PetscFunctionReturn(PETSC_SUCCESS);
    1733: }

    1735: /*@
    1736:    BVGetCachedBV - Returns a `BV` object stored internally that holds the
    1737:    result of $BX$ after a call to `BVApplyMatrixBV()`.

    1739:    Collective

    1741:    Input Parameter:
    1742: .  bv    - the basis vectors context

    1744:    Output Parameter:
    1745: .  cached - the cached `BV`

    1747:    Note:
    1748:    The cached `BV` is created if not available previously.

    1750:    Level: developer

    1752: .seealso: [](sec:bv), `BVApplyMatrixBV()`
    1753: @*/
    1754: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
    1755: {
    1756:   PetscFunctionBegin;
    1758:   PetscAssertPointer(cached,2);
    1759:   BVCheckSizes(bv,1);
    1760:   if (!bv->cached) {
    1761:     PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
    1762:     bv->cached->N = bv->N;
    1763:     bv->cached->n = bv->n;
    1764:     bv->cached->m = bv->m;
    1765:     bv->cached->k = bv->m;
    1766:     PetscCall(BVDuplicate_Private(bv,bv->cached));
    1767:   }
    1768:   *cached = bv->cached;
    1769:   PetscFunctionReturn(PETSC_SUCCESS);
    1770: }

    1772: /*@
    1773:    BVCopy - Copies a basis vector object into another one, $W \leftarrow V$.

    1775:    Logically Collective

    1777:    Input Parameter:
    1778: .  V - basis vectors context

    1780:    Output Parameter:
    1781: .  W - the copy

    1783:    Note:
    1784:    Both `V` and `W` must be distributed in the same manner; local copies are
    1785:    done. Only active columns (excluding the leading ones) are copied.
    1786:    In the destination `W`, columns are overwritten starting from the leading ones.
    1787:    Constraints are not copied.

    1789:    Level: beginner

    1791: .seealso: [](sec:bv), `BVCopyVec()`, `BVCopyColumn()`, `BVDuplicate()`, `BVSetActiveColumns()`
    1792: @*/
    1793: PetscErrorCode BVCopy(BV V,BV W)
    1794: {
    1795:   PetscScalar       *womega;
    1796:   const PetscScalar *vomega;

    1798:   PetscFunctionBegin;
    1801:   BVCheckSizes(V,1);
    1802:   BVCheckOp(V,1,copy);
    1805:   BVCheckSizes(W,2);
    1806:   PetscCheckSameTypeAndComm(V,1,W,2);
    1807:   PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
    1808:   PetscCheck(V->k-V->l<=W->m-W->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"W has %" PetscInt_FMT " non-leading columns, not enough to store %" PetscInt_FMT " columns",W->m-W->l,V->k-V->l);
    1809:   if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);

    1811:   PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
    1812:   if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
    1813:     /* copy signature */
    1814:     PetscCall(BV_AllocateSignature(W));
    1815:     PetscCall(VecGetArrayRead(V->omega,&vomega));
    1816:     PetscCall(VecGetArray(W->omega,&womega));
    1817:     PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
    1818:     PetscCall(VecRestoreArray(W->omega,&womega));
    1819:     PetscCall(VecRestoreArrayRead(V->omega,&vomega));
    1820:   }
    1821:   PetscUseTypeMethod(V,copy,W);
    1822:   PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
    1823:   PetscCall(PetscObjectStateIncrease((PetscObject)W));
    1824:   PetscFunctionReturn(PETSC_SUCCESS);
    1825: }

    1827: /*@
    1828:    BVCopyVec - Copies one of the columns of a basis vectors object into a `Vec`.

    1830:    Logically Collective

    1832:    Input Parameters:
    1833: +  V - basis vectors context
    1834: -  j - the index of the column to be copied

    1836:    Output Parameter:
    1837: .  w - the copied column

    1839:    Note:
    1840:    Both `V` and `w` must be distributed in the same manner; local copies are done.

    1842:    Level: beginner

    1844: .seealso: [](sec:bv), `BVCopy()`, `BVCopyColumn()`, `BVInsertVec()`
    1845: @*/
    1846: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
    1847: {
    1848:   PetscInt       n,N;
    1849:   Vec            z;

    1851:   PetscFunctionBegin;
    1854:   BVCheckSizes(V,1);
    1857:   PetscCheckSameComm(V,1,w,3);

    1859:   PetscCall(VecGetSize(w,&N));
    1860:   PetscCall(VecGetLocalSize(w,&n));
    1861:   PetscCheck(N==V->N && n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ") do not match BV sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ")",N,n,V->N,V->n);

    1863:   PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
    1864:   PetscCall(BVGetColumn(V,j,&z));
    1865:   PetscCall(VecCopy(z,w));
    1866:   PetscCall(BVRestoreColumn(V,j,&z));
    1867:   PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
    1868:   PetscFunctionReturn(PETSC_SUCCESS);
    1869: }

    1871: /*@
    1872:    BVCopyColumn - Copies the values from one of the columns to another one.

    1874:    Logically Collective

    1876:    Input Parameters:
    1877: +  V - basis vectors context
    1878: .  j - the index of the source column
    1879: -  i - the index of the destination column

    1881:    Level: beginner

    1883: .seealso: [](sec:bv), `BVCopy()`, `BVCopyVec()`
    1884: @*/
    1885: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
    1886: {
    1887:   PetscScalar *omega;

    1889:   PetscFunctionBegin;
    1892:   BVCheckSizes(V,1);
    1895:   if (j==i) PetscFunctionReturn(PETSC_SUCCESS);

    1897:   PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
    1898:   if (V->omega) {
    1899:     PetscCall(VecGetArray(V->omega,&omega));
    1900:     omega[i] = omega[j];
    1901:     PetscCall(VecRestoreArray(V->omega,&omega));
    1902:   }
    1903:   PetscUseTypeMethod(V,copycolumn,j,i);
    1904:   PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
    1905:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
    1906:   PetscFunctionReturn(PETSC_SUCCESS);
    1907: }

    1909: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
    1910: {
    1911:   PetscInt  ncols;

    1913:   PetscFunctionBegin;
    1914:   ncols = left? bv->nc+bv->l: bv->m-bv->l;
    1915:   if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
    1916:   if (!*split) {
    1917:     PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
    1918:     (*split)->issplit = left? 1: 2;
    1919:     (*split)->splitparent = bv;
    1920:     (*split)->N = bv->N;
    1921:     (*split)->n = bv->n;
    1922:     (*split)->m = bv->m;
    1923:     (*split)->k = bv->m;
    1924:     PetscCall(BVDuplicate_Private(bv,*split));
    1925:   }
    1926:   (*split)->l  = 0;
    1927:   (*split)->k  = left? bv->l: bv->k-bv->l;
    1928:   (*split)->nc = left? bv->nc: 0;
    1929:   (*split)->m  = ncols-(*split)->nc;
    1930:   if ((*split)->nc) {
    1931:     (*split)->ci[0] = -(*split)->nc-1;
    1932:     (*split)->ci[1] = -(*split)->nc-1;
    1933:   }
    1934:   if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
    1935:   else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
    1936:   PetscFunctionReturn(PETSC_SUCCESS);
    1937: }

    1939: /*@
    1940:    BVGetSplit - Splits the `BV` object into two `BV` objects that share the
    1941:    internal data, one of them containing the leading columns and the other
    1942:    one containing the remaining columns.

    1944:    Collective

    1946:    Input Parameter:
    1947: .  bv - the basis vectors context

    1949:    Output Parameters:
    1950: +  L - left `BV` containing leading columns (can be `NULL`)
    1951: -  R - right `BV` containing remaining columns (can be `NULL`)

    1953:    Notes:
    1954:    The columns are split in two sets. The leading columns (including the
    1955:    constraints) are assigned to the left `BV` and the remaining columns
    1956:    are assigned to the right `BV`. The number of leading columns, as
    1957:    specified with `BVSetActiveColumns()`, must be between `1` and `m-1` (to
    1958:    guarantee that both `L` and `R` have at least one column).

    1960:    The returned `BV`s must be seen as references (not copies) of the input
    1961:    `BV`, that is, modifying them will change the entries of `bv` as well.
    1962:    The returned `BV`s must not be destroyed. `BVRestoreSplit()` must be called
    1963:    when they are no longer needed.

    1965:    Pass `NULL` for any of the output `BV`s that is not needed.

    1967:    Level: advanced

    1969: .seealso: [](sec:bv), `BVRestoreSplit()`, `BVSetActiveColumns()`, `BVSetNumConstraints()`, `BVGetSplitRows()`
    1970: @*/
    1971: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
    1972: {
    1973:   PetscFunctionBegin;
    1976:   BVCheckSizes(bv,1);
    1977:   PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
    1978:   PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
    1979:   PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
    1980:   bv->lsplit = bv->nc+bv->l;
    1981:   PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
    1982:   PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
    1983:   if (L) *L = bv->L;
    1984:   if (R) *R = bv->R;
    1985:   PetscFunctionReturn(PETSC_SUCCESS);
    1986: }

    1988: /*@
    1989:    BVRestoreSplit - Restore the BV objects obtained with `BVGetSplit()`.

    1991:    Logically Collective

    1993:    Input Parameters:
    1994: +  bv - the basis vectors context
    1995: .  L  - left BV obtained with BVGetSplit()
    1996: -  R  - right BV obtained with BVGetSplit()

    1998:    Note:
    1999:    The arguments must match the corresponding call to `BVGetSplit()`.

    2001:    Level: advanced

    2003: .seealso: [](sec:bv), `BVGetSplit()`
    2004: @*/
    2005: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
    2006: {
    2007:   PetscFunctionBegin;
    2010:   BVCheckSizes(bv,1);
    2011:   PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
    2012:   PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
    2013:   PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
    2014:   PetscCheck(!L || ((*L)->ci[0]<=(*L)->nc-1 && (*L)->ci[1]<=(*L)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 2 has unrestored columns, use BVRestoreColumn()");
    2015:   PetscCheck(!R || ((*R)->ci[0]<=(*R)->nc-1 && (*R)->ci[1]<=(*R)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 3 has unrestored columns, use BVRestoreColumn()");

    2017:   PetscTryTypeMethod(bv,restoresplit,L,R);
    2018:   bv->lsplit = 0;
    2019:   if (L) *L = NULL;
    2020:   if (R) *R = NULL;
    2021:   PetscFunctionReturn(PETSC_SUCCESS);
    2022: }

    2024: /*
    2025:    Copy all user-provided attributes of V to another BV object W with different layout
    2026:  */
    2027: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
    2028: {
    2029:   PetscFunctionBegin;
    2030:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
    2031:   PetscCall(BVSetVecType(W,V->vtype));
    2032:   W->ld           = V->ld;
    2033:   PetscCall(BVSetType(W,((PetscObject)V)->type_name));
    2034:   W->orthog_type  = V->orthog_type;
    2035:   W->orthog_ref   = V->orthog_ref;
    2036:   W->orthog_eta   = V->orthog_eta;
    2037:   W->orthog_block = V->orthog_block;
    2038:   W->vmm          = V->vmm;
    2039:   W->rrandom      = V->rrandom;
    2040:   W->deftol       = V->deftol;
    2041:   PetscCall(PetscObjectReference((PetscObject)V->rand));
    2042:   W->rand         = V->rand;
    2043:   W->sfocalled    = V->sfocalled;
    2044:   PetscTryTypeMethod(V,duplicate,W);
    2045:   PetscCall(PetscObjectStateIncrease((PetscObject)W));
    2046:   PetscFunctionReturn(PETSC_SUCCESS);
    2047: }

    2049: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
    2050: {
    2051:   PetscInt  rstart,rend,lstart;
    2052:   PetscInt  N,n;
    2053:   PetscBool contig;

    2055:   PetscFunctionBegin;
    2056:   PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
    2057:   PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
    2058:   PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
    2059:   if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
    2060:   else bv->lsplit = -lstart;
    2061:   PetscCall(ISGetSize(is,&N));
    2062:   PetscCall(ISGetLocalSize(is,&n));
    2063:   if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
    2064:   if (!*split) {
    2065:     PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
    2066:     (*split)->issplit = top? -1: -2;
    2067:     (*split)->splitparent = bv;
    2068:     (*split)->N = N;
    2069:     (*split)->n = n;
    2070:     (*split)->m = bv->m;
    2071:     PetscCall(BVDuplicateNewLayout_Private(bv,*split));
    2072:   }
    2073:   (*split)->k  = bv->k;
    2074:   (*split)->l  = bv->l;
    2075:   (*split)->nc = bv->nc;
    2076:   if ((*split)->nc) {
    2077:     (*split)->ci[0] = -(*split)->nc-1;
    2078:     (*split)->ci[1] = -(*split)->nc-1;
    2079:   }
    2080:   if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
    2081:   else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
    2082:   PetscFunctionReturn(PETSC_SUCCESS);
    2083: }

    2085: /*@
    2086:    BVGetSplitRows - Splits the `BV` object into two `BV` objects that share the
    2087:    internal data, using a disjoint horizontal splitting.

    2089:    Collective

    2091:    Input Parameters:
    2092: +  bv   - the basis vectors context
    2093: .  isup - the index set that defines the upper part of the horizontal splitting
    2094: -  islo - the index set that defines the lower part of the horizontal splitting

    2096:    Output Parameters:
    2097: +  U - the resulting `BV` containing the upper rows
    2098: -  L - the resulting `BV` containing the lower rows

    2100:    Notes:
    2101:    The index sets must be such that every MPI process can extract the selected
    2102:    rows from its local part of the input `BV`, and this part must be contiguous.
    2103:    With one process, `isup` will list contiguous indices starting from 0, and `islo`
    2104:    will contain the remaining indices, hence we refer to upper and lower part.
    2105:    However, with several processes the indices will be interleaved because
    2106:    `isup` will refer to the upper part of the local array.

    2108:    The intended use of this function is with matrices of `MATNEST` type, where
    2109:    `MatNestGetISs()` will return the appropriate index sets.

    2111:    The returned `BV`s must be seen as references (not copies) of the input
    2112:    `BV`, that is, modifying them will change the entries of `bv` as well.
    2113:    The returned `BV`s must not be destroyed. `BVRestoreSplitRows()` must be called
    2114:    when they are no longer needed.

    2116:    Pass `NULL` for any of the output `BV`s that is not needed.

    2118:    Level: advanced

    2120: .seealso: [](sec:bv), `BVRestoreSplitRows()`, `BVGetSplit()`
    2121: @*/
    2122: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
    2123: {
    2124:   PetscFunctionBegin;
    2127:   BVCheckSizes(bv,1);
    2130:   PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
    2131:   PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
    2132:   if (U) {
    2133:     PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
    2134:     *U = bv->R;
    2135:   }
    2136:   if (L) {
    2137:     PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
    2138:     *L = bv->L;
    2139:   }
    2140:   PetscFunctionReturn(PETSC_SUCCESS);
    2141: }

    2143: /*@
    2144:    BVRestoreSplitRows - Restore the BV objects obtained with `BVGetSplitRows()`.

    2146:    Logically Collective

    2148:    Input Parameters:
    2149: +  bv - the basis vectors context
    2150: .  isup - the index set that defines the upper part of the horizontal splitting
    2151: .  islo - the index set that defines the lower part of the horizontal splitting
    2152: .  U - the `BV` containing the upper rows
    2153: -  L - the `BV` containing the lower rows

    2155:    Note:
    2156:    The arguments must match the corresponding call to `BVGetSplitRows()`.

    2158:    Level: advanced

    2160: .seealso: [](sec:bv), `BVGetSplitRows()`
    2161: @*/
    2162: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
    2163: {
    2164:   PetscFunctionBegin;
    2167:   BVCheckSizes(bv,1);
    2168:   PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
    2169:   PetscCheck(!U || ((*U)->ci[0]<=(*U)->nc-1 && (*U)->ci[1]<=(*U)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The upper BV has unrestored columns, use BVRestoreColumn()");
    2170:   PetscCheck(!L || ((*L)->ci[0]<=(*L)->nc-1 && (*L)->ci[1]<=(*L)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The lower BV has unrestored columns, use BVRestoreColumn()");
    2171:   PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
    2172:   bv->lsplit = 0;
    2173:   PetscFunctionReturn(PETSC_SUCCESS);
    2174: }

    2176: /*@
    2177:    BVSetDefiniteTolerance - Set the tolerance to be used when checking a
    2178:    definite inner product.

    2180:    Logically Collective

    2182:    Input Parameters:
    2183: +  bv     - basis vectors
    2184: -  deftol - the tolerance

    2186:    Options Database Key:
    2187: .  -bv_definite_tol deftol - the tolerance

    2189:    Notes:
    2190:    When using a non-standard inner product, see `BVSetMatrix()`, the solver needs
    2191:    to compute $\sqrt{z^*B z}$ for various vectors $z$. If the inner product has not
    2192:    been declared indefinite, the value $z^*B z$ must be positive, but due to
    2193:    rounding error a tiny value may become negative. A tolerance is used to
    2194:    detect this situation. Likewise, in complex arithmetic $z^*B z$ should be
    2195:    real, and we use the same tolerance to check whether a nonzero imaginary part
    2196:    can be considered negligible.

    2198:    This function sets this tolerance, which defaults to `10*PETSC_MACHINE_EPSILON`.
    2199:    The default value should be good for most applications.

    2201:    Level: advanced

    2203: .seealso: [](sec:bv), `BVSetMatrix()`
    2204: @*/
    2205: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
    2206: {
    2207:   PetscFunctionBegin;
    2210:   if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
    2211:   else {
    2212:     PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
    2213:     bv->deftol = deftol;
    2214:   }
    2215:   PetscFunctionReturn(PETSC_SUCCESS);
    2216: }

    2218: /*@
    2219:    BVGetDefiniteTolerance - Returns the tolerance for checking a definite
    2220:    inner product.

    2222:    Not Collective

    2224:    Input Parameter:
    2225: .  bv - the basis vectors

    2227:    Output Parameter:
    2228: .  deftol - the tolerance

    2230:    Level: advanced

    2232: .seealso: [](sec:bv), `BVSetDefiniteTolerance()`
    2233: @*/
    2234: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
    2235: {
    2236:   PetscFunctionBegin;
    2238:   PetscAssertPointer(deftol,2);
    2239:   *deftol = bv->deftol;
    2240:   PetscFunctionReturn(PETSC_SUCCESS);
    2241: }

    2243: /*@
    2244:    BVSetLeadingDimension - Set the leading dimension to be used for storing the `BV` data.

    2246:    Not Collective

    2248:    Input Parameters:
    2249: +  bv - basis vectors
    2250: -  ld - the leading dimension

    2252:    Notes:
    2253:    This parameter is relevant for `BVMAT`, though it might be employed in other types
    2254:    as well.

    2256:    When the internal data of the `BV` is stored as a dense matrix, the leading dimension
    2257:    has the same meaning as in `MatDenseSetLDA()`, i.e., the distance in number of
    2258:    elements from one entry of the matrix to the one in the next column at the same
    2259:    row. The leading dimension refers to the local array, and hence can be different
    2260:    in different processes.

    2262:    The user does not need to change this parameter. The default value is equal to the
    2263:    number of local rows, but this value may be increased a little to guarantee alignment
    2264:    (especially in the case of GPU storage).

    2266:    Level: advanced

    2268: .seealso: [](sec:bv), `BVGetLeadingDimension()`
    2269: @*/
    2270: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
    2271: {
    2272:   PetscFunctionBegin;
    2275:   PetscCheck((bv->n<0 && bv->N<0) || !bv->ops->create,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Must call BVSetLeadingDimension() before setting the BV type and sizes");
    2276:   if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
    2277:   else {
    2278:     PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
    2279:     bv->ld = ld;
    2280:   }
    2281:   PetscFunctionReturn(PETSC_SUCCESS);
    2282: }

    2284: /*@
    2285:    BVGetLeadingDimension - Returns the leading dimension of the `BV`.

    2287:    Not Collective

    2289:    Input Parameter:
    2290: .  bv - the basis vectors

    2292:    Output Parameter:
    2293: .  ld - the leading dimension

    2295:    Level: advanced

    2297:    Notes:
    2298:    The returned value may be different in different processes.

    2300:    The leading dimension must be used when accessing the internal array via
    2301:    `BVGetArray()` or `BVGetArrayRead()`.

    2303: .seealso: [](sec:bv), `BVSetLeadingDimension()`, `BVGetArray()`, `BVGetArrayRead()`
    2304: @*/
    2305: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
    2306: {
    2307:   PetscFunctionBegin;
    2309:   PetscAssertPointer(ld,2);
    2310:   *ld = bv->ld;
    2311:   PetscFunctionReturn(PETSC_SUCCESS);
    2312: }