Actual source code: bvbasic.c

slepc-main 2024-11-22
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    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: 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: 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 processor calls this with N of PETSC_DECIDE then all processors must,
104:    otherwise the program will hang.

106:    Level: beginner

108: .seealso: 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:    Level: beginner

164: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
165: @*/
166: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
167: {
168:   PetscInt       ma;
169:   PetscLayout    map;
170:   VecType        vtype;

172:   PetscFunctionBegin;
175:   PetscCheckSameComm(bv,1,t,2);
177:   PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
178:   PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
179:   PetscCall(VecGetType(t,&vtype));
180:   PetscCall(BVSetVecType(bv,vtype));
181:   PetscCall(VecGetLayout(t,&map));
182:   PetscCall(PetscLayoutReference(map,&bv->map));
183:   PetscCall(VecGetSize(t,&bv->N));
184:   PetscCall(VecGetLocalSize(t,&bv->n));
185:   if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
186:     PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
187:     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);
188:   }
189:   bv->m = m;
190:   bv->k = m;
191:   if (bv->ops->create) {
192:     PetscUseTypeMethod(bv,create);
193:     bv->ops->create = NULL;
194:     bv->defersfo = PETSC_FALSE;
195:   }
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

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

202:    Not Collective

204:    Input Parameter:
205: .  bv - the basis vectors

207:    Output Parameters:
208: +  n  - the local size
209: .  N  - the global size
210: -  m  - the number of columns

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

219:    Level: beginner

221: .seealso: BVSetSizes(), BVSetSizesFromVec()
222: @*/
223: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
224: {
225:   PetscFunctionBegin;
226:   if (!bv) {
227:     if (m && !n && !N) *m = 0;
228:     PetscFunctionReturn(PETSC_SUCCESS);
229:   }
231:   if (n || N) BVCheckSizes(bv,1);
232:   if (n) *n = bv->n;
233:   if (N) *N = bv->N;
234:   if (m) *m = bv->m;
235:   if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: /*@
240:    BVSetNumConstraints - Set the number of constraints.

242:    Logically Collective

244:    Input Parameters:
245: +  V  - basis vectors
246: -  nc - number of constraints

248:    Notes:
249:    This function sets the number of constraints to nc and marks all remaining
250:    columns as regular. Normal user would call BVInsertConstraints() instead.

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

256:    Level: developer

258: .seealso: BVInsertConstraints()
259: @*/
260: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
261: {
262:   PetscInt       total,diff,i;
263:   Vec            x,y;

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

273:   diff = nc-V->nc;
274:   if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
275:   total = V->nc+V->m;
276:   PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
277:   if (diff<0) {  /* lessen constraints, shift contents of BV */
278:     for (i=0;i<V->m;i++) {
279:       PetscCall(BVGetColumn(V,i,&x));
280:       PetscCall(BVGetColumn(V,i+diff,&y));
281:       PetscCall(VecCopy(x,y));
282:       PetscCall(BVRestoreColumn(V,i,&x));
283:       PetscCall(BVRestoreColumn(V,i+diff,&y));
284:     }
285:   }
286:   V->nc = nc;
287:   V->ci[0] = -V->nc-1;
288:   V->ci[1] = -V->nc-1;
289:   V->m = total-nc;
290:   V->l = PetscMin(V->l,V->m);
291:   V->k = PetscMin(V->k,V->m);
292:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: /*@
297:    BVGetNumConstraints - Returns the number of constraints.

299:    Not Collective

301:    Input Parameter:
302: .  bv - the basis vectors

304:    Output Parameters:
305: .  nc - the number of constraints

307:    Level: advanced

309: .seealso: BVGetSizes(), BVInsertConstraints()
310: @*/
311: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
312: {
313:   PetscFunctionBegin;
315:   PetscAssertPointer(nc,2);
316:   *nc = bv->nc;
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: /*@
321:    BVResize - Change the number of columns.

323:    Collective

325:    Input Parameters:
326: +  bv   - the basis vectors
327: .  m    - the new number of columns
328: -  copy - a flag indicating whether current values should be kept

330:    Note:
331:    Internal storage is reallocated. If the copy flag is set to true, then
332:    the contents are copied to the leading part of the new space.

334:    Level: advanced

336: .seealso: BVSetSizes(), BVSetSizesFromVec()
337: @*/
338: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
339: {
340:   PetscScalar       *array;
341:   const PetscScalar *omega;
342:   Vec               v;

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

354:   PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
355:   PetscUseTypeMethod(bv,resize,m,copy);
356:   PetscCall(VecDestroy(&bv->buffer));
357:   PetscCall(BVDestroy(&bv->cached));
358:   PetscCall(PetscFree2(bv->h,bv->c));
359:   if (bv->omega) {
360:     if (bv->cuda) {
361: #if defined(PETSC_HAVE_CUDA)
362:       PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
363: #else
364:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
365: #endif
366:     } else if (bv->hip) {
367: #if defined(PETSC_HAVE_HIP)
368:       PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
369: #else
370:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
371: #endif
372:     } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
373:     if (copy) {
374:       PetscCall(VecGetArray(v,&array));
375:       PetscCall(VecGetArrayRead(bv->omega,&omega));
376:       PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
377:       PetscCall(VecRestoreArrayRead(bv->omega,&omega));
378:       PetscCall(VecRestoreArray(v,&array));
379:     } else PetscCall(VecSet(v,1.0));
380:     PetscCall(VecDestroy(&bv->omega));
381:     bv->omega = v;
382:   }
383:   bv->m = m;
384:   bv->k = PetscMin(bv->k,m);
385:   bv->l = PetscMin(bv->l,m);
386:   PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
387:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

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

394:    Logically Collective

396:    Input Parameters:
397: +  bv - the basis vectors context
398: .  l  - number of leading columns
399: -  k  - number of active columns

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

406:    Also in operations such as BVMult() or BVDot(), the first l columns are
407:    normally not included in the computation. See the manpage of each
408:    operation.

410:    In orthogonalization operations, the first l columns are treated
411:    differently, they participate in the orthogonalization but the computed
412:    coefficients are not stored.

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

417:    Level: intermediate

419: .seealso: BVGetActiveColumns(), BVSetSizes()
420: @*/
421: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
422: {
423:   PetscFunctionBegin;
427:   BVCheckSizes(bv,1);
428:   if (PetscUnlikely(k == PETSC_DETERMINE)) {
429:     bv->k = bv->m;
430:   } else if (k != PETSC_CURRENT) {
431:     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);
432:     bv->k = k;
433:   }
434:   if (PetscUnlikely(l == PETSC_DETERMINE)) {
435:     bv->l = 0;
436:   } else if (l != PETSC_CURRENT) {
437:     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);
438:     bv->l = l;
439:   }
440:   PetscFunctionReturn(PETSC_SUCCESS);
441: }

443: /*@
444:    BVGetActiveColumns - Returns the current active dimensions.

446:    Not Collective

448:    Input Parameter:
449: .  bv - the basis vectors context

451:    Output Parameters:
452: +  l  - number of leading columns
453: -  k  - number of active columns

455:    Level: intermediate

457: .seealso: BVSetActiveColumns()
458: @*/
459: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
460: {
461:   PetscFunctionBegin;
463:   if (l) *l = bv->l;
464:   if (k) *k = bv->k;
465:   PetscFunctionReturn(PETSC_SUCCESS);
466: }

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

471:    Collective

473:    Input Parameters:
474: +  bv    - the basis vectors context
475: .  B     - a symmetric matrix (may be NULL)
476: -  indef - a flag indicating if the matrix is indefinite

478:    Notes:
479:    This is used to specify a non-standard inner product, whose matrix
480:    representation is given by B. Then, all inner products required during
481:    orthogonalization are computed as (x,y)_B=y^H*B*x rather than the
482:    standard form (x,y)=y^H*x.

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

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

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

493:    Level: advanced

495: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
496: @*/
497: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
498: {
499:   PetscInt       m,n;

501:   PetscFunctionBegin;
504:   if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
505:     if (B) {
507:       PetscCall(MatGetLocalSize(B,&m,&n));
508:       PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
509:       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);
510:     }
511:     if (B) PetscCall(PetscObjectReference((PetscObject)B));
512:     PetscCall(MatDestroy(&bv->matrix));
513:     bv->matrix = B;
514:     bv->indef  = indef;
515:     PetscCall(PetscObjectStateIncrease((PetscObject)bv));
516:     if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
517:     if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
518:   }
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

522: /*@
523:    BVGetMatrix - Retrieves the matrix representation of the inner product.

525:    Not Collective

527:    Input Parameter:
528: .  bv    - the basis vectors context

530:    Output Parameters:
531: +  B     - the matrix of the inner product (may be NULL)
532: -  indef - the flag indicating if the matrix is indefinite

534:    Level: advanced

536: .seealso: BVSetMatrix()
537: @*/
538: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
539: {
540:   PetscFunctionBegin;
542:   if (B)     *B     = bv->matrix;
543:   if (indef) *indef = bv->indef;
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: /*@
548:    BVApplyMatrix - Multiplies a vector by the matrix representation of the
549:    inner product.

551:    Neighbor-wise Collective

553:    Input Parameters:
554: +  bv - the basis vectors context
555: -  x  - the vector

557:    Output Parameter:
558: .  y  - the result

560:    Note:
561:    If no matrix was specified this function copies the vector.

563:    Level: advanced

565: .seealso: BVSetMatrix(), BVApplyMatrixBV()
566: @*/
567: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
568: {
569:   PetscFunctionBegin;
573:   if (bv->matrix) {
574:     PetscCall(BV_IPMatMult(bv,x));
575:     PetscCall(VecCopy(bv->Bx,y));
576:   } else PetscCall(VecCopy(x,y));
577:   PetscFunctionReturn(PETSC_SUCCESS);
578: }

580: /*@
581:    BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
582:    of the inner product.

584:    Neighbor-wise Collective

586:    Input Parameter:
587: .  X - the basis vectors context

589:    Output Parameter:
590: .  Y - the basis vectors to store the result (optional)

592:    Note:
593:    This function computes Y = B*X, where B is the matrix given with
594:    BVSetMatrix(). This operation is computed as in BVMatMult().
595:    If no matrix was specified, then it just copies Y = X.

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

599:    Level: developer

601: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
602: @*/
603: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
604: {
605:   PetscFunctionBegin;
607:   if (Y) {
609:     if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
610:     else PetscCall(BVCopy(X,Y));
611:   } else PetscCall(BV_IPMatMultBV(X));
612:   PetscFunctionReturn(PETSC_SUCCESS);
613: }

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

618:    Logically Collective

620:    Input Parameters:
621: +  bv    - the basis vectors context
622: -  omega - a vector representing the diagonal of the signature matrix

624:    Note:
625:    The signature matrix Omega = V'*B*V is relevant only for an indefinite B.

627:    Level: developer

629: .seealso: BVSetMatrix(), BVGetSignature()
630: @*/
631: PetscErrorCode BVSetSignature(BV bv,Vec omega)
632: {
633:   PetscInt          i,n;
634:   const PetscScalar *pomega;
635:   PetscScalar       *intern;

637:   PetscFunctionBegin;
639:   BVCheckSizes(bv,1);

643:   PetscCall(VecGetSize(omega,&n));
644:   PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
645:   PetscCall(BV_AllocateSignature(bv));
646:   if (bv->indef) {
647:     PetscCall(VecGetArrayRead(omega,&pomega));
648:     PetscCall(VecGetArray(bv->omega,&intern));
649:     for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
650:     PetscCall(VecRestoreArray(bv->omega,&intern));
651:     PetscCall(VecRestoreArrayRead(omega,&pomega));
652:   } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
653:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
654:   PetscFunctionReturn(PETSC_SUCCESS);
655: }

657: /*@
658:    BVGetSignature - Retrieves the signature matrix from last orthogonalization.

660:    Not Collective

662:    Input Parameter:
663: .  bv    - the basis vectors context

665:    Output Parameter:
666: .  omega - a vector representing the diagonal of the signature matrix

668:    Note:
669:    The signature matrix Omega = V'*B*V is relevant only for an indefinite B.

671:    Level: developer

673: .seealso: BVSetMatrix(), BVSetSignature()
674: @*/
675: PetscErrorCode BVGetSignature(BV bv,Vec omega)
676: {
677:   PetscInt          i,n;
678:   PetscScalar       *pomega;
679:   const PetscScalar *intern;

681:   PetscFunctionBegin;
683:   BVCheckSizes(bv,1);

687:   PetscCall(VecGetSize(omega,&n));
688:   PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
689:   if (bv->indef && bv->omega) {
690:     PetscCall(VecGetArray(omega,&pomega));
691:     PetscCall(VecGetArrayRead(bv->omega,&intern));
692:     for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
693:     PetscCall(VecRestoreArrayRead(bv->omega,&intern));
694:     PetscCall(VecRestoreArray(omega,&pomega));
695:   } else PetscCall(VecSet(omega,1.0));
696:   PetscFunctionReturn(PETSC_SUCCESS);
697: }

699: /*@
700:    BVSetBufferVec - Attach a vector object to be used as buffer space for
701:    several operations.

703:    Collective

705:    Input Parameters:
706: +  bv     - the basis vectors context)
707: -  buffer - the vector

709:    Notes:
710:    Use BVGetBufferVec() to retrieve the vector (for example, to free it
711:    at the end of the computations).

713:    The vector must be sequential of length (nc+m)*m, where m is the number
714:    of columns of bv and nc is the number of constraints.

716:    Level: developer

718: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
719: @*/
720: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
721: {
722:   PetscInt       ld,n;
723:   PetscMPIInt    size;

725:   PetscFunctionBegin;
728:   BVCheckSizes(bv,1);
729:   PetscCall(VecGetSize(buffer,&n));
730:   ld = bv->m+bv->nc;
731:   PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
732:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
733:   PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");

735:   PetscCall(PetscObjectReference((PetscObject)buffer));
736:   PetscCall(VecDestroy(&bv->buffer));
737:   bv->buffer = buffer;
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }

741: /*@
742:    BVGetBufferVec - Obtain the buffer vector associated with the BV object.

744:    Collective

746:    Input Parameters:
747: .  bv - the basis vectors context

749:    Output Parameter:
750: .  buffer - vector

752:    Notes:
753:    The vector is created if not available previously. It is a sequential vector
754:    of length (nc+m)*m, where m is the number of columns of bv and nc is the number
755:    of constraints.

757:    Developer Notes:
758:    The buffer vector is viewed as a column-major matrix with leading dimension
759:    ld=nc+m, and m columns at most. In the most common usage, it has the structure
760: .vb
761:       | | C |
762:       |s|---|
763:       | | H |
764: .ve
765:    where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
766:    related to orthogonalization against constraints (first nc rows), and s is the
767:    first column that contains scratch values computed during Gram-Schmidt
768:    orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
769:    store the coefficients.

771:    Level: developer

773: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
774: @*/
775: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
776: {
777:   PetscInt       ld;

779:   PetscFunctionBegin;
781:   PetscAssertPointer(buffer,2);
782:   BVCheckSizes(bv,1);
783:   if (!bv->buffer) {
784:     ld = bv->m+bv->nc;
785:     PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
786:     PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
787:     PetscCall(VecSetType(bv->buffer,bv->vtype));
788:   }
789:   *buffer = bv->buffer;
790:   PetscFunctionReturn(PETSC_SUCCESS);
791: }

793: /*@
794:    BVSetRandomContext - Sets the PetscRandom object associated with the BV,
795:    to be used in operations that need random numbers.

797:    Collective

799:    Input Parameters:
800: +  bv   - the basis vectors context
801: -  rand - the random number generator context

803:    Level: advanced

805: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
806: @*/
807: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
808: {
809:   PetscFunctionBegin;
812:   PetscCheckSameComm(bv,1,rand,2);
813:   PetscCall(PetscObjectReference((PetscObject)rand));
814:   PetscCall(PetscRandomDestroy(&bv->rand));
815:   bv->rand = rand;
816:   PetscFunctionReturn(PETSC_SUCCESS);
817: }

819: /*@
820:    BVGetRandomContext - Gets the PetscRandom object associated with the BV.

822:    Collective

824:    Input Parameter:
825: .  bv - the basis vectors context

827:    Output Parameter:
828: .  rand - the random number generator context

830:    Level: advanced

832: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
833: @*/
834: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
835: {
836:   PetscFunctionBegin;
838:   PetscAssertPointer(rand,2);
839:   if (!bv->rand) {
840:     PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
841:     if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
842:     if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
843:     if (bv->rrandom) {
844:       PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
845:       PetscCall(PetscRandomSeed(bv->rand));
846:     }
847:   }
848:   *rand = bv->rand;
849:   PetscFunctionReturn(PETSC_SUCCESS);
850: }

852: /*@
853:    BVSetFromOptions - Sets BV options from the options database.

855:    Collective

857:    Input Parameter:
858: .  bv - the basis vectors context

860:    Level: beginner

862: .seealso: BVSetOptionsPrefix()
863: @*/
864: PetscErrorCode BVSetFromOptions(BV bv)
865: {
866:   char               type[256];
867:   PetscBool          flg1,flg2,flg3,flg4;
868:   PetscReal          r;
869:   BVOrthogType       otype;
870:   BVOrthogRefineType orefine;
871:   BVOrthogBlockType  oblock;

873:   PetscFunctionBegin;
875:   PetscCall(BVRegisterAll());
876:   PetscObjectOptionsBegin((PetscObject)bv);
877:     PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
878:     if (flg1) PetscCall(BVSetType(bv,type));
879:     else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));

881:     otype = bv->orthog_type;
882:     PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
883:     orefine = bv->orthog_ref;
884:     PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
885:     r = bv->orthog_eta;
886:     PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
887:     oblock = bv->orthog_block;
888:     PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
889:     if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));

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

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

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

899:     if (bv->ops->create) bv->defersfo = PETSC_TRUE;   /* defer call to setfromoptions */
900:     else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
901:     PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
902:   PetscOptionsEnd();
903:   bv->sfocalled = PETSC_TRUE;
904:   PetscFunctionReturn(PETSC_SUCCESS);
905: }

907: /*@
908:    BVSetOrthogonalization - Specifies the method used for the orthogonalization of
909:    vectors (classical or modified Gram-Schmidt with or without refinement), and
910:    for the block-orthogonalization (simultaneous orthogonalization of a set of
911:    vectors).

913:    Logically Collective

915:    Input Parameters:
916: +  bv     - the basis vectors context
917: .  type   - the method of vector orthogonalization
918: .  refine - type of refinement
919: .  eta    - parameter for selective refinement
920: -  block  - the method of block orthogonalization

922:    Options Database Keys:
923: +  -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
924:                          (default) or mgs for Modified Gram-Schmidt orthogonalization
925: .  -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
926: .  -bv_orthog_eta <eta> -  For setting the value of eta
927: -  -bv_orthog_block <block> - Where <block> is the block-orthogonalization method

929:    Notes:
930:    The default settings work well for most problems.

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

936:    When using several processors, MGS is likely to result in bad scalability.

938:    If the method set for block orthogonalization is GS, then the computation
939:    is done column by column with the vector orthogonalization.

941:    Level: advanced

943: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
944: @*/
945: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
946: {
947:   PetscFunctionBegin;
953:   switch (type) {
954:     case BV_ORTHOG_CGS:
955:     case BV_ORTHOG_MGS:
956:       bv->orthog_type = type;
957:       break;
958:     default:
959:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
960:   }
961:   switch (refine) {
962:     case BV_ORTHOG_REFINE_NEVER:
963:     case BV_ORTHOG_REFINE_IFNEEDED:
964:     case BV_ORTHOG_REFINE_ALWAYS:
965:       bv->orthog_ref = refine;
966:       break;
967:     default:
968:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
969:   }
970:   if (eta == (PetscReal)PETSC_DETERMINE) {
971:     bv->orthog_eta = 0.7071;
972:   } else if (eta != (PetscReal)PETSC_CURRENT) {
973:     PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
974:     bv->orthog_eta = eta;
975:   }
976:   switch (block) {
977:     case BV_ORTHOG_BLOCK_GS:
978:     case BV_ORTHOG_BLOCK_CHOL:
979:     case BV_ORTHOG_BLOCK_TSQR:
980:     case BV_ORTHOG_BLOCK_TSQRCHOL:
981:     case BV_ORTHOG_BLOCK_SVQB:
982:       bv->orthog_block = block;
983:       break;
984:     default:
985:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
986:   }
987:   PetscFunctionReturn(PETSC_SUCCESS);
988: }

990: /*@
991:    BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.

993:    Not Collective

995:    Input Parameter:
996: .  bv - basis vectors context

998:    Output Parameters:
999: +  type   - the method of vector orthogonalization
1000: .  refine - type of refinement
1001: .  eta    - parameter for selective refinement
1002: -  block  - the method of block orthogonalization

1004:    Level: advanced

1006: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
1007: @*/
1008: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1009: {
1010:   PetscFunctionBegin;
1012:   if (type)   *type   = bv->orthog_type;
1013:   if (refine) *refine = bv->orthog_ref;
1014:   if (eta)    *eta    = bv->orthog_eta;
1015:   if (block)  *block  = bv->orthog_block;
1016:   PetscFunctionReturn(PETSC_SUCCESS);
1017: }

1019: /*@
1020:    BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.

1022:    Logically Collective

1024:    Input Parameters:
1025: +  bv     - the basis vectors context
1026: -  method - the method for the BVMatMult() operation

1028:    Options Database Keys:
1029: .  -bv_matmult <meth> - choose one of the methods: vecs, mat

1031:    Notes:
1032:    Allowed values are
1033: +  BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1034: .  BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1035: -  BV_MATMULT_MAT_SAVE - this case is deprecated

1037:    The default is BV_MATMULT_MAT except in the case of BVVECS.

1039:    Level: advanced

1041: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1042: @*/
1043: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1044: {
1045:   PetscFunctionBegin;
1048:   switch (method) {
1049:     case BV_MATMULT_VECS:
1050:     case BV_MATMULT_MAT:
1051:       bv->vmm = method;
1052:       break;
1053:     case BV_MATMULT_MAT_SAVE:
1054:       PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1055:       bv->vmm = BV_MATMULT_MAT;
1056:       break;
1057:     default:
1058:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1059:   }
1060:   PetscFunctionReturn(PETSC_SUCCESS);
1061: }

1063: /*@
1064:    BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.

1066:    Not Collective

1068:    Input Parameter:
1069: .  bv - basis vectors context

1071:    Output Parameter:
1072: .  method - the method for the BVMatMult() operation

1074:    Level: advanced

1076: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1077: @*/
1078: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1079: {
1080:   PetscFunctionBegin;
1082:   PetscAssertPointer(method,2);
1083:   *method = bv->vmm;
1084:   PetscFunctionReturn(PETSC_SUCCESS);
1085: }

1087: /*@
1088:    BVGetColumn - Returns a Vec object that contains the entries of the
1089:    requested column of the basis vectors object.

1091:    Logically Collective

1093:    Input Parameters:
1094: +  bv - the basis vectors context
1095: -  j  - the index of the requested column

1097:    Output Parameter:
1098: .  v  - vector containing the jth column

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

1104:    The returned Vec must not be destroyed. BVRestoreColumn() must be
1105:    called when it is no longer needed. At most, two columns can be fetched,
1106:    that is, this function can only be called twice before the corresponding
1107:    BVRestoreColumn() is invoked.

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

1112:    Level: beginner

1114: .seealso: BVRestoreColumn(), BVInsertConstraints()
1115: @*/
1116: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1117: {
1118:   PetscInt       l;

1120:   PetscFunctionBegin;
1123:   BVCheckSizes(bv,1);
1124:   BVCheckOp(bv,1,getcolumn);
1126:   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);
1127:   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);
1128:   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);
1129:   l = BVAvailableVec;
1130:   PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1131:   PetscUseTypeMethod(bv,getcolumn,j,v);
1132:   bv->ci[l] = j;
1133:   PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
1134:   PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1135:   *v = bv->cv[l];
1136:   PetscFunctionReturn(PETSC_SUCCESS);
1137: }

1139: /*@
1140:    BVRestoreColumn - Restore a column obtained with BVGetColumn().

1142:    Logically Collective

1144:    Input Parameters:
1145: +  bv - the basis vectors context
1146: .  j  - the index of the column
1147: -  v  - vector obtained with BVGetColumn()

1149:    Note:
1150:    The arguments must match the corresponding call to BVGetColumn().

1152:    Level: beginner

1154: .seealso: BVGetColumn()
1155: @*/
1156: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1157: {
1158:   PetscObjectId    id;
1159:   PetscObjectState st;
1160:   PetscInt         l;

1162:   PetscFunctionBegin;
1165:   BVCheckSizes(bv,1);
1167:   PetscAssertPointer(v,3);
1169:   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);
1170:   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);
1171:   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);
1172:   l = (j==bv->ci[0])? 0: 1;
1173:   PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1174:   PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1175:   PetscCall(VecGetState(*v,&st));
1176:   if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1177:   PetscUseTypeMethod(bv,restorecolumn,j,v);
1178:   bv->ci[l] = -bv->nc-1;
1179:   bv->st[l] = -1;
1180:   bv->id[l] = 0;
1181:   *v = NULL;
1182:   PetscFunctionReturn(PETSC_SUCCESS);
1183: }

1185: /*@C
1186:    BVGetArray - Returns a pointer to a contiguous array that contains this
1187:    processor's portion of the BV data.

1189:    Logically Collective

1191:    Input Parameters:
1192: .  bv - the basis vectors context

1194:    Output Parameter:
1195: .  a  - location to put pointer to the array

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

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

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

1209:    Use BVGetArrayRead() for read-only access.

1211:    Level: advanced

1213: .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArrayRead()
1214: @*/
1215: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1216: {
1217:   PetscFunctionBegin;
1220:   BVCheckSizes(bv,1);
1221:   BVCheckOp(bv,1,getarray);
1222:   PetscUseTypeMethod(bv,getarray,a);
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: /*@C
1227:    BVRestoreArray - Restore the BV object after BVGetArray() has been called.

1229:    Logically Collective

1231:    Input Parameters:
1232: +  bv - the basis vectors context
1233: -  a  - location of pointer to array obtained from BVGetArray()

1235:    Note:
1236:    This operation may imply a data copy, for BV types that do not store
1237:    data contiguously in memory.

1239:    Level: advanced

1241: .seealso: BVGetColumn()
1242: @*/
1243: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1244: {
1245:   PetscFunctionBegin;
1248:   BVCheckSizes(bv,1);
1249:   PetscTryTypeMethod(bv,restorearray,a);
1250:   if (a) *a = NULL;
1251:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1252:   PetscFunctionReturn(PETSC_SUCCESS);
1253: }

1255: /*@C
1256:    BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1257:    contains this processor's portion of the BV data.

1259:    Not Collective

1261:    Input Parameters:
1262: .  bv - the basis vectors context

1264:    Output Parameter:
1265: .  a  - location to put pointer to the array

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

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

1275:    Level: advanced

1277: .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArray()
1278: @*/
1279: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1280: {
1281:   PetscFunctionBegin;
1284:   BVCheckSizes(bv,1);
1285:   BVCheckOp(bv,1,getarrayread);
1286:   PetscUseTypeMethod(bv,getarrayread,a);
1287:   PetscFunctionReturn(PETSC_SUCCESS);
1288: }

1290: /*@C
1291:    BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1292:    been called.

1294:    Not Collective

1296:    Input Parameters:
1297: +  bv - the basis vectors context
1298: -  a  - location of pointer to array obtained from BVGetArrayRead()

1300:    Level: advanced

1302: .seealso: BVGetColumn()
1303: @*/
1304: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1305: {
1306:   PetscFunctionBegin;
1309:   BVCheckSizes(bv,1);
1310:   PetscTryTypeMethod(bv,restorearrayread,a);
1311:   if (a) *a = NULL;
1312:   PetscFunctionReturn(PETSC_SUCCESS);
1313: }

1315: /*@
1316:    BVCreateVec - Creates a new Vec object with the same type and dimensions
1317:    as the columns of the basis vectors object.

1319:    Collective

1321:    Input Parameter:
1322: .  bv - the basis vectors context

1324:    Output Parameter:
1325: .  v  - the new vector

1327:    Note:
1328:    The user is responsible of destroying the returned vector.

1330:    Level: beginner

1332: .seealso: BVCreateMat(), BVCreateVecEmpty()
1333: @*/
1334: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1335: {
1336:   PetscFunctionBegin;
1338:   BVCheckSizes(bv,1);
1339:   PetscAssertPointer(v,2);
1340:   PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1341:   PetscCall(VecSetLayout(*v,bv->map));
1342:   PetscCall(VecSetType(*v,bv->vtype));
1343:   PetscCall(VecSetUp(*v));
1344:   PetscFunctionReturn(PETSC_SUCCESS);
1345: }

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

1351:    Collective

1353:    Input Parameter:
1354: .  bv - the basis vectors context

1356:    Output Parameter:
1357: .  v  - the new vector

1359:    Note:
1360:    This works as BVCreateVec(), but the new vector does not have the array allocated,
1361:    so the intended usage is with VecPlaceArray().

1363:    Level: developer

1365: .seealso: BVCreateVec()
1366: @*/
1367: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1368: {
1369:   PetscBool  standard,cuda,hip,mpi;
1370:   PetscInt   N,nloc,bs;

1372:   PetscFunctionBegin;
1374:   BVCheckSizes(bv,1);
1375:   PetscAssertPointer(v,2);

1377:   PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1378:   PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1379:   PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
1380:   if (standard || cuda || hip) {
1381:     PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
1382:     PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1383:     PetscCall(PetscLayoutGetSize(bv->map,&N));
1384:     PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1385:     if (cuda) {
1386: #if defined(PETSC_HAVE_CUDA)
1387:       if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1388:       else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1389: #endif
1390:     } else if (hip) {
1391: #if defined(PETSC_HAVE_HIP)
1392:       if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1393:       else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1394: #endif
1395:     } else {
1396:       if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1397:       else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1398:     }
1399:   } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1400:   PetscFunctionReturn(PETSC_SUCCESS);
1401: }

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

1406:    Collective

1408:    Input Parameters:
1409: +  bv - the basis vectors context
1410: -  vtype - the vector type

1412:    Level: advanced

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

1418: .seealso: BVGetVecType(), BVSetSizesFromVec(), BVSetSizes()
1419: @*/
1420: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1421: {
1422:   PetscBool   std;
1423:   PetscMPIInt size;

1425:   PetscFunctionBegin;
1427:   PetscCall(PetscFree(bv->vtype));
1428:   PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1429:   if (std) {
1430:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1431:     PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1432:   } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1433:   PetscFunctionReturn(PETSC_SUCCESS);
1434: }

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

1439:    Not Collective

1441:    Input Parameter:
1442: .  bv - the basis vectors context

1444:    Output Parameter:
1445: .  vtype - the vector type

1447:    Level: advanced

1449: .seealso: BVSetVecType()
1450: @*/
1451: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1452: {
1453:   PetscFunctionBegin;
1455:   PetscAssertPointer(vtype,2);
1456:   *vtype = bv->vtype;
1457:   PetscFunctionReturn(PETSC_SUCCESS);
1458: }

1460: /*@
1461:    BVCreateMat - Creates a new Mat object of dense type and copies the contents
1462:    of the BV object.

1464:    Collective

1466:    Input Parameter:
1467: .  bv - the basis vectors context

1469:    Output Parameter:
1470: .  A  - the new matrix

1472:    Notes:
1473:    The user is responsible of destroying the returned matrix.

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

1477:    Level: intermediate

1479: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1480: @*/
1481: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1482: {
1483:   PetscInt ksave,lsave;
1484:   Mat      B;

1486:   PetscFunctionBegin;
1488:   BVCheckSizes(bv,1);
1489:   PetscAssertPointer(A,2);

1491:   PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1492:   lsave = bv->l;
1493:   ksave = bv->k;
1494:   bv->l = 0;
1495:   bv->k = bv->m;
1496:   PetscCall(BVGetMat(bv,&B));
1497:   PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1498:   PetscCall(BVRestoreMat(bv,&B));
1499:   bv->l = lsave;
1500:   bv->k = ksave;
1501:   PetscFunctionReturn(PETSC_SUCCESS);
1502: }

1504: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1505: {
1506:   PetscScalar *vv,*aa;
1507:   PetscBool   create=PETSC_FALSE;
1508:   PetscInt    m,cols;

1510:   PetscFunctionBegin;
1511:   m = bv->k-bv->l;
1512:   if (!bv->Aget) create=PETSC_TRUE;
1513:   else {
1514:     PetscCall(MatDenseGetArray(bv->Aget,&aa));
1515:     PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1516:     PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1517:     if (cols!=m) {
1518:       PetscCall(MatDestroy(&bv->Aget));
1519:       create=PETSC_TRUE;
1520:     }
1521:   }
1522:   PetscCall(BVGetArray(bv,&vv));
1523:   if (create) {
1524:     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 */
1525:     PetscCall(MatDenseReplaceArray(bv->Aget,NULL));  /* replace with a null pointer, the value after BVRestoreMat */
1526:   }
1527:   PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld));  /* set the actual pointer */
1528:   *A = bv->Aget;
1529:   PetscFunctionReturn(PETSC_SUCCESS);
1530: }

1532: /*@
1533:    BVGetMat - Returns a Mat object of dense type that shares the memory of
1534:    the BV object.

1536:    Collective

1538:    Input Parameter:
1539: .  bv - the basis vectors context

1541:    Output Parameter:
1542: .  A  - the matrix

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

1549:    This operation implies a call to BVGetArray(), which may result in data
1550:    copies.

1552:    Level: advanced

1554: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1555: @*/
1556: PetscErrorCode BVGetMat(BV bv,Mat *A)
1557: {
1558:   char name[64];

1560:   PetscFunctionBegin;
1562:   BVCheckSizes(bv,1);
1563:   PetscAssertPointer(A,2);
1564:   PetscUseTypeMethod(bv,getmat,A);
1565:   if (((PetscObject)bv)->name) {   /* set A's name based on BV name */
1566:     PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1567:     PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1568:     PetscCall(PetscObjectSetName((PetscObject)*A,name));
1569:   }
1570:   PetscFunctionReturn(PETSC_SUCCESS);
1571: }

1573: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1574: {
1575:   PetscScalar *vv,*aa;

1577:   PetscFunctionBegin;
1578:   PetscCall(MatDenseGetArray(bv->Aget,&aa));
1579:   vv = aa-(bv->nc+bv->l)*bv->ld;
1580:   PetscCall(MatDenseResetArray(bv->Aget));
1581:   PetscCall(BVRestoreArray(bv,&vv));
1582:   *A = NULL;
1583:   PetscFunctionReturn(PETSC_SUCCESS);
1584: }

1586: /*@
1587:    BVRestoreMat - Restores the Mat obtained with BVGetMat().

1589:    Logically Collective

1591:    Input Parameters:
1592: +  bv - the basis vectors context
1593: -  A  - the fetched matrix

1595:    Note:
1596:    A call to this function must match a previous call of BVGetMat().
1597:    The effect is that the contents of the Mat are copied back to the
1598:    BV internal data structures.

1600:    Level: advanced

1602: .seealso: BVGetMat(), BVRestoreArray()
1603: @*/
1604: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1605: {
1606:   PetscFunctionBegin;
1608:   BVCheckSizes(bv,1);
1609:   PetscAssertPointer(A,2);
1610:   PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1611:   PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1612:   PetscUseTypeMethod(bv,restoremat,A);
1613:   PetscFunctionReturn(PETSC_SUCCESS);
1614: }

1616: /*
1617:    Copy all user-provided attributes of V to another BV object W
1618:  */
1619: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1620: {
1621:   PetscFunctionBegin;
1622:   PetscCall(PetscLayoutReference(V->map,&W->map));
1623:   PetscCall(BVSetVecType(W,V->vtype));
1624:   W->ld           = V->ld;
1625:   PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1626:   W->orthog_type  = V->orthog_type;
1627:   W->orthog_ref   = V->orthog_ref;
1628:   W->orthog_eta   = V->orthog_eta;
1629:   W->orthog_block = V->orthog_block;
1630:   if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1631:   W->matrix       = V->matrix;
1632:   W->indef        = V->indef;
1633:   W->vmm          = V->vmm;
1634:   W->rrandom      = V->rrandom;
1635:   W->deftol       = V->deftol;
1636:   if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1637:   W->rand         = V->rand;
1638:   W->sfocalled    = V->sfocalled;
1639:   PetscTryTypeMethod(V,duplicate,W);
1640:   PetscCall(PetscObjectStateIncrease((PetscObject)W));
1641:   PetscFunctionReturn(PETSC_SUCCESS);
1642: }

1644: /*@
1645:    BVDuplicate - Creates a new basis vector object of the same type and
1646:    dimensions as an existing one.

1648:    Collective

1650:    Input Parameter:
1651: .  V - basis vectors context

1653:    Output Parameter:
1654: .  W - location to put the new BV

1656:    Notes:
1657:    The new BV has the same type and dimensions as V. Also, the inner
1658:    product matrix and orthogonalization options are copied.

1660:    BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1661:    for the new basis vectors. Use BVCopy() to copy the contents.

1663:    Level: intermediate

1665: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1666: @*/
1667: PetscErrorCode BVDuplicate(BV V,BV *W)
1668: {
1669:   PetscFunctionBegin;
1672:   BVCheckSizes(V,1);
1673:   PetscAssertPointer(W,2);
1674:   PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1675:   (*W)->N = V->N;
1676:   (*W)->n = V->n;
1677:   (*W)->m = V->m;
1678:   (*W)->k = V->m;
1679:   PetscCall(BVDuplicate_Private(V,*W));
1680:   PetscFunctionReturn(PETSC_SUCCESS);
1681: }

1683: /*@
1684:    BVDuplicateResize - Creates a new basis vector object of the same type and
1685:    dimensions as an existing one, but with possibly different number of columns.

1687:    Collective

1689:    Input Parameters:
1690: +  V - basis vectors context
1691: -  m - the new number of columns

1693:    Output Parameter:
1694: .  W - location to put the new BV

1696:    Note:
1697:    This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1698:    contents of V are not copied to W.

1700:    Level: intermediate

1702: .seealso: BVDuplicate(), BVResize()
1703: @*/
1704: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1705: {
1706:   PetscFunctionBegin;
1709:   BVCheckSizes(V,1);
1711:   PetscAssertPointer(W,3);
1712:   PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1713:   (*W)->N = V->N;
1714:   (*W)->n = V->n;
1715:   (*W)->m = m;
1716:   (*W)->k = m;
1717:   PetscCall(BVDuplicate_Private(V,*W));
1718:   PetscFunctionReturn(PETSC_SUCCESS);
1719: }

1721: /*@
1722:    BVGetCachedBV - Returns a BV object stored internally that holds the
1723:    result of B*X after a call to BVApplyMatrixBV().

1725:    Collective

1727:    Input Parameter:
1728: .  bv    - the basis vectors context

1730:    Output Parameter:
1731: .  cached - the cached BV

1733:    Note:
1734:    The cached BV is created if not available previously.

1736:    Level: developer

1738: .seealso: BVApplyMatrixBV()
1739: @*/
1740: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1741: {
1742:   PetscFunctionBegin;
1744:   PetscAssertPointer(cached,2);
1745:   BVCheckSizes(bv,1);
1746:   if (!bv->cached) {
1747:     PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1748:     bv->cached->N = bv->N;
1749:     bv->cached->n = bv->n;
1750:     bv->cached->m = bv->m;
1751:     bv->cached->k = bv->m;
1752:     PetscCall(BVDuplicate_Private(bv,bv->cached));
1753:   }
1754:   *cached = bv->cached;
1755:   PetscFunctionReturn(PETSC_SUCCESS);
1756: }

1758: /*@
1759:    BVCopy - Copies a basis vector object into another one, W <- V.

1761:    Logically Collective

1763:    Input Parameter:
1764: .  V - basis vectors context

1766:    Output Parameter:
1767: .  W - the copy

1769:    Note:
1770:    Both V and W must be distributed in the same manner; local copies are
1771:    done. Only active columns (excluding the leading ones) are copied.
1772:    In the destination W, columns are overwritten starting from the leading ones.
1773:    Constraints are not copied.

1775:    Level: beginner

1777: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1778: @*/
1779: PetscErrorCode BVCopy(BV V,BV W)
1780: {
1781:   PetscScalar       *womega;
1782:   const PetscScalar *vomega;

1784:   PetscFunctionBegin;
1787:   BVCheckSizes(V,1);
1788:   BVCheckOp(V,1,copy);
1791:   BVCheckSizes(W,2);
1792:   PetscCheckSameTypeAndComm(V,1,W,2);
1793:   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);
1794:   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);
1795:   if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);

1797:   PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1798:   if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1799:     /* copy signature */
1800:     PetscCall(BV_AllocateSignature(W));
1801:     PetscCall(VecGetArrayRead(V->omega,&vomega));
1802:     PetscCall(VecGetArray(W->omega,&womega));
1803:     PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1804:     PetscCall(VecRestoreArray(W->omega,&womega));
1805:     PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1806:   }
1807:   PetscUseTypeMethod(V,copy,W);
1808:   PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1809:   PetscCall(PetscObjectStateIncrease((PetscObject)W));
1810:   PetscFunctionReturn(PETSC_SUCCESS);
1811: }

1813: /*@
1814:    BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.

1816:    Logically Collective

1818:    Input Parameters:
1819: +  V - basis vectors context
1820: -  j - the column number to be copied

1822:    Output Parameter:
1823: .  w - the copied column

1825:    Note:
1826:    Both V and w must be distributed in the same manner; local copies are done.

1828:    Level: beginner

1830: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1831: @*/
1832: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1833: {
1834:   PetscInt       n,N;
1835:   Vec            z;

1837:   PetscFunctionBegin;
1840:   BVCheckSizes(V,1);
1843:   PetscCheckSameComm(V,1,w,3);

1845:   PetscCall(VecGetSize(w,&N));
1846:   PetscCall(VecGetLocalSize(w,&n));
1847:   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);

1849:   PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1850:   PetscCall(BVGetColumn(V,j,&z));
1851:   PetscCall(VecCopy(z,w));
1852:   PetscCall(BVRestoreColumn(V,j,&z));
1853:   PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1854:   PetscFunctionReturn(PETSC_SUCCESS);
1855: }

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

1860:    Logically Collective

1862:    Input Parameters:
1863: +  V - basis vectors context
1864: .  j - the number of the source column
1865: -  i - the number of the destination column

1867:    Level: beginner

1869: .seealso: BVCopy(), BVCopyVec()
1870: @*/
1871: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1872: {
1873:   PetscScalar *omega;

1875:   PetscFunctionBegin;
1878:   BVCheckSizes(V,1);
1881:   if (j==i) PetscFunctionReturn(PETSC_SUCCESS);

1883:   PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1884:   if (V->omega) {
1885:     PetscCall(VecGetArray(V->omega,&omega));
1886:     omega[i] = omega[j];
1887:     PetscCall(VecRestoreArray(V->omega,&omega));
1888:   }
1889:   PetscUseTypeMethod(V,copycolumn,j,i);
1890:   PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1891:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
1892:   PetscFunctionReturn(PETSC_SUCCESS);
1893: }

1895: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1896: {
1897:   PetscInt  ncols;

1899:   PetscFunctionBegin;
1900:   ncols = left? bv->nc+bv->l: bv->m-bv->l;
1901:   if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1902:   if (!*split) {
1903:     PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1904:     (*split)->issplit = left? 1: 2;
1905:     (*split)->splitparent = bv;
1906:     (*split)->N = bv->N;
1907:     (*split)->n = bv->n;
1908:     (*split)->m = bv->m;
1909:     (*split)->k = bv->m;
1910:     PetscCall(BVDuplicate_Private(bv,*split));
1911:   }
1912:   (*split)->l  = 0;
1913:   (*split)->k  = left? bv->l: bv->k-bv->l;
1914:   (*split)->nc = left? bv->nc: 0;
1915:   (*split)->m  = ncols-(*split)->nc;
1916:   if ((*split)->nc) {
1917:     (*split)->ci[0] = -(*split)->nc-1;
1918:     (*split)->ci[1] = -(*split)->nc-1;
1919:   }
1920:   if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1921:   else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1922:   PetscFunctionReturn(PETSC_SUCCESS);
1923: }

1925: /*@
1926:    BVGetSplit - Splits the BV object into two BV objects that share the
1927:    internal data, one of them containing the leading columns and the other
1928:    one containing the remaining columns.

1930:    Collective

1932:    Input Parameter:
1933: .  bv - the basis vectors context

1935:    Output Parameters:
1936: +  L - left BV containing leading columns (can be NULL)
1937: -  R - right BV containing remaining columns (can be NULL)

1939:    Notes:
1940:    The columns are split in two sets. The leading columns (including the
1941:    constraints) are assigned to the left BV and the remaining columns
1942:    are assigned to the right BV. The number of leading columns, as
1943:    specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1944:    guarantee that both L and R have at least one column).

1946:    The returned BV's must be seen as references (not copies) of the input
1947:    BV, that is, modifying them will change the entries of bv as well.
1948:    The returned BV's must not be destroyed. BVRestoreSplit() must be called
1949:    when they are no longer needed.

1951:    Pass NULL for any of the output BV's that is not needed.

1953:    Level: advanced

1955: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints(), BVGetSplitRows()
1956: @*/
1957: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1958: {
1959:   PetscFunctionBegin;
1962:   BVCheckSizes(bv,1);
1963:   PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1964:   PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1965:   PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1966:   bv->lsplit = bv->nc+bv->l;
1967:   PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1968:   PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1969:   if (L) *L = bv->L;
1970:   if (R) *R = bv->R;
1971:   PetscFunctionReturn(PETSC_SUCCESS);
1972: }

1974: /*@
1975:    BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().

1977:    Logically Collective

1979:    Input Parameters:
1980: +  bv - the basis vectors context
1981: .  L  - left BV obtained with BVGetSplit()
1982: -  R  - right BV obtained with BVGetSplit()

1984:    Note:
1985:    The arguments must match the corresponding call to BVGetSplit().

1987:    Level: advanced

1989: .seealso: BVGetSplit()
1990: @*/
1991: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1992: {
1993:   PetscFunctionBegin;
1996:   BVCheckSizes(bv,1);
1997:   PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1998:   PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1999:   PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
2000:   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()");
2001:   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()");

2003:   PetscTryTypeMethod(bv,restoresplit,L,R);
2004:   bv->lsplit = 0;
2005:   if (L) *L = NULL;
2006:   if (R) *R = NULL;
2007:   PetscFunctionReturn(PETSC_SUCCESS);
2008: }

2010: /*
2011:    Copy all user-provided attributes of V to another BV object W with different layout
2012:  */
2013: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
2014: {
2015:   PetscFunctionBegin;
2016:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2017:   PetscCall(BVSetVecType(W,V->vtype));
2018:   W->ld           = V->ld;
2019:   PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2020:   W->orthog_type  = V->orthog_type;
2021:   W->orthog_ref   = V->orthog_ref;
2022:   W->orthog_eta   = V->orthog_eta;
2023:   W->orthog_block = V->orthog_block;
2024:   W->vmm          = V->vmm;
2025:   W->rrandom      = V->rrandom;
2026:   W->deftol       = V->deftol;
2027:   if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
2028:   W->rand         = V->rand;
2029:   W->sfocalled    = V->sfocalled;
2030:   PetscTryTypeMethod(V,duplicate,W);
2031:   PetscCall(PetscObjectStateIncrease((PetscObject)W));
2032:   PetscFunctionReturn(PETSC_SUCCESS);
2033: }

2035: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2036: {
2037:   PetscInt  rstart,rend,lstart;
2038:   PetscInt  N,n;
2039:   PetscBool contig;

2041:   PetscFunctionBegin;
2042:   PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2043:   PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2044:   PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2045:   if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2046:   else bv->lsplit = -lstart;
2047:   PetscCall(ISGetSize(is,&N));
2048:   PetscCall(ISGetLocalSize(is,&n));
2049:   if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2050:   if (!*split) {
2051:     PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2052:     (*split)->issplit = top? -1: -2;
2053:     (*split)->splitparent = bv;
2054:     (*split)->N = N;
2055:     (*split)->n = n;
2056:     (*split)->m = bv->m;
2057:     PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2058:   }
2059:   (*split)->k  = bv->k;
2060:   (*split)->l  = bv->l;
2061:   (*split)->nc = bv->nc;
2062:   if ((*split)->nc) {
2063:     (*split)->ci[0] = -(*split)->nc-1;
2064:     (*split)->ci[1] = -(*split)->nc-1;
2065:   }
2066:   if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2067:   else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2068:   PetscFunctionReturn(PETSC_SUCCESS);
2069: }

2071: /*@
2072:    BVGetSplitRows - Splits the BV object into two BV objects that share the
2073:    internal data, using a disjoint horizontal splitting.

2075:    Collective

2077:    Input Parameters:
2078: +  bv   - the basis vectors context
2079: .  isup - the index set that defines the upper part of the horizontal splitting
2080: -  islo - the index set that defines the lower part of the horizontal splitting

2082:    Output Parameters:
2083: +  U - the resulting BV containing the upper rows
2084: -  L - the resulting BV containing the lower rows

2086:    Notes:
2087:    The index sets must be such that every MPI process can extract the selected
2088:    rows from its local part of the input BV, and this part must be contiguous.
2089:    With one process, isup will list contiguous indices starting from 0, and islo
2090:    will contain the remaining indices, hence we refer to upper and lower part.
2091:    However, with several processes the indices will be interleaved because
2092:    isup will refer to the upper part of the local array.

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

2097:    The returned BV's must be seen as references (not copies) of the input
2098:    BV, that is, modifying them will change the entries of bv as well.
2099:    The returned BV's must not be destroyed. BVRestoreSplitRows() must be called
2100:    when they are no longer needed.

2102:    Pass NULL for any of the output BV's that is not needed.

2104:    Level: advanced

2106: .seealso: BVRestoreSplitRows(), BVGetSplit()
2107: @*/
2108: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2109: {
2110:   PetscFunctionBegin;
2113:   BVCheckSizes(bv,1);
2116:   PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2117:   PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2118:   if (U) {
2119:     PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2120:     *U = bv->R;
2121:   }
2122:   if (L) {
2123:     PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2124:     *L = bv->L;
2125:   }
2126:   PetscFunctionReturn(PETSC_SUCCESS);
2127: }

2129: /*@
2130:    BVRestoreSplitRows - Restore the BV objects obtained with BVGetSplitRows().

2132:    Logically Collective

2134:    Input Parameters:
2135: +  bv - the basis vectors context
2136: .  isup - the index set that defines the upper part of the horizontal splitting
2137: .  islo - the index set that defines the lower part of the horizontal splitting
2138: .  U - the BV containing the upper rows
2139: -  L - the BV containing the lower rows

2141:    Note:
2142:    The arguments must match the corresponding call to BVGetSplitRows().

2144:    Level: advanced

2146: .seealso: BVGetSplitRows()
2147: @*/
2148: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2149: {
2150:   PetscFunctionBegin;
2153:   BVCheckSizes(bv,1);
2154:   PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2155:   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()");
2156:   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()");
2157:   PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2158:   bv->lsplit = 0;
2159:   PetscFunctionReturn(PETSC_SUCCESS);
2160: }

2162: /*@
2163:    BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2164:    definite inner product.

2166:    Logically Collective

2168:    Input Parameters:
2169: +  bv     - basis vectors
2170: -  deftol - the tolerance

2172:    Options Database Key:
2173: .  -bv_definite_tol <deftol> - the tolerance

2175:    Notes:
2176:    When using a non-standard inner product, see BVSetMatrix(), the solver needs
2177:    to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
2178:    been declared indefinite, the value z'*B*z must be positive, but due to
2179:    rounding error a tiny value may become negative. A tolerance is used to
2180:    detect this situation. Likewise, in complex arithmetic z'*B*z should be
2181:    real, and we use the same tolerance to check whether a nonzero imaginary part
2182:    can be considered negligible.

2184:    This function sets this tolerance, which defaults to 10*eps, where eps is
2185:    the machine epsilon. The default value should be good for most applications.

2187:    Level: advanced

2189: .seealso: BVSetMatrix()
2190: @*/
2191: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2192: {
2193:   PetscFunctionBegin;
2196:   if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2197:   else {
2198:     PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2199:     bv->deftol = deftol;
2200:   }
2201:   PetscFunctionReturn(PETSC_SUCCESS);
2202: }

2204: /*@
2205:    BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2206:    inner product.

2208:    Not Collective

2210:    Input Parameter:
2211: .  bv - the basis vectors

2213:    Output Parameter:
2214: .  deftol - the tolerance

2216:    Level: advanced

2218: .seealso: BVSetDefiniteTolerance()
2219: @*/
2220: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2221: {
2222:   PetscFunctionBegin;
2224:   PetscAssertPointer(deftol,2);
2225:   *deftol = bv->deftol;
2226:   PetscFunctionReturn(PETSC_SUCCESS);
2227: }

2229: /*@
2230:    BVSetLeadingDimension - Set the leading dimension to be used for storing the BV data.

2232:    Not Collective

2234:    Input Parameters:
2235: +  bv - basis vectors
2236: -  ld - the leading dimension

2238:    Notes:
2239:    This parameter is relevant for BVMAT, though it might be employed in other types
2240:    as well.

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

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

2252:    Level: advanced

2254: .seealso: BVGetLeadingDimension()
2255: @*/
2256: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2257: {
2258:   PetscFunctionBegin;
2261:   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");
2262:   if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
2263:   else {
2264:     PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2265:     bv->ld = ld;
2266:   }
2267:   PetscFunctionReturn(PETSC_SUCCESS);
2268: }

2270: /*@
2271:    BVGetLeadingDimension - Returns the leading dimension of the BV.

2273:    Not Collective

2275:    Input Parameter:
2276: .  bv - the basis vectors

2278:    Output Parameter:
2279: .  ld - the leading dimension

2281:    Level: advanced

2283:    Notes:
2284:    The returned value may be different in different processes.

2286:    The leading dimension must be used when accessing the internal array via
2287:    BVGetArray() or BVGetArrayRead().

2289: .seealso: BVSetLeadingDimension(), BVGetArray(), BVGetArrayRead()
2290: @*/
2291: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2292: {
2293:   PetscFunctionBegin;
2295:   PetscAssertPointer(ld,2);
2296:   *ld = bv->ld;
2297:   PetscFunctionReturn(PETSC_SUCCESS);
2298: }