Actual source code: bvbasic.c

slepc-3.17.2 2022-08-09
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 = 0;

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

 22:    Logically Collective on bv

 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;


 43:   PetscObjectTypeCompare((PetscObject)bv,type,&match);
 44:   if (match) PetscFunctionReturn(0);
 45:   PetscStrcmp(type,BVTENSOR,&match);

 48:   PetscFunctionListFind(BVList,type,&r);

 51:   if (bv->ops->destroy) (*bv->ops->destroy)(bv);
 52:   PetscMemzero(bv->ops,sizeof(struct _BVOps));

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

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

 68:    Not Collective

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

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

 76:    Level: intermediate

 78: .seealso: BVSetType()
 79: @*/
 80: PetscErrorCode BVGetType(BV bv,BVType *type)
 81: {
 84:   *type = ((PetscObject)bv)->type_name;
 85:   PetscFunctionReturn(0);
 86: }

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

 91:    Collective on bv

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

 99:    Notes:
100:    n and N cannot be both PETSC_DECIDE.
101:    If one processor calls this with N of PETSC_DECIDE then all processors must,
102:    otherwise the program will hang.

104:    Level: beginner

106: .seealso: BVSetSizesFromVec(), BVGetSizes(), BVResize()
107: @*/
108: PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)
109: {
110:   PetscInt       ma;

119:   bv->n = n;
120:   bv->N = N;
121:   bv->m = m;
122:   bv->k = m;
123:   if (!bv->t) {  /* create template vector and get actual dimensions */
124:     VecCreate(PetscObjectComm((PetscObject)bv),&bv->t);
125:     VecSetSizes(bv->t,bv->n,bv->N);
126:     VecSetFromOptions(bv->t);
127:     VecGetSize(bv->t,&bv->N);
128:     VecGetLocalSize(bv->t,&bv->n);
129:     if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
130:       MatGetLocalSize(bv->matrix,&ma,NULL);
132:     }
133:   }
134:   if (bv->ops->create) {
135:     PetscLogEventBegin(BV_Create,bv,0,0,0);
136:     (*bv->ops->create)(bv);
137:     PetscLogEventEnd(BV_Create,bv,0,0,0);
138:     bv->ops->create = 0;
139:     bv->defersfo = PETSC_FALSE;
140:   }
141:   PetscFunctionReturn(0);
142: }

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

148:    Collective on bv

150:    Input Parameters:
151: +  bv - the basis vectors
152: .  t  - the template vector
153: -  m  - the number of columns

155:    Level: beginner

157: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
158: @*/
159: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
160: {
161:   PetscInt       ma;

169:   VecGetSize(t,&bv->N);
170:   VecGetLocalSize(t,&bv->n);
171:   if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
172:     MatGetLocalSize(bv->matrix,&ma,NULL);
174:   }
175:   bv->m = m;
176:   bv->k = m;
177:   bv->t = t;
178:   PetscObjectReference((PetscObject)t);
179:   if (bv->ops->create) {
180:     (*bv->ops->create)(bv);
181:     bv->ops->create = 0;
182:     bv->defersfo = PETSC_FALSE;
183:   }
184:   PetscFunctionReturn(0);
185: }

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

190:    Not Collective

192:    Input Parameter:
193: .  bv - the basis vectors

195:    Output Parameters:
196: +  n  - the local size
197: .  N  - the global size
198: -  m  - the number of columns

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

207:    Level: beginner

209: .seealso: BVSetSizes(), BVSetSizesFromVec()
210: @*/
211: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
212: {
213:   if (!bv) {
214:     if (m && !n && !N) *m = 0;
215:     PetscFunctionReturn(0);
216:   }
218:   if (n || N) BVCheckSizes(bv,1);
219:   if (n) *n = bv->n;
220:   if (N) *N = bv->N;
221:   if (m) *m = bv->m;
222:   if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
223:   PetscFunctionReturn(0);
224: }

226: /*@
227:    BVSetNumConstraints - Set the number of constraints.

229:    Logically Collective on V

231:    Input Parameters:
232: +  V  - basis vectors
233: -  nc - number of constraints

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

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

243:    Level: developer

245: .seealso: BVInsertConstraints()
246: @*/
247: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
248: {
249:   PetscInt       total,diff,i;
250:   Vec            x,y;

256:   BVCheckSizes(V,1);

259:   diff = nc-V->nc;
260:   if (!diff) PetscFunctionReturn(0);
261:   total = V->nc+V->m;
263:   if (diff<0) {  /* lessen constraints, shift contents of BV */
264:     for (i=0;i<V->m;i++) {
265:       BVGetColumn(V,i,&x);
266:       BVGetColumn(V,i+diff,&y);
267:       VecCopy(x,y);
268:       BVRestoreColumn(V,i,&x);
269:       BVRestoreColumn(V,i+diff,&y);
270:     }
271:   }
272:   V->nc = nc;
273:   V->ci[0] = -V->nc-1;
274:   V->ci[1] = -V->nc-1;
275:   V->m = total-nc;
276:   V->l = PetscMin(V->l,V->m);
277:   V->k = PetscMin(V->k,V->m);
278:   PetscObjectStateIncrease((PetscObject)V);
279:   PetscFunctionReturn(0);
280: }

282: /*@
283:    BVGetNumConstraints - Returns the number of constraints.

285:    Not Collective

287:    Input Parameter:
288: .  bv - the basis vectors

290:    Output Parameters:
291: .  nc - the number of constraints

293:    Level: advanced

295: .seealso: BVGetSizes(), BVInsertConstraints()
296: @*/
297: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
298: {
301:   *nc = bv->nc;
302:   PetscFunctionReturn(0);
303: }

305: /*@
306:    BVResize - Change the number of columns.

308:    Collective on bv

310:    Input Parameters:
311: +  bv   - the basis vectors
312: .  m    - the new number of columns
313: -  copy - a flag indicating whether current values should be kept

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

319:    Level: advanced

321: .seealso: BVSetSizes(), BVSetSizesFromVec()
322: @*/
323: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
324: {
325:   PetscScalar       *array;
326:   const PetscScalar *omega;
327:   Vec               v;

335:   if (bv->m == m) PetscFunctionReturn(0);
336:   BVCheckOp(bv,1,resize);

338:   PetscLogEventBegin(BV_Create,bv,0,0,0);
339:   (*bv->ops->resize)(bv,m,copy);
340:   VecDestroy(&bv->buffer);
341:   BVDestroy(&bv->cached);
342:   PetscFree2(bv->h,bv->c);
343:   if (bv->omega) {
344:     if (bv->cuda) {
345: #if defined(PETSC_HAVE_CUDA)
346:       VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v);
347: #else
348:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
349: #endif
350:     } else VecCreateSeq(PETSC_COMM_SELF,m,&v);
351:     PetscLogObjectParent((PetscObject)bv,(PetscObject)v);
352:     if (copy) {
353:       VecGetArray(v,&array);
354:       VecGetArrayRead(bv->omega,&omega);
355:       PetscArraycpy(array,omega,PetscMin(m,bv->m));
356:       VecRestoreArrayRead(bv->omega,&omega);
357:       VecRestoreArray(v,&array);
358:     } else VecSet(v,1.0);
359:     VecDestroy(&bv->omega);
360:     bv->omega = v;
361:   }
362:   bv->m = m;
363:   bv->k = PetscMin(bv->k,m);
364:   bv->l = PetscMin(bv->l,m);
365:   PetscLogEventEnd(BV_Create,bv,0,0,0);
366:   PetscObjectStateIncrease((PetscObject)bv);
367:   PetscFunctionReturn(0);
368: }

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

373:    Logically Collective on bv

375:    Input Parameters:
376: +  bv - the basis vectors context
377: .  l  - number of leading columns
378: -  k  - number of active columns

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

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

389:    In orthogonalization operations, the first l columns are treated
390:    differently, they participate in the orthogonalization but the computed
391:    coefficients are not stored.

393:    Level: intermediate

395: .seealso: BVGetActiveColumns(), BVSetSizes()
396: @*/
397: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
398: {
402:   BVCheckSizes(bv,1);
403:   if (PetscUnlikely(k==PETSC_DECIDE || k==PETSC_DEFAULT)) {
404:     bv->k = bv->m;
405:   } else {
407:     bv->k = k;
408:   }
409:   if (PetscUnlikely(l==PETSC_DECIDE || l==PETSC_DEFAULT)) {
410:     bv->l = 0;
411:   } else {
413:     bv->l = l;
414:   }
415:   PetscFunctionReturn(0);
416: }

418: /*@
419:    BVGetActiveColumns - Returns the current active dimensions.

421:    Not Collective

423:    Input Parameter:
424: .  bv - the basis vectors context

426:    Output Parameters:
427: +  l  - number of leading columns
428: -  k  - number of active columns

430:    Level: intermediate

432: .seealso: BVSetActiveColumns()
433: @*/
434: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
435: {
437:   if (l) *l = bv->l;
438:   if (k) *k = bv->k;
439:   PetscFunctionReturn(0);
440: }

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

445:    Collective on bv

447:    Input Parameters:
448: +  bv    - the basis vectors context
449: .  B     - a symmetric matrix (may be NULL)
450: -  indef - a flag indicating if the matrix is indefinite

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

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

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

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

467:    Level: advanced

469: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
470: @*/
471: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
472: {
473:   PetscInt       m,n;

477:   if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
478:     if (B) {
480:       MatGetLocalSize(B,&m,&n);
483:     }
484:     if (B) PetscObjectReference((PetscObject)B);
485:     MatDestroy(&bv->matrix);
486:     bv->matrix = B;
487:     bv->indef  = indef;
488:     PetscObjectStateIncrease((PetscObject)bv);
489:     if (bv->Bx) PetscObjectStateIncrease((PetscObject)bv->Bx);
490:     if (bv->cached) PetscObjectStateIncrease((PetscObject)bv->cached);
491:   }
492:   PetscFunctionReturn(0);
493: }

495: /*@
496:    BVGetMatrix - Retrieves the matrix representation of the inner product.

498:    Not collective, though a parallel Mat may be returned

500:    Input Parameter:
501: .  bv    - the basis vectors context

503:    Output Parameters:
504: +  B     - the matrix of the inner product (may be NULL)
505: -  indef - the flag indicating if the matrix is indefinite

507:    Level: advanced

509: .seealso: BVSetMatrix()
510: @*/
511: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
512: {
514:   if (B)     *B     = bv->matrix;
515:   if (indef) *indef = bv->indef;
516:   PetscFunctionReturn(0);
517: }

519: /*@
520:    BVApplyMatrix - Multiplies a vector by the matrix representation of the
521:    inner product.

523:    Neighbor-wise Collective on bv

525:    Input Parameters:
526: +  bv - the basis vectors context
527: -  x  - the vector

529:    Output Parameter:
530: .  y  - the result

532:    Note:
533:    If no matrix was specified this function copies the vector.

535:    Level: advanced

537: .seealso: BVSetMatrix(), BVApplyMatrixBV()
538: @*/
539: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
540: {
544:   if (bv->matrix) {
545:     BV_IPMatMult(bv,x);
546:     VecCopy(bv->Bx,y);
547:   } else VecCopy(x,y);
548:   PetscFunctionReturn(0);
549: }

551: /*@
552:    BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
553:    of the inner product.

555:    Neighbor-wise Collective on X

557:    Input Parameter:
558: .  X - the basis vectors context

560:    Output Parameter:
561: .  Y - the basis vectors to store the result (optional)

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

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

570:    Level: developer

572: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
573: @*/
574: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
575: {
577:   if (Y) {
579:     if (X->matrix) BVMatMult(X,X->matrix,Y);
580:     else BVCopy(X,Y);
581:   } else BV_IPMatMultBV(X);
582:   PetscFunctionReturn(0);
583: }

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

588:    Logically Collective on bv

590:    Input Parameters:
591: +  bv    - the basis vectors context
592: -  omega - a vector representing the diagonal of the signature matrix

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

597:    Level: developer

599: .seealso: BVSetMatrix(), BVGetSignature()
600: @*/
601: PetscErrorCode BVSetSignature(BV bv,Vec omega)
602: {
603:   PetscInt          i,n;
604:   const PetscScalar *pomega;
605:   PetscScalar       *intern;

608:   BVCheckSizes(bv,1);

612:   VecGetSize(omega,&n);
614:   BV_AllocateSignature(bv);
615:   if (bv->indef) {
616:     VecGetArrayRead(omega,&pomega);
617:     VecGetArray(bv->omega,&intern);
618:     for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
619:     VecRestoreArray(bv->omega,&intern);
620:     VecRestoreArrayRead(omega,&pomega);
621:   } else PetscInfo(bv,"Ignoring signature because BV is not indefinite\n");
622:   PetscObjectStateIncrease((PetscObject)bv);
623:   PetscFunctionReturn(0);
624: }

626: /*@
627:    BVGetSignature - Retrieves the signature matrix from last orthogonalization.

629:    Not collective

631:    Input Parameter:
632: .  bv    - the basis vectors context

634:    Output Parameter:
635: .  omega - a vector representing the diagonal of the signature matrix

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

640:    Level: developer

642: .seealso: BVSetMatrix(), BVSetSignature()
643: @*/
644: PetscErrorCode BVGetSignature(BV bv,Vec omega)
645: {
646:   PetscInt          i,n;
647:   PetscScalar       *pomega;
648:   const PetscScalar *intern;

651:   BVCheckSizes(bv,1);

655:   VecGetSize(omega,&n);
657:   if (bv->indef && bv->omega) {
658:     VecGetArray(omega,&pomega);
659:     VecGetArrayRead(bv->omega,&intern);
660:     for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
661:     VecRestoreArrayRead(bv->omega,&intern);
662:     VecRestoreArray(omega,&pomega);
663:   } else VecSet(omega,1.0);
664:   PetscFunctionReturn(0);
665: }

667: /*@
668:    BVSetBufferVec - Attach a vector object to be used as buffer space for
669:    several operations.

671:    Collective on bv

673:    Input Parameters:
674: +  bv     - the basis vectors context)
675: -  buffer - the vector

677:    Notes:
678:    Use BVGetBufferVec() to retrieve the vector (for example, to free it
679:    at the end of the computations).

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

684:    Level: developer

686: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
687: @*/
688: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
689: {
690:   PetscInt       ld,n;
691:   PetscMPIInt    size;

695:   BVCheckSizes(bv,1);
696:   VecGetSize(buffer,&n);
697:   ld = bv->m+bv->nc;
699:   MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size);

702:   PetscObjectReference((PetscObject)buffer);
703:   VecDestroy(&bv->buffer);
704:   bv->buffer = buffer;
705:   PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
706:   PetscFunctionReturn(0);
707: }

709: /*@
710:    BVGetBufferVec - Obtain the buffer vector associated with the BV object.

712:    Not Collective, but Vec returned is parallel if BV is parallel

714:    Input Parameters:
715: .  bv - the basis vectors context

717:    Output Parameter:
718: .  buffer - vector

720:    Notes:
721:    The vector is created if not available previously. It is a sequential vector
722:    of length (nc+m)*m, where m is the number of columns of bv and nc is the number
723:    of constraints.

725:    Developer Notes:
726:    The buffer vector is viewed as a column-major matrix with leading dimension
727:    ld=nc+m, and m columns at most. In the most common usage, it has the structure
728: .vb
729:       | | C |
730:       |s|---|
731:       | | H |
732: .ve
733:    where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
734:    related to orthogonalization against constraints (first nc rows), and s is the
735:    first column that contains scratch values computed during Gram-Schmidt
736:    orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
737:    store the coefficients.

739:    Level: developer

741: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
742: @*/
743: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
744: {
745:   PetscInt       ld;

749:   BVCheckSizes(bv,1);
750:   if (!bv->buffer) {
751:     ld = bv->m+bv->nc;
752:     VecCreate(PETSC_COMM_SELF,&bv->buffer);
753:     VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m);
754:     VecSetType(bv->buffer,((PetscObject)bv->t)->type_name);
755:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
756:   }
757:   *buffer = bv->buffer;
758:   PetscFunctionReturn(0);
759: }

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

765:    Collective on bv

767:    Input Parameters:
768: +  bv   - the basis vectors context
769: -  rand - the random number generator context

771:    Level: advanced

773: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
774: @*/
775: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
776: {
780:   PetscObjectReference((PetscObject)rand);
781:   PetscRandomDestroy(&bv->rand);
782:   bv->rand = rand;
783:   PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
784:   PetscFunctionReturn(0);
785: }

787: /*@
788:    BVGetRandomContext - Gets the PetscRandom object associated with the BV.

790:    Not Collective

792:    Input Parameter:
793: .  bv - the basis vectors context

795:    Output Parameter:
796: .  rand - the random number generator context

798:    Level: advanced

800: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
801: @*/
802: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
803: {
806:   if (!bv->rand) {
807:     PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand);
808:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
809:     if (bv->cuda) PetscRandomSetType(bv->rand,PETSCCURAND);
810:     if (bv->sfocalled) PetscRandomSetFromOptions(bv->rand);
811:     if (bv->rrandom) {
812:       PetscRandomSetSeed(bv->rand,0x12345678);
813:       PetscRandomSeed(bv->rand);
814:     }
815:   }
816:   *rand = bv->rand;
817:   PetscFunctionReturn(0);
818: }

820: /*@
821:    BVSetFromOptions - Sets BV options from the options database.

823:    Collective on bv

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

828:    Level: beginner

830: .seealso: BVSetOptionsPrefix()
831: @*/
832: PetscErrorCode BVSetFromOptions(BV bv)
833: {
834:   PetscErrorCode     ierr;
835:   char               type[256];
836:   PetscBool          flg1,flg2,flg3,flg4;
837:   PetscReal          r;
838:   BVOrthogType       otype;
839:   BVOrthogRefineType orefine;
840:   BVOrthogBlockType  oblock;

843:   BVRegisterAll();
844:   ierr = PetscObjectOptionsBegin((PetscObject)bv);
845:     PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVSVEC),type,sizeof(type),&flg1);
846:     if (flg1) BVSetType(bv,type);
847:     else if (!((PetscObject)bv)->type_name) BVSetType(bv,BVSVEC);

849:     otype = bv->orthog_type;
850:     PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1);
851:     orefine = bv->orthog_ref;
852:     PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2);
853:     r = bv->orthog_eta;
854:     PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3);
855:     oblock = bv->orthog_block;
856:     PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4);
857:     if (flg1 || flg2 || flg3 || flg4) BVSetOrthogonalization(bv,otype,orefine,r,oblock);

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

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

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

867:     if (bv->ops->create) bv->defersfo = PETSC_TRUE;   /* defer call to setfromoptions */
868:     else if (bv->ops->setfromoptions) (*bv->ops->setfromoptions)(PetscOptionsObject,bv);
869:     PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)bv);
870:   ierr = PetscOptionsEnd();
871:   bv->sfocalled = PETSC_TRUE;
872:   PetscFunctionReturn(0);
873: }

875: /*@
876:    BVSetOrthogonalization - Specifies the method used for the orthogonalization of
877:    vectors (classical or modified Gram-Schmidt with or without refinement), and
878:    for the block-orthogonalization (simultaneous orthogonalization of a set of
879:    vectors).

881:    Logically Collective on bv

883:    Input Parameters:
884: +  bv     - the basis vectors context
885: .  type   - the method of vector orthogonalization
886: .  refine - type of refinement
887: .  eta    - parameter for selective refinement
888: -  block  - the method of block orthogonalization

890:    Options Database Keys:
891: +  -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
892:                          (default) or mgs for Modified Gram-Schmidt orthogonalization
893: .  -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
894: .  -bv_orthog_eta <eta> -  For setting the value of eta
895: -  -bv_orthog_block <block> - Where <block> is the block-orthogonalization method

897:    Notes:
898:    The default settings work well for most problems.

900:    The parameter eta should be a real value between 0 and 1 (or PETSC_DEFAULT).
901:    The value of eta is used only when the refinement type is "ifneeded".

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

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

908:    Level: advanced

910: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
911: @*/
912: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
913: {
919:   switch (type) {
920:     case BV_ORTHOG_CGS:
921:     case BV_ORTHOG_MGS:
922:       bv->orthog_type = type;
923:       break;
924:     default:
925:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
926:   }
927:   switch (refine) {
928:     case BV_ORTHOG_REFINE_NEVER:
929:     case BV_ORTHOG_REFINE_IFNEEDED:
930:     case BV_ORTHOG_REFINE_ALWAYS:
931:       bv->orthog_ref = refine;
932:       break;
933:     default:
934:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
935:   }
936:   if (eta == PETSC_DEFAULT) {
937:     bv->orthog_eta = 0.7071;
938:   } else {
940:     bv->orthog_eta = eta;
941:   }
942:   switch (block) {
943:     case BV_ORTHOG_BLOCK_GS:
944:     case BV_ORTHOG_BLOCK_CHOL:
945:     case BV_ORTHOG_BLOCK_TSQR:
946:     case BV_ORTHOG_BLOCK_TSQRCHOL:
947:     case BV_ORTHOG_BLOCK_SVQB:
948:       bv->orthog_block = block;
949:       break;
950:     default:
951:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
952:   }
953:   PetscFunctionReturn(0);
954: }

956: /*@
957:    BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.

959:    Not Collective

961:    Input Parameter:
962: .  bv - basis vectors context

964:    Output Parameters:
965: +  type   - the method of vector orthogonalization
966: .  refine - type of refinement
967: .  eta    - parameter for selective refinement
968: -  block  - the method of block orthogonalization

970:    Level: advanced

972: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
973: @*/
974: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
975: {
977:   if (type)   *type   = bv->orthog_type;
978:   if (refine) *refine = bv->orthog_ref;
979:   if (eta)    *eta    = bv->orthog_eta;
980:   if (block)  *block  = bv->orthog_block;
981:   PetscFunctionReturn(0);
982: }

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

987:    Logically Collective on bv

989:    Input Parameters:
990: +  bv     - the basis vectors context
991: -  method - the method for the BVMatMult() operation

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

996:    Notes:
997:    Allowed values are
998: +  BV_MATMULT_VECS - perform a matrix-vector multiply per each column
999: .  BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1000: -  BV_MATMULT_MAT_SAVE - this case is deprecated

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

1004:    Level: advanced

1006: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1007: @*/
1008: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1009: {
1012:   switch (method) {
1013:     case BV_MATMULT_VECS:
1014:     case BV_MATMULT_MAT:
1015:       bv->vmm = method;
1016:       break;
1017:     case BV_MATMULT_MAT_SAVE:
1018:       PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n");
1019:       bv->vmm = BV_MATMULT_MAT;
1020:       break;
1021:     default:
1022:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1023:   }
1024:   PetscFunctionReturn(0);
1025: }

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

1030:    Not Collective

1032:    Input Parameter:
1033: .  bv - basis vectors context

1035:    Output Parameter:
1036: .  method - the method for the BVMatMult() operation

1038:    Level: advanced

1040: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1041: @*/
1042: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1043: {
1046:   *method = bv->vmm;
1047:   PetscFunctionReturn(0);
1048: }

1050: /*@
1051:    BVGetColumn - Returns a Vec object that contains the entries of the
1052:    requested column of the basis vectors object.

1054:    Logically Collective on bv

1056:    Input Parameters:
1057: +  bv - the basis vectors context
1058: -  j  - the index of the requested column

1060:    Output Parameter:
1061: .  v  - vector containing the jth column

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

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

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

1075:    Level: beginner

1077: .seealso: BVRestoreColumn(), BVInsertConstraints()
1078: @*/
1079: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1080: {
1081:   PetscInt       l;

1085:   BVCheckSizes(bv,1);
1086:   BVCheckOp(bv,1,getcolumn);
1091:   l = BVAvailableVec;
1093:   (*bv->ops->getcolumn)(bv,j,v);
1094:   bv->ci[l] = j;
1095:   PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]);
1096:   PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]);
1097:   *v = bv->cv[l];
1098:   PetscFunctionReturn(0);
1099: }

1101: /*@
1102:    BVRestoreColumn - Restore a column obtained with BVGetColumn().

1104:    Logically Collective on bv

1106:    Input Parameters:
1107: +  bv - the basis vectors context
1108: .  j  - the index of the column
1109: -  v  - vector obtained with BVGetColumn()

1111:    Note:
1112:    The arguments must match the corresponding call to BVGetColumn().

1114:    Level: beginner

1116: .seealso: BVGetColumn()
1117: @*/
1118: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1119: {
1120:   PetscObjectId    id;
1121:   PetscObjectState st;
1122:   PetscInt         l;

1126:   BVCheckSizes(bv,1);
1133:   l = (j==bv->ci[0])? 0: 1;
1134:   PetscObjectGetId((PetscObject)*v,&id);
1136:   PetscObjectStateGet((PetscObject)*v,&st);
1137:   if (st!=bv->st[l]) PetscObjectStateIncrease((PetscObject)bv);
1138:   if (bv->ops->restorecolumn) (*bv->ops->restorecolumn)(bv,j,v);
1139:   else bv->cv[l] = NULL;
1140:   bv->ci[l] = -bv->nc-1;
1141:   bv->st[l] = -1;
1142:   bv->id[l] = 0;
1143:   *v = NULL;
1144:   PetscFunctionReturn(0);
1145: }

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

1151:    Logically Collective on bv

1153:    Input Parameters:
1154: .  bv - the basis vectors context

1156:    Output Parameter:
1157: .  a  - location to put pointer to the array

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

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

1167:    Level: advanced

1169: .seealso: BVRestoreArray(), BVInsertConstraints()
1170: @*/
1171: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1172: {
1175:   BVCheckSizes(bv,1);
1176:   BVCheckOp(bv,1,getarray);
1177:   (*bv->ops->getarray)(bv,a);
1178:   PetscFunctionReturn(0);
1179: }

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

1184:    Logically Collective on bv

1186:    Input Parameters:
1187: +  bv - the basis vectors context
1188: -  a  - location of pointer to array obtained from BVGetArray()

1190:    Note:
1191:    This operation may imply a data copy, for BV types that do not store
1192:    data contiguously in memory.

1194:    Level: advanced

1196: .seealso: BVGetColumn()
1197: @*/
1198: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1199: {
1202:   BVCheckSizes(bv,1);
1203:   if (bv->ops->restorearray) (*bv->ops->restorearray)(bv,a);
1204:   if (a) *a = NULL;
1205:   PetscObjectStateIncrease((PetscObject)bv);
1206:   PetscFunctionReturn(0);
1207: }

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

1213:    Not Collective

1215:    Input Parameters:
1216: .  bv - the basis vectors context

1218:    Output Parameter:
1219: .  a  - location to put pointer to the array

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

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

1229:    Level: advanced

1231: .seealso: BVRestoreArray(), BVInsertConstraints()
1232: @*/
1233: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1234: {
1237:   BVCheckSizes(bv,1);
1238:   BVCheckOp(bv,1,getarrayread);
1239:   (*bv->ops->getarrayread)(bv,a);
1240:   PetscFunctionReturn(0);
1241: }

1243: /*@C
1244:    BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1245:    been called.

1247:    Not Collective

1249:    Input Parameters:
1250: +  bv - the basis vectors context
1251: -  a  - location of pointer to array obtained from BVGetArrayRead()

1253:    Level: advanced

1255: .seealso: BVGetColumn()
1256: @*/
1257: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1258: {
1261:   BVCheckSizes(bv,1);
1262:   if (bv->ops->restorearrayread) (*bv->ops->restorearrayread)(bv,a);
1263:   if (a) *a = NULL;
1264:   PetscFunctionReturn(0);
1265: }

1267: /*@
1268:    BVCreateVec - Creates a new Vec object with the same type and dimensions
1269:    as the columns of the basis vectors object.

1271:    Collective on bv

1273:    Input Parameter:
1274: .  bv - the basis vectors context

1276:    Output Parameter:
1277: .  v  - the new vector

1279:    Note:
1280:    The user is responsible of destroying the returned vector.

1282:    Level: beginner

1284: .seealso: BVCreateMat()
1285: @*/
1286: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1287: {
1289:   BVCheckSizes(bv,1);
1291:   VecDuplicate(bv->t,v);
1292:   PetscFunctionReturn(0);
1293: }

1295: /*@
1296:    BVCreateMat - Creates a new Mat object of dense type and copies the contents
1297:    of the BV object.

1299:    Collective on bv

1301:    Input Parameter:
1302: .  bv - the basis vectors context

1304:    Output Parameter:
1305: .  A  - the new matrix

1307:    Notes:
1308:    The user is responsible of destroying the returned matrix.

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

1312:    Level: intermediate

1314: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1315: @*/
1316: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1317: {
1318:   PetscScalar       *aa;
1319:   const PetscScalar *vv;

1322:   BVCheckSizes(bv,1);

1325:   MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,bv->N,bv->m,NULL,A);
1326:   MatDenseGetArrayWrite(*A,&aa);
1327:   BVGetArrayRead(bv,&vv);
1328:   PetscArraycpy(aa,vv,bv->m*bv->n);
1329:   BVRestoreArrayRead(bv,&vv);
1330:   MatDenseRestoreArrayWrite(*A,&aa);
1331:   MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1332:   MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1333:   PetscFunctionReturn(0);
1334: }

1336: /*@
1337:    BVGetMat - Returns a Mat object of dense type that shares the memory of
1338:    the BV object.

1340:    Collective on bv

1342:    Input Parameter:
1343: .  bv - the basis vectors context

1345:    Output Parameter:
1346: .  A  - the matrix

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

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

1356:    Level: advanced

1358: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1359: @*/
1360: PetscErrorCode BVGetMat(BV bv,Mat *A)
1361: {
1362:   PetscScalar    *vv,*aa;
1363:   PetscBool      create=PETSC_FALSE;
1364:   PetscInt       m,cols;

1367:   BVCheckSizes(bv,1);
1369:   if (bv->ops->getmat) (*bv->ops->getmat)(bv,A);
1370:   else {
1371:     m = bv->k-bv->l;
1372:     if (!bv->Aget) create=PETSC_TRUE;
1373:     else {
1374:       MatDenseGetArray(bv->Aget,&aa);
1376:       MatGetSize(bv->Aget,NULL,&cols);
1377:       if (cols!=m) {
1378:         MatDestroy(&bv->Aget);
1379:         create=PETSC_TRUE;
1380:       }
1381:     }
1382:     BVGetArray(bv,&vv);
1383:     if (create) {
1384:       MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,vv,&bv->Aget); /* pass a pointer to avoid allocation of storage */
1385:       MatDenseReplaceArray(bv->Aget,NULL);  /* replace with a null pointer, the value after BVRestoreMat */
1386:       PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Aget);
1387:     }
1388:     MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->n);  /* set the actual pointer */
1389:     *A = bv->Aget;
1390:   }
1391:   PetscFunctionReturn(0);
1392: }

1394: /*@
1395:    BVRestoreMat - Restores the Mat obtained with BVGetMat().

1397:    Logically Collective on bv

1399:    Input Parameters:
1400: +  bv - the basis vectors context
1401: -  A  - the fetched matrix

1403:    Note:
1404:    A call to this function must match a previous call of BVGetMat().
1405:    The effect is that the contents of the Mat are copied back to the
1406:    BV internal data structures.

1408:    Level: advanced

1410: .seealso: BVGetMat(), BVRestoreArray()
1411: @*/
1412: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1413: {
1414:   PetscScalar    *vv,*aa;

1417:   BVCheckSizes(bv,1);
1421:   if (bv->ops->restoremat) (*bv->ops->restoremat)(bv,A);
1422:   else {
1423:     MatDenseGetArray(bv->Aget,&aa);
1424:     vv = aa-(bv->nc+bv->l)*bv->n;
1425:     MatDenseResetArray(bv->Aget);
1426:     BVRestoreArray(bv,&vv);
1427:     *A = NULL;
1428:   }
1429:   PetscFunctionReturn(0);
1430: }

1432: /*
1433:    Copy all user-provided attributes of V to another BV object W
1434:  */
1435: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1436: {
1437:   BVSetType(W,((PetscObject)V)->type_name);
1438:   W->orthog_type  = V->orthog_type;
1439:   W->orthog_ref   = V->orthog_ref;
1440:   W->orthog_eta   = V->orthog_eta;
1441:   W->orthog_block = V->orthog_block;
1442:   if (V->matrix) PetscObjectReference((PetscObject)V->matrix);
1443:   W->matrix       = V->matrix;
1444:   W->indef        = V->indef;
1445:   W->vmm          = V->vmm;
1446:   W->rrandom      = V->rrandom;
1447:   W->deftol       = V->deftol;
1448:   if (V->rand) PetscObjectReference((PetscObject)V->rand);
1449:   W->rand         = V->rand;
1450:   W->sfocalled    = V->sfocalled;
1451:   if (V->ops->duplicate) (*V->ops->duplicate)(V,W);
1452:   PetscObjectStateIncrease((PetscObject)W);
1453:   PetscFunctionReturn(0);
1454: }

1456: /*@
1457:    BVDuplicate - Creates a new basis vector object of the same type and
1458:    dimensions as an existing one.

1460:    Collective on V

1462:    Input Parameter:
1463: .  V - basis vectors context

1465:    Output Parameter:
1466: .  W - location to put the new BV

1468:    Notes:
1469:    The new BV has the same type and dimensions as V, and it shares the same
1470:    template vector. Also, the inner product matrix and orthogonalization
1471:    options are copied.

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

1476:    Level: intermediate

1478: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1479: @*/
1480: PetscErrorCode BVDuplicate(BV V,BV *W)
1481: {
1484:   BVCheckSizes(V,1);
1486:   BVCreate(PetscObjectComm((PetscObject)V),W);
1487:   BVSetSizesFromVec(*W,V->t,V->m);
1488:   BVDuplicate_Private(V,*W);
1489:   PetscFunctionReturn(0);
1490: }

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

1496:    Collective on V

1498:    Input Parameters:
1499: +  V - basis vectors context
1500: -  m - the new number of columns

1502:    Output Parameter:
1503: .  W - location to put the new BV

1505:    Note:
1506:    This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1507:    contents of V are not copied to W.

1509:    Level: intermediate

1511: .seealso: BVDuplicate(), BVResize()
1512: @*/
1513: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1514: {
1517:   BVCheckSizes(V,1);
1520:   BVCreate(PetscObjectComm((PetscObject)V),W);
1521:   BVSetSizesFromVec(*W,V->t,m);
1522:   BVDuplicate_Private(V,*W);
1523:   PetscFunctionReturn(0);
1524: }

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

1530:    Not collective

1532:    Input Parameter:
1533: .  bv    - the basis vectors context

1535:    Output Parameter:
1536: .  cached - the cached BV

1538:    Note:
1539:    The cached BV is created if not available previously.

1541:    Level: developer

1543: .seealso: BVApplyMatrixBV()
1544: @*/
1545: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1546: {
1549:   BVCheckSizes(bv,1);
1550:   if (!bv->cached) {
1551:     BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached);
1552:     BVSetSizesFromVec(bv->cached,bv->t,bv->m);
1553:     BVDuplicate_Private(bv,bv->cached);
1554:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->cached);
1555:   }
1556:   *cached = bv->cached;
1557:   PetscFunctionReturn(0);
1558: }

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

1563:    Logically Collective on V

1565:    Input Parameter:
1566: .  V - basis vectors context

1568:    Output Parameter:
1569: .  W - the copy

1571:    Note:
1572:    Both V and W must be distributed in the same manner; local copies are
1573:    done. Only active columns (excluding the leading ones) are copied.
1574:    In the destination W, columns are overwritten starting from the leading ones.
1575:    Constraints are not copied.

1577:    Level: beginner

1579: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1580: @*/
1581: PetscErrorCode BVCopy(BV V,BV W)
1582: {
1583:   PetscScalar       *womega;
1584:   const PetscScalar *vomega;

1588:   BVCheckSizes(V,1);
1589:   BVCheckOp(V,1,copy);
1592:   BVCheckSizes(W,2);
1596:   if (V==W || !V->n) PetscFunctionReturn(0);

1598:   PetscLogEventBegin(BV_Copy,V,W,0,0);
1599:   if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1600:     /* copy signature */
1601:     BV_AllocateSignature(W);
1602:     VecGetArrayRead(V->omega,&vomega);
1603:     VecGetArray(W->omega,&womega);
1604:     PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l);
1605:     VecRestoreArray(W->omega,&womega);
1606:     VecRestoreArrayRead(V->omega,&vomega);
1607:   }
1608:   (*V->ops->copy)(V,W);
1609:   PetscLogEventEnd(BV_Copy,V,W,0,0);
1610:   PetscObjectStateIncrease((PetscObject)W);
1611:   PetscFunctionReturn(0);
1612: }

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

1617:    Logically Collective on V

1619:    Input Parameters:
1620: +  V - basis vectors context
1621: -  j - the column number to be copied

1623:    Output Parameter:
1624: .  w - the copied column

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

1629:    Level: beginner

1631: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1632: @*/
1633: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1634: {
1635:   PetscInt       n,N;
1636:   Vec            z;

1640:   BVCheckSizes(V,1);

1645:   VecGetSize(w,&N);
1646:   VecGetLocalSize(w,&n);

1649:   PetscLogEventBegin(BV_Copy,V,w,0,0);
1650:   BVGetColumn(V,j,&z);
1651:   VecCopy(z,w);
1652:   BVRestoreColumn(V,j,&z);
1653:   PetscLogEventEnd(BV_Copy,V,w,0,0);
1654:   PetscFunctionReturn(0);
1655: }

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

1660:    Logically Collective on V

1662:    Input Parameters:
1663: +  V - basis vectors context
1664: .  j - the number of the source column
1665: -  i - the number of the destination column

1667:    Level: beginner

1669: .seealso: BVCopy(), BVCopyVec()
1670: @*/
1671: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1672: {
1673:   Vec            z,w;
1674:   PetscScalar    *omega;

1678:   BVCheckSizes(V,1);
1681:   if (j==i) PetscFunctionReturn(0);

1683:   PetscLogEventBegin(BV_Copy,V,0,0,0);
1684:   if (V->omega) {
1685:     VecGetArray(V->omega,&omega);
1686:     omega[i] = omega[j];
1687:     VecRestoreArray(V->omega,&omega);
1688:   }
1689:   if (V->ops->copycolumn) (*V->ops->copycolumn)(V,j,i);
1690:   else {
1691:     BVGetColumn(V,j,&z);
1692:     BVGetColumn(V,i,&w);
1693:     VecCopy(z,w);
1694:     BVRestoreColumn(V,j,&z);
1695:     BVRestoreColumn(V,i,&w);
1696:   }
1697:   PetscLogEventEnd(BV_Copy,V,0,0,0);
1698:   PetscObjectStateIncrease((PetscObject)V);
1699:   PetscFunctionReturn(0);
1700: }

1702: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1703: {
1704:   PetscInt       ncols;

1706:   ncols = left? bv->nc+bv->l: bv->m-bv->l;
1707:   if (*split && ncols!=(*split)->m) BVDestroy(split);
1708:   if (!*split) {
1709:     BVCreate(PetscObjectComm((PetscObject)bv),split);
1710:     PetscLogObjectParent((PetscObject)bv,(PetscObject)*split);
1711:     (*split)->issplit = left? 1: 2;
1712:     (*split)->splitparent = bv;
1713:     BVSetSizesFromVec(*split,bv->t,ncols);
1714:     BVDuplicate_Private(bv,*split);
1715:   }
1716:   (*split)->l  = 0;
1717:   (*split)->k  = left? bv->l: bv->k-bv->l;
1718:   (*split)->nc = left? bv->nc: 0;
1719:   (*split)->m  = ncols-(*split)->nc;
1720:   if ((*split)->nc) {
1721:     (*split)->ci[0] = -(*split)->nc-1;
1722:     (*split)->ci[1] = -(*split)->nc-1;
1723:   }
1724:   if (left) PetscObjectStateGet((PetscObject)*split,&bv->lstate);
1725:   else PetscObjectStateGet((PetscObject)*split,&bv->rstate);
1726:   PetscFunctionReturn(0);
1727: }

1729: /*@
1730:    BVGetSplit - Splits the BV object into two BV objects that share the
1731:    internal data, one of them containing the leading columns and the other
1732:    one containing the remaining columns.

1734:    Logically Collective on bv

1736:    Input Parameter:
1737: .  bv - the basis vectors context

1739:    Output Parameters:
1740: +  L - left BV containing leading columns (can be NULL)
1741: -  R - right BV containing remaining columns (can be NULL)

1743:    Notes:
1744:    The columns are split in two sets. The leading columns (including the
1745:    constraints) are assigned to the left BV and the remaining columns
1746:    are assigned to the right BV. The number of leading columns, as
1747:    specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1748:    guarantee that both L and R have at least one column).

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

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

1757:    Level: advanced

1759: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1760: @*/
1761: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1762: {
1765:   BVCheckSizes(bv,1);
1768:   bv->lsplit = bv->nc+bv->l;
1769:   BVGetSplit_Private(bv,PETSC_TRUE,&bv->L);
1770:   BVGetSplit_Private(bv,PETSC_FALSE,&bv->R);
1771:   if (L) *L = bv->L;
1772:   if (R) *R = bv->R;
1773:   PetscFunctionReturn(0);
1774: }

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

1779:    Logically Collective on bv

1781:    Input Parameters:
1782: +  bv - the basis vectors context
1783: .  L  - left BV obtained with BVGetSplit()
1784: -  R  - right BV obtained with BVGetSplit()

1786:    Note:
1787:    The arguments must match the corresponding call to BVGetSplit().

1789:    Level: advanced

1791: .seealso: BVGetSplit()
1792: @*/
1793: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1794: {
1797:   BVCheckSizes(bv,1);

1804:   if (bv->ops->restoresplit) (*bv->ops->restoresplit)(bv,L,R);
1805:   bv->lsplit = 0;
1806:   if (L) *L = NULL;
1807:   if (R) *R = NULL;
1808:   PetscFunctionReturn(0);
1809: }

1811: /*@
1812:    BVSetDefiniteTolerance - Set the tolerance to be used when checking a
1813:    definite inner product.

1815:    Logically Collective on bv

1817:    Input Parameters:
1818: +  bv     - basis vectors
1819: -  deftol - the tolerance

1821:    Options Database Key:
1822: .  -bv_definite_tol <deftol> - the tolerance

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

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

1836:    Level: advanced

1838: .seealso: BVSetMatrix()
1839: @*/
1840: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
1841: {
1844:   if (deftol == PETSC_DEFAULT) bv->deftol = 10*PETSC_MACHINE_EPSILON;
1845:   else {
1847:     bv->deftol = deftol;
1848:   }
1849:   PetscFunctionReturn(0);
1850: }

1852: /*@
1853:    BVGetDefiniteTolerance - Returns the tolerance for checking a definite
1854:    inner product.

1856:    Not Collective

1858:    Input Parameter:
1859: .  bv - the basis vectors

1861:    Output Parameters:
1862: .  deftol - the tolerance

1864:    Level: advanced

1866: .seealso: BVSetDefiniteTolerance()
1867: @*/
1868: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
1869: {
1872:   *deftol = bv->deftol;
1873:   PetscFunctionReturn(0);
1874: }