Actual source code: bvbasic.c

slepc-3.18.0 2022-10-01
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: /*@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) return 0;
 45:   PetscStrcmp(type,BVTENSOR,&match);

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

 51:   PetscTryTypeMethod(bv,destroy);
 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:   return 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:   return 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:     PetscUseTypeMethod(bv,create);
137:     PetscLogEventEnd(BV_Create,bv,0,0,0);
138:     bv->ops->create = NULL;
139:     bv->defersfo = PETSC_FALSE;
140:   }
141:   return 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:     PetscUseTypeMethod(bv,create);
181:     bv->ops->create = NULL;
182:     bv->defersfo = PETSC_FALSE;
183:   }
184:   return 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:     return 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:   return 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) return 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:   return 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:   return 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) return 0;
336:   BVCheckOp(bv,1,resize);

338:   PetscLogEventBegin(BV_Create,bv,0,0,0);
339:   PetscUseTypeMethod(bv,resize,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:     if (copy) {
352:       VecGetArray(v,&array);
353:       VecGetArrayRead(bv->omega,&omega);
354:       PetscArraycpy(array,omega,PetscMin(m,bv->m));
355:       VecRestoreArrayRead(bv->omega,&omega);
356:       VecRestoreArray(v,&array);
357:     } else VecSet(v,1.0);
358:     VecDestroy(&bv->omega);
359:     bv->omega = v;
360:   }
361:   bv->m = m;
362:   bv->k = PetscMin(bv->k,m);
363:   bv->l = PetscMin(bv->l,m);
364:   PetscLogEventEnd(BV_Create,bv,0,0,0);
365:   PetscObjectStateIncrease((PetscObject)bv);
366:   return 0;
367: }

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

372:    Logically Collective on bv

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

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

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

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

392:    Level: intermediate

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

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

420:    Not Collective

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

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

429:    Level: intermediate

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

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

444:    Collective on bv

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

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

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

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

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

466:    Level: advanced

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

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

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

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

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

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

506:    Level: advanced

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

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

522:    Neighbor-wise Collective on bv

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

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

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

534:    Level: advanced

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

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

554:    Neighbor-wise Collective on X

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

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

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

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

569:    Level: developer

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

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

587:    Logically Collective on bv

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

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

596:    Level: developer

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

607:   BVCheckSizes(bv,1);

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

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

628:    Not collective

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

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

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

639:    Level: developer

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

650:   BVCheckSizes(bv,1);

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

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

670:    Collective on bv

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

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

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

683:    Level: developer

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

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

701:   PetscObjectReference((PetscObject)buffer);
702:   VecDestroy(&bv->buffer);
703:   bv->buffer = buffer;
704:   return 0;
705: }

707: /*@
708:    BVGetBufferVec - Obtain the buffer vector associated with the BV object.

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

712:    Input Parameters:
713: .  bv - the basis vectors context

715:    Output Parameter:
716: .  buffer - vector

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

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

737:    Level: developer

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

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

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

762:    Collective on bv

764:    Input Parameters:
765: +  bv   - the basis vectors context
766: -  rand - the random number generator context

768:    Level: advanced

770: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
771: @*/
772: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
773: {
777:   PetscObjectReference((PetscObject)rand);
778:   PetscRandomDestroy(&bv->rand);
779:   bv->rand = rand;
780:   return 0;
781: }

783: /*@
784:    BVGetRandomContext - Gets the PetscRandom object associated with the BV.

786:    Not Collective

788:    Input Parameter:
789: .  bv - the basis vectors context

791:    Output Parameter:
792: .  rand - the random number generator context

794:    Level: advanced

796: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
797: @*/
798: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
799: {
802:   if (!bv->rand) {
803:     PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand);
804:     if (bv->cuda) PetscRandomSetType(bv->rand,PETSCCURAND);
805:     if (bv->sfocalled) PetscRandomSetFromOptions(bv->rand);
806:     if (bv->rrandom) {
807:       PetscRandomSetSeed(bv->rand,0x12345678);
808:       PetscRandomSeed(bv->rand);
809:     }
810:   }
811:   *rand = bv->rand;
812:   return 0;
813: }

815: /*@
816:    BVSetFromOptions - Sets BV options from the options database.

818:    Collective on bv

820:    Input Parameter:
821: .  bv - the basis vectors context

823:    Level: beginner

825: .seealso: BVSetOptionsPrefix()
826: @*/
827: PetscErrorCode BVSetFromOptions(BV bv)
828: {
829:   char               type[256];
830:   PetscBool          flg1,flg2,flg3,flg4;
831:   PetscReal          r;
832:   BVOrthogType       otype;
833:   BVOrthogRefineType orefine;
834:   BVOrthogBlockType  oblock;

837:   BVRegisterAll();
838:   PetscObjectOptionsBegin((PetscObject)bv);
839:     PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVSVEC),type,sizeof(type),&flg1);
840:     if (flg1) BVSetType(bv,type);
841:     else if (!((PetscObject)bv)->type_name) BVSetType(bv,BVSVEC);

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

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

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

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

861:     if (bv->ops->create) bv->defersfo = PETSC_TRUE;   /* defer call to setfromoptions */
862:     else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
863:     PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject);
864:   PetscOptionsEnd();
865:   bv->sfocalled = PETSC_TRUE;
866:   return 0;
867: }

869: /*@
870:    BVSetOrthogonalization - Specifies the method used for the orthogonalization of
871:    vectors (classical or modified Gram-Schmidt with or without refinement), and
872:    for the block-orthogonalization (simultaneous orthogonalization of a set of
873:    vectors).

875:    Logically Collective on bv

877:    Input Parameters:
878: +  bv     - the basis vectors context
879: .  type   - the method of vector orthogonalization
880: .  refine - type of refinement
881: .  eta    - parameter for selective refinement
882: -  block  - the method of block orthogonalization

884:    Options Database Keys:
885: +  -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
886:                          (default) or mgs for Modified Gram-Schmidt orthogonalization
887: .  -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
888: .  -bv_orthog_eta <eta> -  For setting the value of eta
889: -  -bv_orthog_block <block> - Where <block> is the block-orthogonalization method

891:    Notes:
892:    The default settings work well for most problems.

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

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

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

902:    Level: advanced

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

950: /*@
951:    BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.

953:    Not Collective

955:    Input Parameter:
956: .  bv - basis vectors context

958:    Output Parameters:
959: +  type   - the method of vector orthogonalization
960: .  refine - type of refinement
961: .  eta    - parameter for selective refinement
962: -  block  - the method of block orthogonalization

964:    Level: advanced

966: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
967: @*/
968: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
969: {
971:   if (type)   *type   = bv->orthog_type;
972:   if (refine) *refine = bv->orthog_ref;
973:   if (eta)    *eta    = bv->orthog_eta;
974:   if (block)  *block  = bv->orthog_block;
975:   return 0;
976: }

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

981:    Logically Collective on bv

983:    Input Parameters:
984: +  bv     - the basis vectors context
985: -  method - the method for the BVMatMult() operation

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

990:    Notes:
991:    Allowed values are
992: +  BV_MATMULT_VECS - perform a matrix-vector multiply per each column
993: .  BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
994: -  BV_MATMULT_MAT_SAVE - this case is deprecated

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

998:    Level: advanced

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

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

1024:    Not Collective

1026:    Input Parameter:
1027: .  bv - basis vectors context

1029:    Output Parameter:
1030: .  method - the method for the BVMatMult() operation

1032:    Level: advanced

1034: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1035: @*/
1036: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1037: {
1040:   *method = bv->vmm;
1041:   return 0;
1042: }

1044: /*@
1045:    BVGetColumn - Returns a Vec object that contains the entries of the
1046:    requested column of the basis vectors object.

1048:    Logically Collective on bv

1050:    Input Parameters:
1051: +  bv - the basis vectors context
1052: -  j  - the index of the requested column

1054:    Output Parameter:
1055: .  v  - vector containing the jth column

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

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

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

1069:    Level: beginner

1071: .seealso: BVRestoreColumn(), BVInsertConstraints()
1072: @*/
1073: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1074: {
1075:   PetscInt       l;

1079:   BVCheckSizes(bv,1);
1080:   BVCheckOp(bv,1,getcolumn);
1085:   l = BVAvailableVec;
1087:   PetscUseTypeMethod(bv,getcolumn,j,v);
1088:   bv->ci[l] = j;
1089:   PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]);
1090:   PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]);
1091:   *v = bv->cv[l];
1092:   return 0;
1093: }

1095: /*@
1096:    BVRestoreColumn - Restore a column obtained with BVGetColumn().

1098:    Logically Collective on bv

1100:    Input Parameters:
1101: +  bv - the basis vectors context
1102: .  j  - the index of the column
1103: -  v  - vector obtained with BVGetColumn()

1105:    Note:
1106:    The arguments must match the corresponding call to BVGetColumn().

1108:    Level: beginner

1110: .seealso: BVGetColumn()
1111: @*/
1112: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1113: {
1114:   PetscObjectId    id;
1115:   PetscObjectState st;
1116:   PetscInt         l;

1120:   BVCheckSizes(bv,1);
1127:   l = (j==bv->ci[0])? 0: 1;
1128:   PetscObjectGetId((PetscObject)*v,&id);
1130:   PetscObjectStateGet((PetscObject)*v,&st);
1131:   if (st!=bv->st[l]) PetscObjectStateIncrease((PetscObject)bv);
1132:   PetscUseTypeMethod(bv,restorecolumn,j,v);
1133:   bv->ci[l] = -bv->nc-1;
1134:   bv->st[l] = -1;
1135:   bv->id[l] = 0;
1136:   *v = NULL;
1137:   return 0;
1138: }

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

1144:    Logically Collective on bv

1146:    Input Parameters:
1147: .  bv - the basis vectors context

1149:    Output Parameter:
1150: .  a  - location to put pointer to the array

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

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

1160:    Level: advanced

1162: .seealso: BVRestoreArray(), BVInsertConstraints()
1163: @*/
1164: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1165: {
1168:   BVCheckSizes(bv,1);
1169:   BVCheckOp(bv,1,getarray);
1170:   PetscUseTypeMethod(bv,getarray,a);
1171:   return 0;
1172: }

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

1177:    Logically Collective on bv

1179:    Input Parameters:
1180: +  bv - the basis vectors context
1181: -  a  - location of pointer to array obtained from BVGetArray()

1183:    Note:
1184:    This operation may imply a data copy, for BV types that do not store
1185:    data contiguously in memory.

1187:    Level: advanced

1189: .seealso: BVGetColumn()
1190: @*/
1191: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1192: {
1195:   BVCheckSizes(bv,1);
1196:   PetscTryTypeMethod(bv,restorearray,a);
1197:   if (a) *a = NULL;
1198:   PetscObjectStateIncrease((PetscObject)bv);
1199:   return 0;
1200: }

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

1206:    Not Collective

1208:    Input Parameters:
1209: .  bv - the basis vectors context

1211:    Output Parameter:
1212: .  a  - location to put pointer to the array

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

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

1222:    Level: advanced

1224: .seealso: BVRestoreArray(), BVInsertConstraints()
1225: @*/
1226: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1227: {
1230:   BVCheckSizes(bv,1);
1231:   BVCheckOp(bv,1,getarrayread);
1232:   PetscUseTypeMethod(bv,getarrayread,a);
1233:   return 0;
1234: }

1236: /*@C
1237:    BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1238:    been called.

1240:    Not Collective

1242:    Input Parameters:
1243: +  bv - the basis vectors context
1244: -  a  - location of pointer to array obtained from BVGetArrayRead()

1246:    Level: advanced

1248: .seealso: BVGetColumn()
1249: @*/
1250: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1251: {
1254:   BVCheckSizes(bv,1);
1255:   PetscTryTypeMethod(bv,restorearrayread,a);
1256:   if (a) *a = NULL;
1257:   return 0;
1258: }

1260: /*@
1261:    BVCreateVec - Creates a new Vec object with the same type and dimensions
1262:    as the columns of the basis vectors object.

1264:    Collective on bv

1266:    Input Parameter:
1267: .  bv - the basis vectors context

1269:    Output Parameter:
1270: .  v  - the new vector

1272:    Note:
1273:    The user is responsible of destroying the returned vector.

1275:    Level: beginner

1277: .seealso: BVCreateMat()
1278: @*/
1279: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1280: {
1282:   BVCheckSizes(bv,1);
1284:   VecDuplicate(bv->t,v);
1285:   return 0;
1286: }

1288: /*@
1289:    BVCreateMat - Creates a new Mat object of dense type and copies the contents
1290:    of the BV object.

1292:    Collective on bv

1294:    Input Parameter:
1295: .  bv - the basis vectors context

1297:    Output Parameter:
1298: .  A  - the new matrix

1300:    Notes:
1301:    The user is responsible of destroying the returned matrix.

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

1305:    Level: intermediate

1307: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1308: @*/
1309: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1310: {
1311:   PetscScalar       *aa;
1312:   const PetscScalar *vv;

1315:   BVCheckSizes(bv,1);

1318:   MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,bv->N,bv->m,NULL,A);
1319:   MatDenseGetArrayWrite(*A,&aa);
1320:   BVGetArrayRead(bv,&vv);
1321:   PetscArraycpy(aa,vv,bv->m*bv->n);
1322:   BVRestoreArrayRead(bv,&vv);
1323:   MatDenseRestoreArrayWrite(*A,&aa);
1324:   return 0;
1325: }

1327: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1328: {
1329:   PetscScalar *vv,*aa;
1330:   PetscBool   create=PETSC_FALSE;
1331:   PetscInt    m,cols;

1333:   m = bv->k-bv->l;
1334:   if (!bv->Aget) create=PETSC_TRUE;
1335:   else {
1336:     MatDenseGetArray(bv->Aget,&aa);
1338:     MatGetSize(bv->Aget,NULL,&cols);
1339:     if (cols!=m) {
1340:       MatDestroy(&bv->Aget);
1341:       create=PETSC_TRUE;
1342:     }
1343:   }
1344:   BVGetArray(bv,&vv);
1345:   if (create) {
1346:     MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,vv,&bv->Aget); /* pass a pointer to avoid allocation of storage */
1347:     MatDenseReplaceArray(bv->Aget,NULL);  /* replace with a null pointer, the value after BVRestoreMat */
1348:   }
1349:   MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->n);  /* set the actual pointer */
1350:   *A = bv->Aget;
1351:   return 0;
1352: }

1354: /*@
1355:    BVGetMat - Returns a Mat object of dense type that shares the memory of
1356:    the BV object.

1358:    Collective on bv

1360:    Input Parameter:
1361: .  bv - the basis vectors context

1363:    Output Parameter:
1364: .  A  - the matrix

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

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

1374:    Level: advanced

1376: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1377: @*/
1378: PetscErrorCode BVGetMat(BV bv,Mat *A)
1379: {
1381:   BVCheckSizes(bv,1);
1383:   PetscUseTypeMethod(bv,getmat,A);
1384:   return 0;
1385: }

1387: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1388: {
1389:   PetscScalar *vv,*aa;

1391:   MatDenseGetArray(bv->Aget,&aa);
1392:   vv = aa-(bv->nc+bv->l)*bv->n;
1393:   MatDenseResetArray(bv->Aget);
1394:   BVRestoreArray(bv,&vv);
1395:   *A = NULL;
1396:   return 0;
1397: }

1399: /*@
1400:    BVRestoreMat - Restores the Mat obtained with BVGetMat().

1402:    Logically Collective on bv

1404:    Input Parameters:
1405: +  bv - the basis vectors context
1406: -  A  - the fetched matrix

1408:    Note:
1409:    A call to this function must match a previous call of BVGetMat().
1410:    The effect is that the contents of the Mat are copied back to the
1411:    BV internal data structures.

1413:    Level: advanced

1415: .seealso: BVGetMat(), BVRestoreArray()
1416: @*/
1417: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1418: {
1420:   BVCheckSizes(bv,1);
1424:   PetscUseTypeMethod(bv,restoremat,A);
1425:   return 0;
1426: }

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

1452: /*@
1453:    BVDuplicate - Creates a new basis vector object of the same type and
1454:    dimensions as an existing one.

1456:    Collective on V

1458:    Input Parameter:
1459: .  V - basis vectors context

1461:    Output Parameter:
1462: .  W - location to put the new BV

1464:    Notes:
1465:    The new BV has the same type and dimensions as V, and it shares the same
1466:    template vector. Also, the inner product matrix and orthogonalization
1467:    options are copied.

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

1472:    Level: intermediate

1474: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1475: @*/
1476: PetscErrorCode BVDuplicate(BV V,BV *W)
1477: {
1480:   BVCheckSizes(V,1);
1482:   BVCreate(PetscObjectComm((PetscObject)V),W);
1483:   BVSetSizesFromVec(*W,V->t,V->m);
1484:   BVDuplicate_Private(V,*W);
1485:   return 0;
1486: }

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

1492:    Collective on V

1494:    Input Parameters:
1495: +  V - basis vectors context
1496: -  m - the new number of columns

1498:    Output Parameter:
1499: .  W - location to put the new BV

1501:    Note:
1502:    This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1503:    contents of V are not copied to W.

1505:    Level: intermediate

1507: .seealso: BVDuplicate(), BVResize()
1508: @*/
1509: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1510: {
1513:   BVCheckSizes(V,1);
1516:   BVCreate(PetscObjectComm((PetscObject)V),W);
1517:   BVSetSizesFromVec(*W,V->t,m);
1518:   BVDuplicate_Private(V,*W);
1519:   return 0;
1520: }

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

1526:    Not collective

1528:    Input Parameter:
1529: .  bv    - the basis vectors context

1531:    Output Parameter:
1532: .  cached - the cached BV

1534:    Note:
1535:    The cached BV is created if not available previously.

1537:    Level: developer

1539: .seealso: BVApplyMatrixBV()
1540: @*/
1541: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1542: {
1545:   BVCheckSizes(bv,1);
1546:   if (!bv->cached) {
1547:     BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached);
1548:     BVSetSizesFromVec(bv->cached,bv->t,bv->m);
1549:     BVDuplicate_Private(bv,bv->cached);
1550:   }
1551:   *cached = bv->cached;
1552:   return 0;
1553: }

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

1558:    Logically Collective on V

1560:    Input Parameter:
1561: .  V - basis vectors context

1563:    Output Parameter:
1564: .  W - the copy

1566:    Note:
1567:    Both V and W must be distributed in the same manner; local copies are
1568:    done. Only active columns (excluding the leading ones) are copied.
1569:    In the destination W, columns are overwritten starting from the leading ones.
1570:    Constraints are not copied.

1572:    Level: beginner

1574: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1575: @*/
1576: PetscErrorCode BVCopy(BV V,BV W)
1577: {
1578:   PetscScalar       *womega;
1579:   const PetscScalar *vomega;

1583:   BVCheckSizes(V,1);
1584:   BVCheckOp(V,1,copy);
1587:   BVCheckSizes(W,2);
1591:   if (V==W || !V->n) return 0;

1593:   PetscLogEventBegin(BV_Copy,V,W,0,0);
1594:   if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1595:     /* copy signature */
1596:     BV_AllocateSignature(W);
1597:     VecGetArrayRead(V->omega,&vomega);
1598:     VecGetArray(W->omega,&womega);
1599:     PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l);
1600:     VecRestoreArray(W->omega,&womega);
1601:     VecRestoreArrayRead(V->omega,&vomega);
1602:   }
1603:   PetscUseTypeMethod(V,copy,W);
1604:   PetscLogEventEnd(BV_Copy,V,W,0,0);
1605:   PetscObjectStateIncrease((PetscObject)W);
1606:   return 0;
1607: }

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

1612:    Logically Collective on V

1614:    Input Parameters:
1615: +  V - basis vectors context
1616: -  j - the column number to be copied

1618:    Output Parameter:
1619: .  w - the copied column

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

1624:    Level: beginner

1626: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1627: @*/
1628: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1629: {
1630:   PetscInt       n,N;
1631:   Vec            z;

1635:   BVCheckSizes(V,1);

1640:   VecGetSize(w,&N);
1641:   VecGetLocalSize(w,&n);

1644:   PetscLogEventBegin(BV_Copy,V,w,0,0);
1645:   BVGetColumn(V,j,&z);
1646:   VecCopy(z,w);
1647:   BVRestoreColumn(V,j,&z);
1648:   PetscLogEventEnd(BV_Copy,V,w,0,0);
1649:   return 0;
1650: }

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

1655:    Logically Collective on V

1657:    Input Parameters:
1658: +  V - basis vectors context
1659: .  j - the number of the source column
1660: -  i - the number of the destination column

1662:    Level: beginner

1664: .seealso: BVCopy(), BVCopyVec()
1665: @*/
1666: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1667: {
1668:   PetscScalar *omega;

1672:   BVCheckSizes(V,1);
1675:   if (j==i) return 0;

1677:   PetscLogEventBegin(BV_Copy,V,0,0,0);
1678:   if (V->omega) {
1679:     VecGetArray(V->omega,&omega);
1680:     omega[i] = omega[j];
1681:     VecRestoreArray(V->omega,&omega);
1682:   }
1683:   PetscUseTypeMethod(V,copycolumn,j,i);
1684:   PetscLogEventEnd(BV_Copy,V,0,0,0);
1685:   PetscObjectStateIncrease((PetscObject)V);
1686:   return 0;
1687: }

1689: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1690: {
1691:   PetscInt       ncols;

1693:   ncols = left? bv->nc+bv->l: bv->m-bv->l;
1694:   if (*split && ncols!=(*split)->m) BVDestroy(split);
1695:   if (!*split) {
1696:     BVCreate(PetscObjectComm((PetscObject)bv),split);
1697:     (*split)->issplit = left? 1: 2;
1698:     (*split)->splitparent = bv;
1699:     BVSetSizesFromVec(*split,bv->t,ncols);
1700:     BVDuplicate_Private(bv,*split);
1701:   }
1702:   (*split)->l  = 0;
1703:   (*split)->k  = left? bv->l: bv->k-bv->l;
1704:   (*split)->nc = left? bv->nc: 0;
1705:   (*split)->m  = ncols-(*split)->nc;
1706:   if ((*split)->nc) {
1707:     (*split)->ci[0] = -(*split)->nc-1;
1708:     (*split)->ci[1] = -(*split)->nc-1;
1709:   }
1710:   if (left) PetscObjectStateGet((PetscObject)*split,&bv->lstate);
1711:   else PetscObjectStateGet((PetscObject)*split,&bv->rstate);
1712:   return 0;
1713: }

1715: /*@
1716:    BVGetSplit - Splits the BV object into two BV objects that share the
1717:    internal data, one of them containing the leading columns and the other
1718:    one containing the remaining columns.

1720:    Logically Collective on bv

1722:    Input Parameter:
1723: .  bv - the basis vectors context

1725:    Output Parameters:
1726: +  L - left BV containing leading columns (can be NULL)
1727: -  R - right BV containing remaining columns (can be NULL)

1729:    Notes:
1730:    The columns are split in two sets. The leading columns (including the
1731:    constraints) are assigned to the left BV and the remaining columns
1732:    are assigned to the right BV. The number of leading columns, as
1733:    specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1734:    guarantee that both L and R have at least one column).

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

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

1743:    Level: advanced

1745: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1746: @*/
1747: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1748: {
1751:   BVCheckSizes(bv,1);
1754:   bv->lsplit = bv->nc+bv->l;
1755:   BVGetSplit_Private(bv,PETSC_TRUE,&bv->L);
1756:   BVGetSplit_Private(bv,PETSC_FALSE,&bv->R);
1757:   if (L) *L = bv->L;
1758:   if (R) *R = bv->R;
1759:   return 0;
1760: }

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

1765:    Logically Collective on bv

1767:    Input Parameters:
1768: +  bv - the basis vectors context
1769: .  L  - left BV obtained with BVGetSplit()
1770: -  R  - right BV obtained with BVGetSplit()

1772:    Note:
1773:    The arguments must match the corresponding call to BVGetSplit().

1775:    Level: advanced

1777: .seealso: BVGetSplit()
1778: @*/
1779: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1780: {
1783:   BVCheckSizes(bv,1);

1790:   PetscTryTypeMethod(bv,restoresplit,L,R);
1791:   bv->lsplit = 0;
1792:   if (L) *L = NULL;
1793:   if (R) *R = NULL;
1794:   return 0;
1795: }

1797: /*@
1798:    BVSetDefiniteTolerance - Set the tolerance to be used when checking a
1799:    definite inner product.

1801:    Logically Collective on bv

1803:    Input Parameters:
1804: +  bv     - basis vectors
1805: -  deftol - the tolerance

1807:    Options Database Key:
1808: .  -bv_definite_tol <deftol> - the tolerance

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

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

1822:    Level: advanced

1824: .seealso: BVSetMatrix()
1825: @*/
1826: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
1827: {
1830:   if (deftol == PETSC_DEFAULT) bv->deftol = 10*PETSC_MACHINE_EPSILON;
1831:   else {
1833:     bv->deftol = deftol;
1834:   }
1835:   return 0;
1836: }

1838: /*@
1839:    BVGetDefiniteTolerance - Returns the tolerance for checking a definite
1840:    inner product.

1842:    Not Collective

1844:    Input Parameter:
1845: .  bv - the basis vectors

1847:    Output Parameters:
1848: .  deftol - the tolerance

1850:    Level: advanced

1852: .seealso: BVSetDefiniteTolerance()
1853: @*/
1854: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
1855: {
1858:   *deftol = bv->deftol;
1859:   return 0;
1860: }