Actual source code: bvbasic.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>      /*I "slepcbv.h" I*/

 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 Parameter:
 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 ierr,(*r)(BV);
 38:   PetscBool      match;


 44:   PetscObjectTypeCompare((PetscObject)bv,type,&match);
 45:   if (match) return(0);
 46:   PetscStrcmp(type,BVTENSOR,&match);
 47:   if (match) SETERRQ(PetscObjectComm((PetscObject)bv),1,"Use BVCreateTensor() to create a BV of type tensor");

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

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

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

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

 69:    Not Collective

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

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

 77:    Level: intermediate

 79: .seealso: BVSetType()
 80: @*/
 81: PetscErrorCode BVGetType(BV bv,BVType *type)
 82: {
 86:   *type = ((PetscObject)bv)->type_name;
 87:   return(0);
 88: }

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

 93:    Collective on BV

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

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

106:    Level: beginner

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

119:   if (N >= 0 && n > N) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local size %D cannot be larger than global size %D",n,N);
120:   if (m <= 0) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %D must be positive",m);
121:   if ((bv->n >= 0 || bv->N >= 0) && (bv->n != n || bv->N != N)) SETERRQ4(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change/reset vector sizes to %D local %D global after previously setting them to %D local %D global",n,N,bv->n,bv->N);
122:   if (bv->m > 0 && bv->m != m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change the number of columns to %D after previously setting it to %D; use BVResize()",m,bv->m);
123:   bv->n = n;
124:   bv->N = N;
125:   bv->m = m;
126:   bv->k = m;
127:   if (!bv->t) {  /* create template vector and get actual dimensions */
128:     VecCreate(PetscObjectComm((PetscObject)bv),&bv->t);
129:     VecSetSizes(bv->t,bv->n,bv->N);
130:     VecSetFromOptions(bv->t);
131:     VecGetSize(bv->t,&bv->N);
132:     VecGetLocalSize(bv->t,&bv->n);
133:     if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
134:       MatGetLocalSize(bv->matrix,&ma,NULL);
135:       if (bv->n!=ma) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %D does not match that of matrix given at BVSetMatrix %D",bv->n,ma);
136:     }
137:   }
138:   if (bv->ops->create) {
139:     PetscLogEventBegin(BV_Create,bv,0,0,0);
140:     (*bv->ops->create)(bv);
141:     PetscLogEventEnd(BV_Create,bv,0,0,0);
142:     bv->ops->create = 0;
143:     bv->defersfo = PETSC_FALSE;
144:   }
145:   return(0);
146: }

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

152:    Collective on BV

154:    Input Parameters:
155: +  bv - the basis vectors
156: .  t  - the template vector
157: -  m  - the number of columns

159:    Level: beginner

161: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
162: @*/
163: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
164: {
166:   PetscInt       ma;

173:   if (m <= 0) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %D must be positive",m);
174:   if (bv->t) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Template vector was already set by a previous call to BVSetSizes/FromVec");
175:   VecGetSize(t,&bv->N);
176:   VecGetLocalSize(t,&bv->n);
177:   if (bv->matrix) {  /* check compatible dimensions of user-provided matrix */
178:     MatGetLocalSize(bv->matrix,&ma,NULL);
179:     if (bv->n!=ma) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %D does not match that of matrix given at BVSetMatrix %D",bv->n,ma);
180:   }
181:   bv->m = m;
182:   bv->k = m;
183:   bv->t = t;
184:   PetscObjectReference((PetscObject)t);
185:   if (bv->ops->create) {
186:     (*bv->ops->create)(bv);
187:     bv->ops->create = 0;
188:     bv->defersfo = PETSC_FALSE;
189:   }
190:   return(0);
191: }

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

196:    Not Collective

198:    Input Parameter:
199: .  bv - the basis vectors

201:    Output Parameters:
202: +  n  - the local size
203: .  N  - the global size
204: -  m  - the number of columns

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

213:    Level: beginner

215: .seealso: BVSetSizes(), BVSetSizesFromVec()
216: @*/
217: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
218: {
220:   if (!bv) {
221:     if (m && !n && !N) *m = 0;
222:     return(0);
223:   }
225:   if (n || N) BVCheckSizes(bv,1);
226:   if (n) *n = bv->n;
227:   if (N) *N = bv->N;
228:   if (m) *m = bv->m;
229:   if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
230:   return(0);
231: }

233: /*@
234:    BVSetNumConstraints - Set the number of constraints.

236:    Logically Collective on BV

238:    Input Parameters:
239: +  V  - basis vectors
240: -  nc - number of constraints

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

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

250:    Level: developer

252: .seealso: BVInsertConstraints()
253: @*/
254: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
255: {
257:   PetscInt       total,diff,i;
258:   Vec            x,y;

263:   if (nc<0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %D) cannot be negative",nc);
265:   BVCheckSizes(V,1);
266:   if (V->ci[0]!=-V->nc-1 || V->ci[1]!=-V->nc-1) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");

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

291: /*@
292:    BVGetNumConstraints - Returns the number of constraints.

294:    Not Collective

296:    Input Parameter:
297: .  bv - the basis vectors

299:    Output Parameters:
300: .  nc - the number of constraints

302:    Level: advanced

304: .seealso: BVGetSizes(), BVInsertConstraints()
305: @*/
306: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
307: {
311:   *nc = bv->nc;
312:   return(0);
313: }

315: /*@
316:    BVResize - Change the number of columns.

318:    Collective on BV

320:    Input Parameters:
321: +  bv   - the basis vectors
322: .  m    - the new number of columns
323: -  copy - a flag indicating whether current values should be kept

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

329:    Level: advanced

331: .seealso: BVSetSizes(), BVSetSizesFromVec()
332: @*/
333: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
334: {
335:   PetscErrorCode    ierr;
336:   PetscScalar       *array;
337:   const PetscScalar *omega;
338:   Vec               v;

345:   if (m <= 0) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %D must be positive",m);
346:   if (bv->nc && !bv->issplit) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
347:   if (bv->m == m) return(0);

349:   PetscLogEventBegin(BV_Create,bv,0,0,0);
350:   (*bv->ops->resize)(bv,m,copy);
351:   VecDestroy(&bv->buffer);
352:   BVDestroy(&bv->cached);
353:   PetscFree2(bv->h,bv->c);
354:   if (bv->omega) {
355:     if (bv->cuda) {
356: #if defined(PETSC_HAVE_CUDA)
357:       VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v);
358: #else
359:       SETERRQ(PetscObjectComm((PetscObject)bv),1,"Something wrong happened");
360: #endif
361:     } else {
362:       VecCreateSeq(PETSC_COMM_SELF,m,&v);
363:     }
364:     PetscLogObjectParent((PetscObject)bv,(PetscObject)v);
365:     if (copy) {
366:       VecGetArray(v,&array);
367:       VecGetArrayRead(bv->omega,&omega);
368:       PetscMemcpy(array,omega,PetscMin(m,bv->m)*sizeof(PetscScalar));
369:       VecRestoreArrayRead(bv->omega,&omega);
370:       VecRestoreArray(v,&array);
371:     } else {
372:       VecSet(v,1.0);
373:     }
374:     VecDestroy(&bv->omega);
375:     bv->omega = v;
376:   }
377:   bv->m = m;
378:   bv->k = PetscMin(bv->k,m);
379:   bv->l = PetscMin(bv->l,m);
380:   PetscLogEventEnd(BV_Create,bv,0,0,0);
381:   PetscObjectStateIncrease((PetscObject)bv);
382:   return(0);
383: }

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

388:    Logically Collective on BV

390:    Input Parameters:
391: +  bv - the basis vectors context
392: .  l  - number of leading columns
393: -  k  - number of active columns

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

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

404:    In orthogonalization operations, the first l columns are treated
405:    differently: they participate in the orthogonalization but the computed
406:    coefficients are not stored.

408:    Level: intermediate

410: .seealso: BVGetActiveColumns(), BVSetSizes()
411: @*/
412: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
413: {
418:   BVCheckSizes(bv,1);
419:   if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
420:     bv->k = bv->m;
421:   } else {
422:     if (k<0 || k>bv->m) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and m");
423:     bv->k = k;
424:   }
425:   if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
426:     bv->l = 0;
427:   } else {
428:     if (l<0 || l>bv->k) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and k");
429:     bv->l = l;
430:   }
431:   return(0);
432: }

434: /*@
435:    BVGetActiveColumns - Returns the current active dimensions.

437:    Not Collective

439:    Input Parameter:
440: .  bv - the basis vectors context

442:    Output Parameter:
443: +  l  - number of leading columns
444: -  k  - number of active columns

446:    Level: intermediate

448: .seealso: BVSetActiveColumns()
449: @*/
450: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
451: {
454:   if (l) *l = bv->l;
455:   if (k) *k = bv->k;
456:   return(0);
457: }

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

462:    Collective on BV

464:    Input Parameters:
465: +  bv    - the basis vectors context
466: .  B     - a symmetric matrix (may be NULL)
467: -  indef - a flag indicating if the matrix is indefinite

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

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

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

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

484:    Level: advanced

486: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize()
487: @*/
488: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
489: {
491:   PetscInt       m,n;

496:   if (B) {
498:     MatGetLocalSize(B,&m,&n);
499:     if (m!=n) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
500:     if (bv->m && bv->n!=n) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension BV %D, Mat %D",bv->n,n);
501:   }
502:   if (B) { PetscObjectReference((PetscObject)B); }
503:   MatDestroy(&bv->matrix);
504:   bv->matrix = B;
505:   bv->indef  = indef;
506:   PetscObjectStateIncrease((PetscObject)bv);
507:   return(0);
508: }

510: /*@
511:    BVGetMatrix - Retrieves the matrix representation of the inner product.

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

515:    Input Parameter:
516: .  bv    - the basis vectors context

518:    Output Parameter:
519: +  B     - the matrix of the inner product (may be NULL)
520: -  indef - the flag indicating if the matrix is indefinite

522:    Level: advanced

524: .seealso: BVSetMatrix()
525: @*/
526: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
527: {
530:   if (B)     *B     = bv->matrix;
531:   if (indef) *indef = bv->indef;
532:   return(0);
533: }

535: /*@
536:    BVApplyMatrix - Multiplies a vector by the matrix representation of the
537:    inner product.

539:    Neighbor-wise Collective on BV and Vec

541:    Input Parameter:
542: +  bv - the basis vectors context
543: -  x  - the vector

545:    Output Parameter:
546: .  y  - the result

548:    Note:
549:    If no matrix was specified this function copies the vector.

551:    Level: advanced

553: .seealso: BVSetMatrix(), BVApplyMatrixBV()
554: @*/
555: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
556: {

563:   if (bv->matrix) {
564:     BV_IPMatMult(bv,x);
565:     VecCopy(bv->Bx,y);
566:   } else {
567:     VecCopy(x,y);
568:   }
569:   return(0);
570: }

572: /*@
573:    BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
574:    of the inner product.

576:    Neighbor-wise Collective on BV

578:    Input Parameter:
579: .  X - the basis vectors context

581:    Output Parameter:
582: .  Y - the basis vectors to store the result (optional)

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

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

591:    Level: developer

593: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
594: @*/
595: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
596: {

601:   if (Y) {
603:     if (X->matrix) {
604:       BVMatMult(X,X->matrix,Y);
605:     } else {
606:       BVCopy(X,Y);
607:     }
608:   } else {
609:     BV_IPMatMultBV(X);
610:   }
611:   return(0);
612: }

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

617:    Logically Collective on BV

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

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

626:    Level: developer

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

639:   BVCheckSizes(bv,1);

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

659: /*@
660:    BVGetSignature - Retrieves the signature matrix from last orthogonalization.

662:    Not collective

664:    Input Parameter:
665: .  bv    - the basis vectors context

667:    Output Parameter:
668: .  omega - a vector representing the diagonal of the signature matrix

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

673:    Level: developer

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

686:   BVCheckSizes(bv,1);

690:   VecGetSize(omega,&n);
691:   if (n!=bv->k) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %D elements, should be %D",n,bv->k);
692:   if (bv->indef && bv->omega) {
693:     VecGetArray(omega,&pomega);
694:     VecGetArrayRead(bv->omega,&intern);
695:     for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
696:     VecRestoreArrayRead(bv->omega,&intern);
697:     VecRestoreArray(omega,&pomega);
698:   } else {
699:     VecSet(omega,1.0);
700:   }
701:   return(0);
702: }

704: /*@
705:    BVSetBufferVec - Attach a vector object to be used as buffer space for
706:    several operations.

708:    Collective on BV

710:    Input Parameters:
711: +  bv     - the basis vectors context)
712: -  buffer - the vector

714:    Notes:
715:    Use BVGetBufferVec() to retrieve the vector (for example, to free it
716:    at the end of the computations).

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

721:    Level: developer

723: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
724: @*/
725: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
726: {
728:   PetscInt       ld,n;
729:   PetscMPIInt    size;

734:   BVCheckSizes(bv,1);
735:   VecGetSize(buffer,&n);
736:   ld = bv->m+bv->nc;
737:   if (n != ld*bv->m) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %d",ld*bv->m);
738:   MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size);
739:   if (size>1) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
740:   
741:   PetscObjectReference((PetscObject)buffer);
742:   VecDestroy(&bv->buffer);
743:   bv->buffer = buffer;
744:   PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
745:   return(0);
746: }

748: /*@
749:    BVGetBufferVec - Obtain the buffer vector associated with the BV object.

751:    Not Collective

753:    Input Parameters:
754: .  bv - the basis vectors context

756:    Output Parameter:
757: .  buffer - vector

759:    Notes:
760:    The vector is created if not available previously. It is a sequential vector
761:    of length (nc+m)*m, where m is the number of columns of bv and nc is the number
762:    of constraints.

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

778:    Level: developer

780: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
781: @*/
782: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
783: {
785:   PetscInt       ld;

790:   BVCheckSizes(bv,1);
791:   if (!bv->buffer) {
792:     ld = bv->m+bv->nc;
793:     VecCreate(PETSC_COMM_SELF,&bv->buffer);
794:     VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m);
795:     VecSetType(bv->buffer,((PetscObject)bv->t)->type_name);
796:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
797:   }
798:   *buffer = bv->buffer;
799:   return(0);
800: }

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

806:    Collective on BV

808:    Input Parameters:
809: +  bv   - the basis vectors context
810: -  rand - the random number generator context

812:    Level: advanced

814: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomColumn()
815: @*/
816: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
817: {

824:   PetscObjectReference((PetscObject)rand);
825:   PetscRandomDestroy(&bv->rand);
826:   bv->rand = rand;
827:   PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
828:   return(0);
829: }

831: /*@
832:    BVGetRandomContext - Gets the PetscRandom object associated with the BV.

834:    Not Collective

836:    Input Parameter:
837: .  bv - the basis vectors context

839:    Output Parameter:
840: .  rand - the random number generator context

842:    Level: advanced

844: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn()
845: @*/
846: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
847: {

853:   if (!bv->rand) {
854:     PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand);
855:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
856:     if (bv->rrandom) {
857:       PetscRandomSetSeed(bv->rand,0x12345678);
858:       PetscRandomSeed(bv->rand);
859:     }
860:   }
861:   *rand = bv->rand;
862:   return(0);
863: }

865: /*@
866:    BVSetFromOptions - Sets BV options from the options database.

868:    Collective on BV

870:    Input Parameter:
871: .  bv - the basis vectors context

873:    Level: beginner
874: @*/
875: PetscErrorCode BVSetFromOptions(BV bv)
876: {
877:   PetscErrorCode     ierr;
878:   char               type[256];
879:   PetscBool          flg1,flg2,flg3,flg4;
880:   PetscReal          r;
881:   BVOrthogType       otype;
882:   BVOrthogRefineType orefine;
883:   BVOrthogBlockType  oblock;

887:   BVRegisterAll();
888:   PetscObjectOptionsBegin((PetscObject)bv);
889:     PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVSVEC),type,256,&flg1);
890:     if (flg1) {
891:       BVSetType(bv,type);
892:     } else if (!((PetscObject)bv)->type_name) {
893:       BVSetType(bv,BVSVEC);
894:     }

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

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

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

911:     if (bv->ops->create) bv->defersfo = PETSC_TRUE;   /* defer call to setfromoptions */
912:     else if (bv->ops->setfromoptions) {
913:       (*bv->ops->setfromoptions)(PetscOptionsObject,bv);
914:     }
915:     PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)bv);
916:   PetscOptionsEnd();

918:   if (!bv->rand) { BVGetRandomContext(bv,&bv->rand); }
919:   PetscRandomSetFromOptions(bv->rand);
920:   return(0);
921: }

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

929:    Logically Collective on BV

931:    Input Parameters:
932: +  bv     - the basis vectors context
933: .  type   - the method of vector orthogonalization
934: .  refine - type of refinement
935: .  eta    - parameter for selective refinement
936: -  block  - the method of block orthogonalization

938:    Options Database Keys:
939: +  -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
940:                          (default) or mgs for Modified Gram-Schmidt orthogonalization
941: .  -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
942: .  -bv_orthog_eta <eta> -  For setting the value of eta
943: -  -bv_orthog_block <block> - Where <block> is the block-orthogonalization method

945:    Notes:
946:    The default settings work well for most problems.

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

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

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

956:    Level: advanced

958: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
959: @*/
960: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
961: {
968:   switch (type) {
969:     case BV_ORTHOG_CGS:
970:     case BV_ORTHOG_MGS:
971:       bv->orthog_type = type;
972:       break;
973:     default:
974:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
975:   }
976:   switch (refine) {
977:     case BV_ORTHOG_REFINE_NEVER:
978:     case BV_ORTHOG_REFINE_IFNEEDED:
979:     case BV_ORTHOG_REFINE_ALWAYS:
980:       bv->orthog_ref = refine;
981:       break;
982:     default:
983:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
984:   }
985:   if (eta == PETSC_DEFAULT) {
986:     bv->orthog_eta = 0.7071;
987:   } else {
988:     if (eta <= 0.0 || eta > 1.0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
989:     bv->orthog_eta = eta;
990:   }
991:   switch (block) {
992:     case BV_ORTHOG_BLOCK_GS:
993:     case BV_ORTHOG_BLOCK_CHOL:
994:     case BV_ORTHOG_BLOCK_TSQR:
995:     case BV_ORTHOG_BLOCK_TSQRCHOL:
996:     case BV_ORTHOG_BLOCK_SVQB:
997:       bv->orthog_block = block;
998:       break;
999:     default:
1000:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
1001:   }
1002:   return(0);
1003: }

1005: /*@
1006:    BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.

1008:    Not Collective

1010:    Input Parameter:
1011: .  bv - basis vectors context

1013:    Output Parameter:
1014: +  type   - the method of vector orthogonalization
1015: .  refine - type of refinement
1016: .  eta    - parameter for selective refinement
1017: -  block  - the method of block orthogonalization

1019:    Level: advanced

1021: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
1022: @*/
1023: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1024: {
1027:   if (type)   *type   = bv->orthog_type;
1028:   if (refine) *refine = bv->orthog_ref;
1029:   if (eta)    *eta    = bv->orthog_eta;
1030:   if (block)  *block  = bv->orthog_block;
1031:   return(0);
1032: }

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

1037:    Logically Collective on BV

1039:    Input Parameters:
1040: +  bv     - the basis vectors context
1041: -  method - the method for the BVMatMult() operation

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

1046:    Notes:
1047:    Allowed values are
1048: +  BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1049: .  BV_MATMULT_MAT - carry out a MatMatMult() product with a dense matrix
1050: -  BV_MATMULT_MAT_SAVE - this case is deprecated

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

1054:    Level: advanced

1056: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1057: @*/
1058: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1059: {

1065:   switch (method) {
1066:     case BV_MATMULT_VECS:
1067:     case BV_MATMULT_MAT:
1068:       bv->vmm = method;
1069:       break;
1070:     case BV_MATMULT_MAT_SAVE:
1071:       PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n");
1072:       bv->vmm = BV_MATMULT_MAT;
1073:       break;
1074:     default:
1075:       SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1076:   }
1077:   return(0);
1078: }

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

1083:    Not Collective

1085:    Input Parameter:
1086: .  bv - basis vectors context

1088:    Output Parameter:
1089: .  method - the method for the BVMatMult() operation

1091:    Level: advanced

1093: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1094: @*/
1095: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1096: {
1100:   *method = bv->vmm;
1101:   return(0);
1102: }

1104: /*@
1105:    BVGetColumn - Returns a Vec object that contains the entries of the
1106:    requested column of the basis vectors object.

1108:    Logically Collective on BV

1110:    Input Parameters:
1111: +  bv - the basis vectors context
1112: -  j  - the index of the requested column

1114:    Output Parameter:
1115: .  v  - vector containing the jth column

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

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

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

1129:    Level: beginner

1131: .seealso: BVRestoreColumn(), BVInsertConstraints()
1132: @*/
1133: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1134: {
1136:   PetscInt       l;

1141:   BVCheckSizes(bv,1);
1143:   if (j<0 && -j>bv->nc) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %D but only %D are available",-j,bv->nc);
1144:   if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %D but only %D are available",j,bv->m);
1145:   if (j==bv->ci[0] || j==bv->ci[1]) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Column %D already fetched in a previous call to BVGetColumn",j);
1146:   l = BVAvailableVec;
1147:   if (l==-1) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1148:   (*bv->ops->getcolumn)(bv,j,v);
1149:   bv->ci[l] = j;
1150:   PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]);
1151:   PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]);
1152:   *v = bv->cv[l];
1153:   return(0);
1154: }

1156: /*@
1157:    BVRestoreColumn - Restore a column obtained with BVGetColumn().

1159:    Logically Collective on BV

1161:    Input Parameters:
1162: +  bv - the basis vectors context
1163: .  j  - the index of the column
1164: -  v  - vector obtained with BVGetColumn()

1166:    Note:
1167:    The arguments must match the corresponding call to BVGetColumn().

1169:    Level: beginner

1171: .seealso: BVGetColumn()
1172: @*/
1173: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1174: {
1175:   PetscErrorCode   ierr;
1176:   PetscObjectId    id;
1177:   PetscObjectState st;
1178:   PetscInt         l;

1183:   BVCheckSizes(bv,1);
1187:   if (j<0 && -j>bv->nc) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %D but only %D are available",-j,bv->nc);
1188:   if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %D but only %D are available",j,bv->m);
1189:   if (j!=bv->ci[0] && j!=bv->ci[1]) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Column %D has not been fetched with a call to BVGetColumn",j);
1190:   l = (j==bv->ci[0])? 0: 1;
1191:   PetscObjectGetId((PetscObject)*v,&id);
1192:   if (id!=bv->id[l]) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1193:   PetscObjectStateGet((PetscObject)*v,&st);
1194:   if (st!=bv->st[l]) {
1195:     PetscObjectStateIncrease((PetscObject)bv);
1196:   }
1197:   if (bv->ops->restorecolumn) {
1198:     (*bv->ops->restorecolumn)(bv,j,v);
1199:   } else bv->cv[l] = NULL;
1200:   bv->ci[l] = -bv->nc-1;
1201:   bv->st[l] = -1;
1202:   bv->id[l] = 0;
1203:   *v = NULL;
1204:   return(0);
1205: }

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

1211:    Logically Collective on BV

1213:    Input Parameters:
1214: .  bv - the basis vectors context

1216:    Output Parameter:
1217: .  a  - location to put pointer to the array

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

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

1227:    Level: advanced

1229: .seealso: BVRestoreArray(), BVInsertConstraints()
1230: @*/
1231: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1232: {

1238:   BVCheckSizes(bv,1);
1239:   (*bv->ops->getarray)(bv,a);
1240:   return(0);
1241: }

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

1246:    Logically Collective on BV

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

1252:    Note:
1253:    This operation may imply a data copy, for BV types that do not store
1254:    data contiguously in memory.

1256:    Level: advanced

1258: .seealso: BVGetColumn()
1259: @*/
1260: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1261: {

1267:   BVCheckSizes(bv,1);
1268:   if (bv->ops->restorearray) {
1269:     (*bv->ops->restorearray)(bv,a);
1270:   }
1271:   if (a) *a = NULL;
1272:   PetscObjectStateIncrease((PetscObject)bv);
1273:   return(0);
1274: }

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

1280:    Not Collective

1282:    Input Parameters:
1283: .  bv - the basis vectors context

1285:    Output Parameter:
1286: .  a  - location to put pointer to the array

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

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

1296:    Level: advanced

1298: .seealso: BVRestoreArray(), BVInsertConstraints()
1299: @*/
1300: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1301: {

1307:   BVCheckSizes(bv,1);
1308:   (*bv->ops->getarrayread)(bv,a);
1309:   return(0);
1310: }

1312: /*@C
1313:    BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1314:    been called.

1316:    Not Collective

1318:    Input Parameters:
1319: +  bv - the basis vectors context
1320: -  a  - location of pointer to array obtained from BVGetArrayRead()

1322:    Level: advanced

1324: .seealso: BVGetColumn()
1325: @*/
1326: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1327: {

1333:   BVCheckSizes(bv,1);
1334:   if (bv->ops->restorearrayread) {
1335:     (*bv->ops->restorearrayread)(bv,a);
1336:   }
1337:   if (a) *a = NULL;
1338:   return(0);
1339: }

1341: /*@
1342:    BVCreateVec - Creates a new Vec object with the same type and dimensions
1343:    as the columns of the basis vectors object.

1345:    Collective on BV

1347:    Input Parameter:
1348: .  bv - the basis vectors context

1350:    Output Parameter:
1351: .  v  - the new vector

1353:    Note:
1354:    The user is responsible of destroying the returned vector.

1356:    Level: beginner

1358: .seealso: BVCreateMat()
1359: @*/
1360: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1361: {

1366:   BVCheckSizes(bv,1);
1368:   VecDuplicate(bv->t,v);
1369:   return(0);
1370: }

1372: /*@
1373:    BVCreateMat - Creates a new Mat object of dense type and copies the contents
1374:    of the BV object.

1376:    Collective on BV

1378:    Input Parameter:
1379: .  bv - the basis vectors context

1381:    Output Parameter:
1382: .  A  - the new matrix

1384:    Notes:
1385:    The user is responsible of destroying the returned matrix.

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

1389:    Level: intermediate

1391: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1392: @*/
1393: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1394: {
1395:   PetscErrorCode    ierr;
1396:   PetscScalar       *aa;
1397:   const PetscScalar *vv;

1401:   BVCheckSizes(bv,1);

1404:   MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,bv->N,bv->m,NULL,A);
1405:   MatDenseGetArray(*A,&aa);
1406:   BVGetArrayRead(bv,&vv);
1407:   PetscMemcpy(aa,vv,bv->m*bv->n*sizeof(PetscScalar));
1408:   BVRestoreArrayRead(bv,&vv);
1409:   MatDenseRestoreArray(*A,&aa);
1410:   MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1411:   MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1412:   return(0);
1413: }

1415: /*@
1416:    BVGetMat - Returns a Mat object of dense type that shares the memory of
1417:    the BV object.

1419:    Collective on BV

1421:    Input Parameter:
1422: .  bv - the basis vectors context

1424:    Output Parameter:
1425: .  A  - the matrix

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

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

1435:    Level: advanced

1437: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1438: @*/
1439: PetscErrorCode BVGetMat(BV bv,Mat *A)
1440: {
1442:   PetscScalar    *vv,*aa;
1443:   PetscBool      create=PETSC_FALSE;
1444:   PetscInt       m,cols;

1448:   BVCheckSizes(bv,1);
1450:   m = bv->k-bv->l;
1451:   if (!bv->Aget) create=PETSC_TRUE;
1452:   else {
1453:     MatDenseGetArray(bv->Aget,&aa);
1454:     if (aa) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1455:     MatGetSize(bv->Aget,NULL,&cols);
1456:     if (cols!=m) {
1457:       MatDestroy(&bv->Aget);
1458:       create=PETSC_TRUE;
1459:     }
1460:   }
1461:   BVGetArray(bv,&vv);
1462:   if (create) {
1463:     MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,vv,&bv->Aget); /* pass a pointer to avoid allocation of storage */
1464:     MatDensePlaceArray(bv->Aget,NULL);  /* replace with a null pointer, the value after BVRestoreMat */
1465:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Aget);
1466:     MatAssemblyBegin(bv->Aget,MAT_FINAL_ASSEMBLY);
1467:     MatAssemblyEnd(bv->Aget,MAT_FINAL_ASSEMBLY);
1468:   }
1469:   MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->n);  /* set the actual pointer */
1470:   *A = bv->Aget;
1471:   return(0);
1472: }

1474: /*@
1475:    BVRestoreMat - Restores the Mat obtained with BVGetMat().

1477:    Collective on BV

1479:    Input Parameters:
1480: +  bv - the basis vectors context
1481: -  A  - the fetched matrix

1483:    Note:
1484:    A call to this function must match a previous call of BVGetMat().
1485:    The effect is that the contents of the Mat are copied back to the
1486:    BV internal data structures.

1488:    Level: advanced

1490: .seealso: BVGetMat(), BVRestoreArray()
1491: @*/
1492: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1493: {
1495:   PetscScalar    *vv,*aa;

1499:   BVCheckSizes(bv,1);
1501:   if (!bv->Aget) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1502:   if (bv->Aget!=*A) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1503:   MatDenseGetArray(bv->Aget,&aa);
1504:   vv = aa-(bv->nc+bv->l)*bv->n;
1505:   MatDenseResetArray(bv->Aget);
1506:   BVRestoreArray(bv,&vv);
1507:   *A = NULL;
1508:   return(0);
1509: }

1511: /*
1512:    Copy all user-provided attributes of V to another BV object W
1513:  */
1514: PETSC_STATIC_INLINE PetscErrorCode BVDuplicate_Private(BV V,BV W)
1515: {

1519:   BVSetType(W,((PetscObject)V)->type_name);
1520:   W->orthog_type  = V->orthog_type;
1521:   W->orthog_ref   = V->orthog_ref;
1522:   W->orthog_eta   = V->orthog_eta;
1523:   W->orthog_block = V->orthog_block;
1524:   if (V->matrix) { PetscObjectReference((PetscObject)V->matrix); }
1525:   W->matrix       = V->matrix;
1526:   W->indef        = V->indef;
1527:   W->vmm          = V->vmm;
1528:   W->rrandom      = V->rrandom;
1529:   if (V->ops->duplicate) { (*V->ops->duplicate)(V,W); }
1530:   PetscObjectStateIncrease((PetscObject)W);
1531:   return(0);
1532: }

1534: /*@
1535:    BVDuplicate - Creates a new basis vector object of the same type and
1536:    dimensions as an existing one.

1538:    Collective on BV

1540:    Input Parameter:
1541: .  V - basis vectors context

1543:    Output Parameter:
1544: .  W - location to put the new BV

1546:    Notes:
1547:    The new BV has the same type and dimensions as V, and it shares the same
1548:    template vector. Also, the inner product matrix and orthogonalization
1549:    options are copied.

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

1554:    Level: intermediate

1556: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1557: @*/
1558: PetscErrorCode BVDuplicate(BV V,BV *W)
1559: {

1565:   BVCheckSizes(V,1);
1567:   BVCreate(PetscObjectComm((PetscObject)V),W);
1568:   BVSetSizesFromVec(*W,V->t,V->m);
1569:   BVDuplicate_Private(V,*W);
1570:   return(0);
1571: }

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

1577:    Collective on BV

1579:    Input Parameter:
1580: +  V - basis vectors context
1581: -  m - the new number of columns

1583:    Output Parameter:
1584: .  W - location to put the new BV

1586:    Note:
1587:    This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1588:    contents of V are not copied to W.

1590:    Level: intermediate

1592: .seealso: BVDuplicate(), BVResize()
1593: @*/
1594: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1595: {

1601:   BVCheckSizes(V,1);
1604:   BVCreate(PetscObjectComm((PetscObject)V),W);
1605:   BVSetSizesFromVec(*W,V->t,m);
1606:   BVDuplicate_Private(V,*W);
1607:   return(0);
1608: }

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

1614:    Not collective

1616:    Input Parameter:
1617: .  bv    - the basis vectors context

1619:    Output Parameter:
1620: .  cached - the cached BV

1622:    Note:
1623:    The cached BV is created if not available previously.

1625:    Level: developer

1627: .seealso: BVApplyMatrixBV()
1628: @*/
1629: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1630: {

1636:   BVCheckSizes(bv,1);
1637:   if (!bv->cached) {
1638:     BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached);
1639:     BVSetSizesFromVec(bv->cached,bv->t,bv->m);
1640:     BVDuplicate_Private(bv,bv->cached);
1641:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->cached);
1642:   }
1643:   *cached = bv->cached;
1644:   return(0);
1645: }

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

1650:    Logically Collective on BV

1652:    Input Parameter:
1653: .  V - basis vectors context

1655:    Output Parameter:
1656: .  W - the copy

1658:    Note:
1659:    Both V and W must be distributed in the same manner; local copies are
1660:    done. Only active columns (excluding the leading ones) are copied.
1661:    In the destination W, columns are overwritten starting from the leading ones.
1662:    Constraints are not copied.

1664:    Level: beginner

1666: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1667: @*/
1668: PetscErrorCode BVCopy(BV V,BV W)
1669: {
1670:   PetscErrorCode    ierr;
1671:   PetscScalar       *womega;
1672:   const PetscScalar *vomega;

1677:   BVCheckSizes(V,1);
1680:   BVCheckSizes(W,2);
1682:   if (V->n!=W->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, W %D",V->n,W->n);
1683:   if (V->k-V->l>W->m-W->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"W has %D non-leading columns, not enough to store %D columns",W->m-W->l,V->k-V->l);
1684:   if (!V->n) return(0);

1686:   PetscLogEventBegin(BV_Copy,V,W,0,0);
1687:   if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1688:     /* copy signature */
1689:     BV_AllocateSignature(W);
1690:     VecGetArrayRead(V->omega,&vomega);
1691:     VecGetArray(W->omega,&womega);
1692:     PetscMemcpy(womega+W->nc+W->l,vomega+V->nc+V->l,(V->k-V->l)*sizeof(PetscScalar));
1693:     VecRestoreArray(W->omega,&womega);
1694:     VecRestoreArrayRead(V->omega,&vomega);
1695:   }
1696:   (*V->ops->copy)(V,W);
1697:   PetscLogEventEnd(BV_Copy,V,W,0,0);
1698:   PetscObjectStateIncrease((PetscObject)W);
1699:   return(0);
1700: }

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

1705:    Logically Collective on BV

1707:    Input Parameter:
1708: +  V - basis vectors context
1709: -  j - the column number to be copied

1711:    Output Parameter:
1712: .  w - the copied column

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

1717:    Level: beginner

1719: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1720: @*/
1721: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1722: {
1724:   PetscInt       n,N;
1725:   Vec            z;

1730:   BVCheckSizes(V,1);

1735:   VecGetSize(w,&N);
1736:   VecGetLocalSize(w,&n);
1737:   if (N!=V->N || n!=V->n) SETERRQ4(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %D, local %D) do not match BV sizes (global %D, local %D)",N,n,V->N,V->n);

1739:   PetscLogEventBegin(BV_Copy,V,w,0,0);
1740:   BVGetColumn(V,j,&z);
1741:   VecCopy(z,w);
1742:   BVRestoreColumn(V,j,&z);
1743:   PetscLogEventEnd(BV_Copy,V,w,0,0);
1744:   return(0);
1745: }

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

1750:    Logically Collective on BV

1752:    Input Parameter:
1753: +  V - basis vectors context
1754: .  j - the number of the source column
1755: -  i - the number of the destination column

1757:    Level: beginner

1759: .seealso: BVCopy(), BVCopyVec()
1760: @*/
1761: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1762: {
1764:   Vec            z,w;
1765:   PetscScalar    *omega;

1770:   BVCheckSizes(V,1);
1773:   if (j==i) return(0);

1775:   PetscLogEventBegin(BV_Copy,V,0,0,0);
1776:   if (V->omega) {
1777:     VecGetArray(V->omega,&omega);
1778:     omega[i] = omega[j];
1779:     VecRestoreArray(V->omega,&omega);
1780:   }
1781:   if (V->ops->copycolumn) {
1782:     (*V->ops->copycolumn)(V,j,i);
1783:   } else {
1784:     BVGetColumn(V,j,&z);
1785:     BVGetColumn(V,i,&w);
1786:     VecCopy(z,w);
1787:     BVRestoreColumn(V,j,&z);
1788:     BVRestoreColumn(V,i,&w);
1789:   }
1790:   PetscLogEventEnd(BV_Copy,V,0,0,0);
1791:   PetscObjectStateIncrease((PetscObject)V);
1792:   return(0);
1793: }

1795: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1796: {
1798:   PetscInt       ncols;

1801:   ncols = left? bv->nc+bv->l: bv->m-bv->l;
1802:   if (!*split) {
1803:     BVCreate(PetscObjectComm((PetscObject)bv),split);
1804:     PetscLogObjectParent((PetscObject)bv,(PetscObject)*split);
1805:     (*split)->issplit = left? 1: 2;
1806:     (*split)->splitparent = bv;
1807:     BVSetSizesFromVec(*split,bv->t,ncols);
1808:     BVDuplicate_Private(bv,*split);
1809:   } else {
1810:     BVResize(*split,ncols,PETSC_FALSE);
1811:   }
1812:   (*split)->l  = 0;
1813:   (*split)->k  = left? bv->l: bv->k-bv->l;
1814:   (*split)->nc = left? bv->nc: 0;
1815:   (*split)->m  = ncols-(*split)->nc;
1816:   if ((*split)->nc) {
1817:     (*split)->ci[0] = -(*split)->nc-1;
1818:     (*split)->ci[1] = -(*split)->nc-1;
1819:   }
1820:   if (left) {
1821:     PetscObjectStateGet((PetscObject)*split,&bv->lstate);
1822:   } else {
1823:     PetscObjectStateGet((PetscObject)*split,&bv->rstate);
1824:   }
1825:   return(0);
1826: }

1828: /*@
1829:    BVGetSplit - Splits the BV object into two BV objects that share the
1830:    internal data, one of them containing the leading columns and the other
1831:    one containing the remaining columns.

1833:    Logically Collective on BV

1835:    Input Parameters:
1836: .  bv - the basis vectors context

1838:    Output Parameters:
1839: +  L - left BV containing leading columns (can be NULL)
1840: -  R - right BV containing remaining columns (can be NULL)

1842:    Notes:
1843:    The columns are split in two sets. The leading columns (including the
1844:    constraints) are assigned to the left BV and the remaining columns
1845:    are assigned to the right BV. The number of leading columns, as
1846:    specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1847:    guarantee that both L and R have at least one column).

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

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

1856:    Level: advanced

1858: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1859: @*/
1860: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1861: {

1867:   BVCheckSizes(bv,1);
1868:   if (!bv->l) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1869:   if (bv->lsplit) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1870:   bv->lsplit = bv->nc+bv->l;
1871:   BVGetSplit_Private(bv,PETSC_TRUE,&bv->L);
1872:   BVGetSplit_Private(bv,PETSC_FALSE,&bv->R);
1873:   if (L) *L = bv->L;
1874:   if (R) *R = bv->R;
1875:   return(0);
1876: }

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

1881:    Logically Collective on BV

1883:    Input Parameters:
1884: +  bv - the basis vectors context
1885: .  L  - left BV obtained with BVGetSplit()
1886: -  R  - right BV obtained with BVGetSplit()

1888:    Note:
1889:    The arguments must match the corresponding call to BVGetSplit().

1891:    Level: advanced

1893: .seealso: BVGetSplit()
1894: @*/
1895: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1896: {

1902:   BVCheckSizes(bv,1);
1903:   if (!bv->lsplit) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1904:   if (L && *L!=bv->L) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1905:   if (R && *R!=bv->R) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
1906:   if (L && ((*L)->ci[0]>(*L)->nc-1 || (*L)->ci[1]>(*L)->nc-1)) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 2 has unrestored columns, use BVRestoreColumn()");
1907:   if (R && ((*R)->ci[0]>(*R)->nc-1 || (*R)->ci[1]>(*R)->nc-1)) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 3 has unrestored columns, use BVRestoreColumn()");

1909:   if (bv->ops->restoresplit) {
1910:     (*bv->ops->restoresplit)(bv,L,R);
1911:   }
1912:   bv->lsplit = 0;
1913:   if (L) *L = NULL;
1914:   if (R) *R = NULL;
1915:   return(0);
1916: }