Actual source code: bvbasic.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Basic BV routines
12: */
14: #include <slepc/private/bvimpl.h>
16: PetscBool BVRegisterAllCalled = PETSC_FALSE;
17: PetscFunctionList BVList = NULL;
19: /*@
20: BVSetType - Selects the type for the `BV` object.
22: Logically Collective
24: Input Parameters:
25: + bv - the basis vectors context
26: - type - a known type
28: Options Database Key:
29: . -bv_type <type> - Sets `BV` type
31: Level: intermediate
33: .seealso: [](sec:bv), `BVGetType()`
34: @*/
35: PetscErrorCode BVSetType(BV bv,BVType type)
36: {
37: PetscErrorCode (*r)(BV);
38: PetscBool match;
40: PetscFunctionBegin;
42: PetscAssertPointer(type,2);
44: PetscCall(PetscObjectTypeCompare((PetscObject)bv,type,&match));
45: if (match) PetscFunctionReturn(PETSC_SUCCESS);
46: PetscCall(PetscStrcmp(type,BVTENSOR,&match));
47: PetscCheck(!match,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Use BVCreateTensor() to create a BV of type tensor");
49: PetscCall(PetscFunctionListFind(BVList,type,&r));
50: PetscCheck(r,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested BV type %s",type);
52: PetscTryTypeMethod(bv,destroy);
53: PetscCall(PetscMemzero(bv->ops,sizeof(struct _BVOps)));
55: PetscCall(PetscObjectChangeTypeName((PetscObject)bv,type));
56: if (bv->n < 0 && bv->N < 0) {
57: bv->ops->create = r;
58: } else {
59: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
60: PetscCall((*r)(bv));
61: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
62: }
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: /*@
67: BVGetType - Gets the `BV` type name (as a string) from the `BV` context.
69: Not Collective
71: Input Parameter:
72: . bv - the basis vectors context
74: Output Parameter:
75: . type - name of the type of basis vectors
77: Level: intermediate
79: .seealso: [](sec:bv), `BVSetType()`
80: @*/
81: PetscErrorCode BVGetType(BV bv,BVType *type)
82: {
83: PetscFunctionBegin;
85: PetscAssertPointer(type,2);
86: *type = ((PetscObject)bv)->type_name;
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: /*@
91: BVSetSizes - Sets the local and global sizes, and the number of columns.
93: Collective
95: Input Parameters:
96: + bv - the basis vectors
97: . n - the local size (or `PETSC_DECIDE` to have it set)
98: . N - the global size (or `PETSC_DECIDE`)
99: - m - the number of columns
101: Notes:
102: `n` and `N` cannot be both `PETSC_DECIDE`.
103: If one process calls this with `N` equal to `PETSC_DECIDE` then all processes must,
104: otherwise the program will hang.
106: Level: beginner
108: .seealso: [](sec:bv), `BVSetSizesFromVec()`, `BVGetSizes()`, `BVResize()`
109: @*/
110: PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)
111: {
112: PetscInt ma;
113: PetscMPIInt size;
115: PetscFunctionBegin;
119: PetscCheck(N<0 || n<=N,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT,n,N);
120: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
121: PetscCheck((bv->n<0 && bv->N<0) || (bv->n==n && bv->N==N),PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change/reset vector sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global",n,N,bv->n,bv->N);
122: PetscCheck(bv->m<=0 || bv->m==m,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change the number of columns to %" PetscInt_FMT " after previously setting it to %" PetscInt_FMT "; use BVResize()",m,bv->m);
123: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
124: bv->n = n;
125: bv->N = N;
126: bv->m = m;
127: bv->k = m;
128: /* create layout and get actual dimensions */
129: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)bv),&bv->map));
130: PetscCall(PetscLayoutSetSize(bv->map,bv->N));
131: PetscCall(PetscLayoutSetLocalSize(bv->map,bv->n));
132: PetscCall(PetscLayoutSetUp(bv->map));
133: PetscCall(PetscLayoutGetSize(bv->map,&bv->N));
134: PetscCall(PetscLayoutGetLocalSize(bv->map,&bv->n));
135: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
136: PetscCall(BVSetVecType(bv,(size==1)?VECSEQ:VECMPI));
137: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
138: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
139: PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
140: }
141: if (bv->ops->create) {
142: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
143: PetscUseTypeMethod(bv,create);
144: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
145: bv->ops->create = NULL;
146: bv->defersfo = PETSC_FALSE;
147: }
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
153: Local and global sizes are specified indirectly by passing a template vector.
155: Collective
157: Input Parameters:
158: + bv - the basis vectors
159: . t - the template vector
160: - m - the number of columns
162: Level: beginner
164: .seealso: [](sec:bv), `BVSetSizes()`, `BVGetSizes()`, `BVResize()`
165: @*/
166: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
167: {
168: PetscInt ma;
169: PetscLayout map;
170: VecType vtype;
172: PetscFunctionBegin;
175: PetscCheckSameComm(bv,1,t,2);
177: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
178: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
179: PetscCall(VecGetType(t,&vtype));
180: PetscCall(BVSetVecType(bv,vtype));
181: PetscCall(VecGetLayout(t,&map));
182: PetscCall(PetscLayoutReference(map,&bv->map));
183: PetscCall(VecGetSize(t,&bv->N));
184: PetscCall(VecGetLocalSize(t,&bv->n));
185: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
186: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
187: PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
188: }
189: bv->m = m;
190: bv->k = m;
191: if (bv->ops->create) {
192: PetscUseTypeMethod(bv,create);
193: bv->ops->create = NULL;
194: bv->defersfo = PETSC_FALSE;
195: }
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: BVGetSizes - Returns the local and global sizes, and the number of columns.
202: Not Collective
204: Input Parameter:
205: . bv - the basis vectors
207: Output Parameters:
208: + n - the local size
209: . N - the global size
210: - m - the number of columns
212: Note:
213: Normal usage requires that `bv` has already been given its sizes, otherwise
214: the call fails. However, this function can also be used to determine if
215: a `BV` object has been initialized completely (sizes and type). For this,
216: call with `n=NULL` and `N=NULL`, then a return value of `m=0` indicates that
217: the `BV` object is not ready for use yet.
219: Level: beginner
221: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
222: @*/
223: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
224: {
225: PetscFunctionBegin;
226: if (!bv) {
227: if (m && !n && !N) *m = 0;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: if (n || N) BVCheckSizes(bv,1);
232: if (n) *n = bv->n;
233: if (N) *N = bv->N;
234: if (m) *m = bv->m;
235: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: BVSetNumConstraints - Set the number of constraints.
242: Logically Collective
244: Input Parameters:
245: + V - basis vectors
246: - nc - number of constraints
248: Notes:
249: This function sets the number of constraints to `nc` and marks all remaining
250: columns as regular. Normal user would call `BVInsertConstraints()` instead.
252: If `nc` is smaller than the previously set value, then some of the constraints
253: are discarded. In particular, using `nc=0` removes all constraints preserving
254: the content of regular columns.
256: Level: developer
258: .seealso: [](sec:bv), `BVInsertConstraints()`
259: @*/
260: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
261: {
262: PetscInt total,diff,i;
263: Vec x,y;
265: PetscFunctionBegin;
268: PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
270: BVCheckSizes(V,1);
271: PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
273: diff = nc-V->nc;
274: if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
275: total = V->nc+V->m;
276: PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
277: if (diff<0) { /* lessen constraints, shift contents of BV */
278: for (i=0;i<V->m;i++) {
279: PetscCall(BVGetColumn(V,i,&x));
280: PetscCall(BVGetColumn(V,i+diff,&y));
281: PetscCall(VecCopy(x,y));
282: PetscCall(BVRestoreColumn(V,i,&x));
283: PetscCall(BVRestoreColumn(V,i+diff,&y));
284: }
285: }
286: V->nc = nc;
287: V->ci[0] = -V->nc-1;
288: V->ci[1] = -V->nc-1;
289: V->m = total-nc;
290: V->l = PetscMin(V->l,V->m);
291: V->k = PetscMin(V->k,V->m);
292: PetscCall(PetscObjectStateIncrease((PetscObject)V));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: BVGetNumConstraints - Returns the number of constraints.
299: Not Collective
301: Input Parameter:
302: . bv - the basis vectors
304: Output Parameter:
305: . nc - the number of constraints
307: Level: advanced
309: .seealso: [](sec:bv), `BVGetSizes()`, `BVInsertConstraints()`
310: @*/
311: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
312: {
313: PetscFunctionBegin;
315: PetscAssertPointer(nc,2);
316: *nc = bv->nc;
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@
321: BVResize - Change the number of columns.
323: Collective
325: Input Parameters:
326: + bv - the basis vectors
327: . m - the new number of columns
328: - copy - a flag indicating whether current values should be kept
330: Note:
331: Internal storage is reallocated. If the `copy` flag is set to true, then
332: the contents are copied to the leading part of the new space.
334: Level: advanced
336: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
337: @*/
338: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
339: {
340: PetscScalar *array;
341: const PetscScalar *omega;
342: Vec v;
344: PetscFunctionBegin;
349: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
350: PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
351: if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
352: BVCheckOp(bv,1,resize);
354: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
355: PetscUseTypeMethod(bv,resize,m,copy);
356: PetscCall(VecDestroy(&bv->buffer));
357: PetscCall(BVDestroy(&bv->cached));
358: PetscCall(PetscFree2(bv->h,bv->c));
359: if (bv->omega) {
360: if (bv->cuda) {
361: #if defined(PETSC_HAVE_CUDA)
362: PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
363: #else
364: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
365: #endif
366: } else if (bv->hip) {
367: #if defined(PETSC_HAVE_HIP)
368: PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
369: #else
370: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
371: #endif
372: } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
373: if (copy) {
374: PetscCall(VecGetArray(v,&array));
375: PetscCall(VecGetArrayRead(bv->omega,&omega));
376: PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
377: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
378: PetscCall(VecRestoreArray(v,&array));
379: } else PetscCall(VecSet(v,1.0));
380: PetscCall(VecDestroy(&bv->omega));
381: bv->omega = v;
382: }
383: bv->m = m;
384: bv->k = PetscMin(bv->k,m);
385: bv->l = PetscMin(bv->l,m);
386: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
387: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: BVSetActiveColumns - Specify the columns that will be involved in operations.
394: Logically Collective
396: Input Parameters:
397: + bv - the basis vectors context
398: . l - number of leading columns
399: - k - number of active columns
401: Notes:
402: In operations such as `BVMult()` or `BVDot()`, only the first `k` columns are
403: considered. This is useful when the `BV` is filled from left to right, so
404: the last `m-k` columns do not have relevant information.
406: Also in operations such as `BVMult()` or `BVDot()`, the first `l` columns are
407: normally not included in the computation. See the manual page of each
408: operation.
410: In orthogonalization operations, the first `l` columns are treated
411: differently, they participate in the orthogonalization but the computed
412: coefficients are not stored.
414: Use `PETSC_CURRENT` to leave any of the values unchanged. Use `PETSC_DETERMINE`
415: to set `l` to the minimum value (`0`) and `k` to the maximum (`m`).
417: Level: intermediate
419: .seealso: [](sec:bv), `BVGetActiveColumns()`, `BVSetSizes()`
420: @*/
421: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
422: {
423: PetscFunctionBegin;
427: BVCheckSizes(bv,1);
428: if (PetscUnlikely(k == PETSC_DETERMINE)) {
429: bv->k = bv->m;
430: } else if (k != PETSC_CURRENT) {
431: PetscCheck(k>=0 && k<=bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k (%" PetscInt_FMT "). Must be between 0 and m (%" PetscInt_FMT ")",k,bv->m);
432: bv->k = k;
433: }
434: if (PetscUnlikely(l == PETSC_DETERMINE)) {
435: bv->l = 0;
436: } else if (l != PETSC_CURRENT) {
437: PetscCheck(l>=0 && l<=bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l (%" PetscInt_FMT "). Must be between 0 and k (%" PetscInt_FMT ")",l,bv->k);
438: bv->l = l;
439: }
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@
444: BVGetActiveColumns - Returns the current active dimensions.
446: Not Collective
448: Input Parameter:
449: . bv - the basis vectors context
451: Output Parameters:
452: + l - number of leading columns
453: - k - number of active columns
455: Level: intermediate
457: .seealso: [](sec:bv), `BVSetActiveColumns()`
458: @*/
459: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
460: {
461: PetscFunctionBegin;
463: if (l) *l = bv->l;
464: if (k) *k = bv->k;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
471: Collective
473: Input Parameters:
474: + bv - the basis vectors context
475: . B - a Hermitian matrix (may be `NULL`)
476: - indef - a flag indicating if the matrix is indefinite
478: Notes:
479: This is used to specify a non-standard inner product, whose matrix
480: representation is given by `B`. Then, all inner products required during
481: orthogonalization are computed as $(x,y)_B=y^*Bx$ rather than the
482: standard form $(x,y)=y^Hx$.
484: Matrix `B` must be real symmetric (or complex Hermitian). A genuine inner
485: product requires that `B` is also positive (semi-)definite. However, we
486: also allow for an indefinite `B` (setting `indef=PETSC_TRUE`), in which
487: case the orthogonalization uses an indefinite inner product.
489: This affects operations `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, and variants.
491: Setting `B=NULL` has the same effect as if the identity matrix was passed.
493: Level: advanced
495: .seealso: [](sec:bv), `BVGetMatrix()`, `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, `BVSetDefiniteTolerance()`
496: @*/
497: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
498: {
499: PetscInt m,n;
501: PetscFunctionBegin;
504: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
505: if (B) {
507: PetscCall(MatGetLocalSize(B,&m,&n));
508: PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
509: PetscCheck(!bv->m || bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension BV %" PetscInt_FMT ", Mat %" PetscInt_FMT,bv->n,n);
510: }
511: if (B) PetscCall(PetscObjectReference((PetscObject)B));
512: PetscCall(MatDestroy(&bv->matrix));
513: bv->matrix = B;
514: bv->indef = indef;
515: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
516: if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
517: if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
518: }
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: BVGetMatrix - Retrieves the matrix representation of the inner product.
525: Not Collective
527: Input Parameter:
528: . bv - the basis vectors context
530: Output Parameters:
531: + B - the matrix of the inner product (may be `NULL`)
532: - indef - the flag indicating if the matrix is indefinite
534: Level: advanced
536: .seealso: [](sec:bv), `BVSetMatrix()`
537: @*/
538: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
539: {
540: PetscFunctionBegin;
542: if (B) *B = bv->matrix;
543: if (indef) *indef = bv->indef;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: BVApplyMatrix - Multiplies a vector by the matrix representation of the
549: inner product.
551: Neighbor-wise Collective
553: Input Parameters:
554: + bv - the basis vectors context
555: - x - the vector
557: Output Parameter:
558: . y - the result
560: Note:
561: If no matrix was specified this function copies the vector.
563: Level: advanced
565: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrixBV()`
566: @*/
567: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
568: {
569: PetscFunctionBegin;
573: if (bv->matrix) {
574: PetscCall(BV_IPMatMult(bv,x));
575: PetscCall(VecCopy(bv->Bx,y));
576: } else PetscCall(VecCopy(x,y));
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: BVApplyMatrixBV - Multiplies the `BV` vectors by the matrix representation
582: of the inner product.
584: Neighbor-wise Collective
586: Input Parameter:
587: . X - the basis vectors context
589: Output Parameter:
590: . Y - the basis vectors to store the result (optional)
592: Note:
593: This function computes $Y = B X$, where $B$ is the matrix given with
594: `BVSetMatrix()`. This operation is computed as in `BVMatMult()`.
595: If no matrix was specified, then it just copies $Y = X$.
597: If no `Y` is given, the result is stored internally in the cached `BV`.
599: Level: developer
601: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrix()`, `BVMatMult()`, `BVGetCachedBV()`
602: @*/
603: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
604: {
605: PetscFunctionBegin;
607: if (Y) {
609: if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
610: else PetscCall(BVCopy(X,Y));
611: } else PetscCall(BV_IPMatMultBV(X));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@
616: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
618: Logically Collective
620: Input Parameters:
621: + bv - the basis vectors context
622: - omega - a vector representing the diagonal of the signature matrix
624: Note:
625: The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.
627: Level: developer
629: .seealso: [](sec:bv), `BVSetMatrix()`, `BVGetSignature()`
630: @*/
631: PetscErrorCode BVSetSignature(BV bv,Vec omega)
632: {
633: PetscInt i,n;
634: const PetscScalar *pomega;
635: PetscScalar *intern;
637: PetscFunctionBegin;
639: BVCheckSizes(bv,1);
643: PetscCall(VecGetSize(omega,&n));
644: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
645: PetscCall(BV_AllocateSignature(bv));
646: if (bv->indef) {
647: PetscCall(VecGetArrayRead(omega,&pomega));
648: PetscCall(VecGetArray(bv->omega,&intern));
649: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
650: PetscCall(VecRestoreArray(bv->omega,&intern));
651: PetscCall(VecRestoreArrayRead(omega,&pomega));
652: } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
653: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: /*@
658: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
660: Not Collective
662: Input Parameter:
663: . bv - the basis vectors context
665: Output Parameter:
666: . omega - a vector representing the diagonal of the signature matrix
668: Note:
669: The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.
671: Level: developer
673: .seealso: [](sec:bv), `BVSetMatrix()`, `BVSetSignature()`
674: @*/
675: PetscErrorCode BVGetSignature(BV bv,Vec omega)
676: {
677: PetscInt i,n;
678: PetscScalar *pomega;
679: const PetscScalar *intern;
681: PetscFunctionBegin;
683: BVCheckSizes(bv,1);
687: PetscCall(VecGetSize(omega,&n));
688: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
689: if (bv->indef && bv->omega) {
690: PetscCall(VecGetArray(omega,&pomega));
691: PetscCall(VecGetArrayRead(bv->omega,&intern));
692: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
693: PetscCall(VecRestoreArrayRead(bv->omega,&intern));
694: PetscCall(VecRestoreArray(omega,&pomega));
695: } else PetscCall(VecSet(omega,1.0));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: BVSetBufferVec - Attach a vector object to be used as buffer space for
701: several operations.
703: Collective
705: Input Parameters:
706: + bv - the basis vectors context)
707: - buffer - the vector
709: Notes:
710: Use `BVGetBufferVec()` to retrieve the vector (for example, to free it
711: at the end of the computations).
713: The vector must be sequential of length $(n_c+m)m$, where $m$ is the number
714: of columns of `bv` and $n_c$ is the number of constraints.
716: Level: developer
718: .seealso: [](sec:bv), `BVGetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`
719: @*/
720: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
721: {
722: PetscInt ld,n;
723: PetscMPIInt size;
725: PetscFunctionBegin;
728: BVCheckSizes(bv,1);
729: PetscCall(VecGetSize(buffer,&n));
730: ld = bv->m+bv->nc;
731: PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
732: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
733: PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
735: PetscCall(PetscObjectReference((PetscObject)buffer));
736: PetscCall(VecDestroy(&bv->buffer));
737: bv->buffer = buffer;
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: /*@
742: BVGetBufferVec - Obtain the buffer vector associated with the `BV` object.
744: Collective
746: Input Parameter:
747: . bv - the basis vectors context
749: Output Parameter:
750: . buffer - vector
752: Notes:
753: The vector is created if not available previously. It is a sequential vector
754: of length $(n_c+m) m$, where $m$ is the number of columns of `bv` and $n_c$ is the number
755: of constraints.
757: Developer Notes:
758: The buffer vector is viewed as a column-major matrix with leading dimension equal
759: to $n_c+m$, and $m$ columns at most. In the most common usage, it has the structure
760: $$
761: \left[\begin{array}{c|c}
762: s & \begin{array}{c} C \\\hline H \end{array}
763: \end{array}\right]
764: $$
765: where $H$ is an upper Hessenberg matrix of order $m (m-1)$, $C$ contains coefficients
766: related to orthogonalization against constraints (first $n_c$ rows), and $s$ is the
767: first column that contains scratch values computed during Gram-Schmidt
768: orthogonalization. In particular, `BVDotColumn()` and `BVMultColumn()` use $s$ to
769: store the coefficients.
771: Level: developer
773: .seealso: [](sec:bv), `BVSetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`, `BVDotColumn()`, `BVMultColumn()`
774: @*/
775: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
776: {
777: PetscInt ld;
779: PetscFunctionBegin;
781: PetscAssertPointer(buffer,2);
782: BVCheckSizes(bv,1);
783: if (!bv->buffer) {
784: ld = bv->m+bv->nc;
785: PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
786: PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
787: PetscCall(VecSetType(bv->buffer,bv->vtype));
788: }
789: *buffer = bv->buffer;
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: BVSetRandomContext - Sets the `PetscRandom` object associated with the `BV`,
795: to be used in operations that need random numbers.
797: Collective
799: Input Parameters:
800: + bv - the basis vectors context
801: - rand - the random number generator context
803: Level: advanced
805: .seealso: [](sec:bv), `BVGetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`, `PetscRandomCreate()`
806: @*/
807: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
808: {
809: PetscFunctionBegin;
812: PetscCheckSameComm(bv,1,rand,2);
813: PetscCall(PetscObjectReference((PetscObject)rand));
814: PetscCall(PetscRandomDestroy(&bv->rand));
815: bv->rand = rand;
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: /*@
820: BVGetRandomContext - Gets the `PetscRandom` object associated with the `BV`.
822: Collective
824: Input Parameter:
825: . bv - the basis vectors context
827: Output Parameter:
828: . rand - the random number generator context
830: Level: advanced
832: .seealso: [](sec:bv), `BVSetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`
833: @*/
834: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
835: {
836: PetscFunctionBegin;
838: PetscAssertPointer(rand,2);
839: if (!bv->rand) {
840: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
841: if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
842: if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
843: if (bv->rrandom) {
844: PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
845: PetscCall(PetscRandomSeed(bv->rand));
846: }
847: }
848: *rand = bv->rand;
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /*@
853: BVSetFromOptions - Sets `BV` options from the options database.
855: Collective
857: Input Parameter:
858: . bv - the basis vectors context
860: Note:
861: To see all options, run your program with the `-help` option.
863: Level: beginner
865: .seealso: [](sec:bv), `BVSetOptionsPrefix()`
866: @*/
867: PetscErrorCode BVSetFromOptions(BV bv)
868: {
869: char type[256];
870: PetscBool flg1,flg2,flg3,flg4;
871: PetscReal r;
872: BVOrthogType otype;
873: BVOrthogRefineType orefine;
874: BVOrthogBlockType oblock;
876: PetscFunctionBegin;
878: PetscCall(BVRegisterAll());
879: PetscObjectOptionsBegin((PetscObject)bv);
880: PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
881: if (flg1) PetscCall(BVSetType(bv,type));
882: else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));
884: otype = bv->orthog_type;
885: PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
886: orefine = bv->orthog_ref;
887: PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
888: r = bv->orthog_eta;
889: PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
890: oblock = bv->orthog_block;
891: PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
892: if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
894: PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
896: PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
897: if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
899: /* undocumented option to generate random vectors that are independent of the number of processes */
900: PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
902: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
903: else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
904: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
905: PetscOptionsEnd();
906: bv->sfocalled = PETSC_TRUE;
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@
911: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
912: vectors (classical or modified Gram-Schmidt with or without refinement), and
913: for the block-orthogonalization (simultaneous orthogonalization of a set of
914: vectors).
916: Logically Collective
918: Input Parameters:
919: + bv - the basis vectors context
920: . type - the method of vector orthogonalization
921: . refine - type of refinement
922: . eta - parameter for selective refinement
923: - block - the method of block orthogonalization
925: Options Database Keys:
926: + -bv_orthog_type <type> - the vector orthogonalization, either `cgs` or `mgs`
927: . -bv_orthog_refine <ref> - the refinement type, `never`, `ifneeded` (default) or `always`
928: . -bv_orthog_eta <eta> - the value of `eta`
929: - -bv_orthog_block <block> - the block-orthogonalization method
931: Notes:
932: The default settings work well for most problems.
934: The parameter `eta` should be a real value between 0 and 1, that is used only when
935: the refinement type is `ifneeded`. Use `PETSC_DETERMINE` to set a reasonable
936: default value. Use `PETSC_CURRENT` to leave the current value unchanged.
938: When using several processes, `BV_ORTHOG_MGS` is likely to result in bad scalability.
940: If the method set for block orthogonalization is `BV_ORTHOG_BLOCK_GS`, then the computation
941: is done column by column with the vector orthogonalization.
943: Level: advanced
945: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVGetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
946: @*/
947: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
948: {
949: PetscFunctionBegin;
955: switch (type) {
956: case BV_ORTHOG_CGS:
957: case BV_ORTHOG_MGS:
958: bv->orthog_type = type;
959: break;
960: default:
961: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
962: }
963: switch (refine) {
964: case BV_ORTHOG_REFINE_NEVER:
965: case BV_ORTHOG_REFINE_IFNEEDED:
966: case BV_ORTHOG_REFINE_ALWAYS:
967: bv->orthog_ref = refine;
968: break;
969: default:
970: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
971: }
972: if (eta == (PetscReal)PETSC_DETERMINE) {
973: bv->orthog_eta = 0.7071;
974: } else if (eta != (PetscReal)PETSC_CURRENT) {
975: PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
976: bv->orthog_eta = eta;
977: }
978: switch (block) {
979: case BV_ORTHOG_BLOCK_GS:
980: case BV_ORTHOG_BLOCK_CHOL:
981: case BV_ORTHOG_BLOCK_TSQR:
982: case BV_ORTHOG_BLOCK_TSQRCHOL:
983: case BV_ORTHOG_BLOCK_SVQB:
984: bv->orthog_block = block;
985: break;
986: default:
987: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
988: }
989: PetscFunctionReturn(PETSC_SUCCESS);
990: }
992: /*@
993: BVGetOrthogonalization - Gets the orthogonalization settings from the `BV` object.
995: Not Collective
997: Input Parameter:
998: . bv - basis vectors context
1000: Output Parameters:
1001: + type - the method of vector orthogonalization
1002: . refine - type of refinement
1003: . eta - parameter for selective refinement
1004: - block - the method of block orthogonalization
1006: Level: advanced
1008: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVSetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
1009: @*/
1010: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1011: {
1012: PetscFunctionBegin;
1014: if (type) *type = bv->orthog_type;
1015: if (refine) *refine = bv->orthog_ref;
1016: if (eta) *eta = bv->orthog_eta;
1017: if (block) *block = bv->orthog_block;
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: /*@
1022: BVSetMatMultMethod - Specifies the method used for the `BVMatMult()` operation.
1024: Logically Collective
1026: Input Parameters:
1027: + bv - the basis vectors context
1028: - method - the method for the `BVMatMult()` operation
1030: Options Database Key:
1031: . -bv_matmult <meth> - choose one of the methods: `vecs`, `mat`
1033: Note:
1034: The default is `BV_MATMULT_MAT` except in the case of `BVVECS`.
1036: Level: advanced
1038: .seealso: [](sec:bv), `BVMatMult()`, `BVGetMatMultMethod()`, `BVMatMultType`
1039: @*/
1040: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1041: {
1042: PetscFunctionBegin;
1045: switch (method) {
1046: case BV_MATMULT_VECS:
1047: case BV_MATMULT_MAT:
1048: bv->vmm = method;
1049: break;
1050: case BV_MATMULT_MAT_SAVE:
1051: PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1052: bv->vmm = BV_MATMULT_MAT;
1053: break;
1054: default:
1055: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1056: }
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: /*@
1061: BVGetMatMultMethod - Gets the method used for the `BVMatMult()` operation.
1063: Not Collective
1065: Input Parameter:
1066: . bv - basis vectors context
1068: Output Parameter:
1069: . method - the method for the `BVMatMult()` operation
1071: Level: advanced
1073: .seealso: [](sec:bv), `BVMatMult()`, `BVSetMatMultMethod()`, `BVMatMultType`
1074: @*/
1075: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1076: {
1077: PetscFunctionBegin;
1079: PetscAssertPointer(method,2);
1080: *method = bv->vmm;
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: /*@
1085: BVGetColumn - Returns a `Vec` object that contains the entries of the
1086: requested column of the basis vectors object.
1088: Logically Collective
1090: Input Parameters:
1091: + bv - the basis vectors context
1092: - j - the index of the requested column
1094: Output Parameter:
1095: . v - vector containing the jth column
1097: Notes:
1098: The returned `Vec` must be seen as a reference (not a copy) of the `BV`
1099: column, that is, modifying the `Vec` will change the `BV` entries as well.
1101: The returned `Vec` must not be destroyed. `BVRestoreColumn()` must be
1102: called when it is no longer needed. At most, two columns can be fetched,
1103: that is, this function can only be called twice before the corresponding
1104: `BVRestoreColumn()` is invoked.
1106: A negative index `j` selects the `i`-th constraint, where `i=-j`. Constraints
1107: should not be modified.
1109: Level: beginner
1111: .seealso: [](sec:bv), `BVRestoreColumn()`, `BVInsertConstraints()`
1112: @*/
1113: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1114: {
1115: PetscInt l;
1117: PetscFunctionBegin;
1120: BVCheckSizes(bv,1);
1121: BVCheckOp(bv,1,getcolumn);
1123: PetscCheck(j>=0 || -j<=bv->nc,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %" PetscInt_FMT " but only %" PetscInt_FMT " are available",-j,bv->nc);
1124: PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %" PetscInt_FMT " but only %" PetscInt_FMT " are available",j,bv->m);
1125: PetscCheck(j!=bv->ci[0] && j!=bv->ci[1],PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Column %" PetscInt_FMT " already fetched in a previous call to BVGetColumn",j);
1126: l = BVAvailableVec;
1127: PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1128: PetscUseTypeMethod(bv,getcolumn,j,v);
1129: bv->ci[l] = j;
1130: PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
1131: PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1132: *v = bv->cv[l];
1133: PetscFunctionReturn(PETSC_SUCCESS);
1134: }
1136: /*@
1137: BVRestoreColumn - Restore a column obtained with `BVGetColumn()`.
1139: Logically Collective
1141: Input Parameters:
1142: + bv - the basis vectors context
1143: . j - the index of the column
1144: - v - vector obtained with `BVGetColumn()`
1146: Note:
1147: The arguments must match the corresponding call to `BVGetColumn()`.
1149: Level: beginner
1151: .seealso: [](sec:bv), `BVGetColumn()`
1152: @*/
1153: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1154: {
1155: PetscObjectId id;
1156: PetscObjectState st;
1157: PetscInt l;
1159: PetscFunctionBegin;
1162: BVCheckSizes(bv,1);
1164: PetscAssertPointer(v,3);
1166: PetscCheck(j>=0 || -j<=bv->nc,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %" PetscInt_FMT " but only %" PetscInt_FMT " are available",-j,bv->nc);
1167: PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %" PetscInt_FMT " but only %" PetscInt_FMT " are available",j,bv->m);
1168: PetscCheck(j==bv->ci[0] || j==bv->ci[1],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Column %" PetscInt_FMT " has not been fetched with a call to BVGetColumn",j);
1169: l = (j==bv->ci[0])? 0: 1;
1170: PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1171: PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1172: PetscCall(VecGetState(*v,&st));
1173: if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1174: PetscUseTypeMethod(bv,restorecolumn,j,v);
1175: bv->ci[l] = -bv->nc-1;
1176: bv->st[l] = -1;
1177: bv->id[l] = 0;
1178: *v = NULL;
1179: PetscFunctionReturn(PETSC_SUCCESS);
1180: }
1182: /*@C
1183: BVGetArray - Returns a pointer to a contiguous array that contains this
1184: process' portion of the `BV` data.
1186: Logically Collective
1188: Input Parameter:
1189: . bv - the basis vectors context
1191: Output Parameter:
1192: . a - location to put pointer to the array
1194: Notes:
1195: `BVRestoreArray()` must be called when access to the array is no longer needed.
1196: This operation may imply a data copy, for `BV` types that do not store
1197: data contiguously in memory.
1199: The pointer will normally point to the first entry of the first column,
1200: but if the `BV` has constraints then these go before the regular columns.
1202: Note that for manipulating the pointer to the `BV` array, one must take into
1203: account the leading dimension, which might be different from the local
1204: number of rows, see `BVGetLeadingDimension()`.
1206: Use `BVGetArrayRead()` for read-only access.
1208: Level: advanced
1210: .seealso: [](sec:bv), `BVRestoreArray()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArrayRead()`, `BVType`
1211: @*/
1212: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1213: {
1214: PetscFunctionBegin;
1217: BVCheckSizes(bv,1);
1218: BVCheckOp(bv,1,getarray);
1219: PetscUseTypeMethod(bv,getarray,a);
1220: PetscFunctionReturn(PETSC_SUCCESS);
1221: }
1223: /*@C
1224: BVRestoreArray - Restore the `BV` object after `BVGetArray()` has been called.
1226: Logically Collective
1228: Input Parameters:
1229: + bv - the basis vectors context
1230: - a - location of pointer to array obtained from `BVGetArray()`
1232: Note:
1233: This operation may imply a data copy, for `BV` types that do not store
1234: data contiguously in memory.
1236: Level: advanced
1238: .seealso: [](sec:bv), `BVGetArray()`
1239: @*/
1240: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1241: {
1242: PetscFunctionBegin;
1245: BVCheckSizes(bv,1);
1246: PetscTryTypeMethod(bv,restorearray,a);
1247: if (a) *a = NULL;
1248: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1249: PetscFunctionReturn(PETSC_SUCCESS);
1250: }
1252: /*@C
1253: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1254: contains this processor's portion of the `BV` data.
1256: Not Collective
1258: Input Parameter:
1259: . bv - the basis vectors context
1261: Output Parameter:
1262: . a - location to put pointer to the array
1264: Notes:
1265: `BVRestoreArrayRead()` must be called when access to the array is no
1266: longer needed. This operation may imply a data copy, for `BV` types that
1267: do not store data contiguously in memory.
1269: The pointer will normally point to the first entry of the first column,
1270: but if the `BV` has constraints then these go before the regular columns.
1272: Level: advanced
1274: .seealso: [](sec:bv), `BVRestoreArray()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArray()`, `BVType`
1275: @*/
1276: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1277: {
1278: PetscFunctionBegin;
1281: BVCheckSizes(bv,1);
1282: BVCheckOp(bv,1,getarrayread);
1283: PetscUseTypeMethod(bv,getarrayread,a);
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: /*@C
1288: BVRestoreArrayRead - Restore the `BV` object after `BVGetArrayRead()` has
1289: been called.
1291: Not Collective
1293: Input Parameters:
1294: + bv - the basis vectors context
1295: - a - location of pointer to array obtained from `BVGetArrayRead()`
1297: Level: advanced
1299: .seealso: [](sec:bv), `BVGetArrayRead()`
1300: @*/
1301: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1302: {
1303: PetscFunctionBegin;
1306: BVCheckSizes(bv,1);
1307: PetscTryTypeMethod(bv,restorearrayread,a);
1308: if (a) *a = NULL;
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: /*@
1313: BVCreateVec - Creates a new `Vec` object with the same type and dimensions
1314: as the columns of the basis vectors object.
1316: Collective
1318: Input Parameter:
1319: . bv - the basis vectors context
1321: Output Parameter:
1322: . v - the new vector
1324: Note:
1325: The user is responsible for destroying the returned vector.
1327: Level: beginner
1329: .seealso: [](sec:bv), `BVCreateMat()`, `BVCreateVecEmpty()`
1330: @*/
1331: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1332: {
1333: PetscFunctionBegin;
1335: BVCheckSizes(bv,1);
1336: PetscAssertPointer(v,2);
1337: PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1338: PetscCall(VecSetLayout(*v,bv->map));
1339: PetscCall(VecSetType(*v,bv->vtype));
1340: PetscCall(VecSetUp(*v));
1341: PetscFunctionReturn(PETSC_SUCCESS);
1342: }
1344: /*@
1345: BVCreateVecEmpty - Creates a new `Vec` object with the same type and dimensions
1346: as the columns of the basis vectors object, but without internal array.
1348: Collective
1350: Input Parameter:
1351: . bv - the basis vectors context
1353: Output Parameter:
1354: . v - the new vector
1356: Note:
1357: This works as `BVCreateVec()`, but the new vector does not have the array allocated,
1358: so the intended usage is with `VecPlaceArray()`.
1360: Level: developer
1362: .seealso: [](sec:bv), `BVCreateVec()`
1363: @*/
1364: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1365: {
1366: PetscBool standard,cuda,hip,mpi;
1367: PetscInt N,nloc,bs;
1369: PetscFunctionBegin;
1371: BVCheckSizes(bv,1);
1372: PetscAssertPointer(v,2);
1374: PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1375: PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1376: PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
1377: if (standard || cuda || hip) {
1378: PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
1379: PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1380: PetscCall(PetscLayoutGetSize(bv->map,&N));
1381: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1382: if (cuda) {
1383: #if defined(PETSC_HAVE_CUDA)
1384: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1385: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1386: #endif
1387: } else if (hip) {
1388: #if defined(PETSC_HAVE_HIP)
1389: if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1390: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1391: #endif
1392: } else {
1393: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1394: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1395: }
1396: } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1397: PetscFunctionReturn(PETSC_SUCCESS);
1398: }
1400: /*@
1401: BVSetVecType - Set the vector type to be used when creating vectors via `BVCreateVec()`.
1403: Collective
1405: Input Parameters:
1406: + bv - the basis vectors context
1407: - vtype - the vector type
1409: Level: advanced
1411: Note:
1412: This is not needed if the `BV` object is set up with `BVSetSizesFromVec()`, but may be
1413: required in the case of `BVSetSizes()` if one wants to work with non-standard vectors.
1415: .seealso: [](sec:bv), `BVGetVecType()`, `BVSetSizesFromVec()`, `BVSetSizes()`
1416: @*/
1417: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1418: {
1419: PetscBool std;
1420: PetscMPIInt size;
1422: PetscFunctionBegin;
1424: PetscCall(PetscFree(bv->vtype));
1425: PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1426: if (std) {
1427: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1428: PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1429: } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1430: PetscFunctionReturn(PETSC_SUCCESS);
1431: }
1433: /*@
1434: BVGetVecType - Get the vector type to be used when creating vectors via `BVCreateVec()`.
1436: Not Collective
1438: Input Parameter:
1439: . bv - the basis vectors context
1441: Output Parameter:
1442: . vtype - the vector type
1444: Level: advanced
1446: .seealso: [](sec:bv), `BVSetVecType()`, `BVCreateVec()`
1447: @*/
1448: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1449: {
1450: PetscFunctionBegin;
1452: PetscAssertPointer(vtype,2);
1453: *vtype = bv->vtype;
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }
1457: /*@
1458: BVCreateMat - Creates a new `Mat` object of dense type and copies the contents
1459: of the `BV` object.
1461: Collective
1463: Input Parameter:
1464: . bv - the basis vectors context
1466: Output Parameter:
1467: . A - the new matrix
1469: Notes:
1470: The user is responsible for destroying the returned matrix.
1472: The matrix contains all columns of the `BV`, not just the active columns.
1474: Level: intermediate
1476: .seealso: [](sec:bv), `BVCreateFromMat()`, `BVCreateVec()`, `BVGetMat()`
1477: @*/
1478: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1479: {
1480: PetscInt ksave,lsave;
1481: Mat B;
1483: PetscFunctionBegin;
1485: BVCheckSizes(bv,1);
1486: PetscAssertPointer(A,2);
1488: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1489: lsave = bv->l;
1490: ksave = bv->k;
1491: bv->l = 0;
1492: bv->k = bv->m;
1493: PetscCall(BVGetMat(bv,&B));
1494: PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1495: PetscCall(BVRestoreMat(bv,&B));
1496: bv->l = lsave;
1497: bv->k = ksave;
1498: PetscFunctionReturn(PETSC_SUCCESS);
1499: }
1501: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1502: {
1503: PetscScalar *vv,*aa;
1504: PetscBool create=PETSC_FALSE;
1505: PetscInt m,cols;
1507: PetscFunctionBegin;
1508: m = bv->k-bv->l;
1509: if (!bv->Aget) create=PETSC_TRUE;
1510: else {
1511: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1512: PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1513: PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1514: if (cols!=m) {
1515: PetscCall(MatDestroy(&bv->Aget));
1516: create=PETSC_TRUE;
1517: }
1518: }
1519: PetscCall(BVGetArray(bv,&vv));
1520: if (create) {
1521: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,m,bv->ld,vv,&bv->Aget)); /* pass a pointer to avoid allocation of storage */
1522: PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1523: }
1524: PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld)); /* set the actual pointer */
1525: *A = bv->Aget;
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: /*@
1530: BVGetMat - Returns a `Mat` object of dense type that shares the memory of
1531: the `BV` object.
1533: Collective
1535: Input Parameter:
1536: . bv - the basis vectors context
1538: Output Parameter:
1539: . A - the matrix
1541: Notes:
1542: The returned matrix contains only the active columns. If the content of
1543: the `Mat` is modified, these changes are also done in the `BV` object. The
1544: user must call `BVRestoreMat()` when no longer needed.
1546: This operation implies a call to `BVGetArray()`, which may result in data
1547: copies.
1549: Level: advanced
1551: .seealso: [](sec:bv), `BVRestoreMat()`, `BVCreateMat()`, `BVGetArray()`
1552: @*/
1553: PetscErrorCode BVGetMat(BV bv,Mat *A)
1554: {
1555: char name[64];
1557: PetscFunctionBegin;
1559: BVCheckSizes(bv,1);
1560: PetscAssertPointer(A,2);
1561: PetscUseTypeMethod(bv,getmat,A);
1562: if (((PetscObject)bv)->name) { /* set A's name based on BV name */
1563: PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1564: PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1565: PetscCall(PetscObjectSetName((PetscObject)*A,name));
1566: }
1567: PetscFunctionReturn(PETSC_SUCCESS);
1568: }
1570: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1571: {
1572: PetscScalar *vv,*aa;
1574: PetscFunctionBegin;
1575: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1576: vv = aa-(bv->nc+bv->l)*bv->ld;
1577: PetscCall(MatDenseResetArray(bv->Aget));
1578: PetscCall(BVRestoreArray(bv,&vv));
1579: *A = NULL;
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@
1584: BVRestoreMat - Restores the Mat obtained with `BVGetMat()`.
1586: Logically Collective
1588: Input Parameters:
1589: + bv - the basis vectors context
1590: - A - the fetched matrix
1592: Note:
1593: A call to this function must match a previous call of `BVGetMat()`.
1594: The effect is that the contents of the `Mat` are copied back to the
1595: `BV` internal data structures.
1597: Level: advanced
1599: .seealso: [](sec:bv), `BVGetMat()`, `BVRestoreArray()`
1600: @*/
1601: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1602: {
1603: PetscFunctionBegin;
1605: BVCheckSizes(bv,1);
1606: PetscAssertPointer(A,2);
1607: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1608: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1609: PetscUseTypeMethod(bv,restoremat,A);
1610: PetscFunctionReturn(PETSC_SUCCESS);
1611: }
1613: /*
1614: Copy all user-provided attributes of V to another BV object W
1615: */
1616: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1617: {
1618: PetscFunctionBegin;
1619: PetscCall(PetscLayoutReference(V->map,&W->map));
1620: PetscCall(BVSetVecType(W,V->vtype));
1621: W->ld = V->ld;
1622: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1623: W->orthog_type = V->orthog_type;
1624: W->orthog_ref = V->orthog_ref;
1625: W->orthog_eta = V->orthog_eta;
1626: W->orthog_block = V->orthog_block;
1627: if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1628: W->matrix = V->matrix;
1629: W->indef = V->indef;
1630: W->vmm = V->vmm;
1631: W->rrandom = V->rrandom;
1632: W->deftol = V->deftol;
1633: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1634: W->rand = V->rand;
1635: W->sfocalled = V->sfocalled;
1636: PetscTryTypeMethod(V,duplicate,W);
1637: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: /*@
1642: BVDuplicate - Creates a new basis vectors object of the same type and
1643: dimensions as an existing one.
1645: Collective
1647: Input Parameter:
1648: . V - basis vectors context
1650: Output Parameter:
1651: . W - location to put the new `BV`
1653: Notes:
1654: The new `BV` has the same type and dimensions as `V`. Also, the inner
1655: product matrix and orthogonalization options are copied.
1657: `BVDuplicate()` DOES NOT COPY the entries, but rather allocates storage
1658: for the new basis vectors. Use `BVCopy()` to copy the content.
1660: Level: intermediate
1662: .seealso: [](sec:bv), `BVDuplicateResize()`, `BVCreate()`, `BVSetSizesFromVec()`, `BVCopy()`
1663: @*/
1664: PetscErrorCode BVDuplicate(BV V,BV *W)
1665: {
1666: PetscFunctionBegin;
1669: BVCheckSizes(V,1);
1670: PetscAssertPointer(W,2);
1671: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1672: (*W)->N = V->N;
1673: (*W)->n = V->n;
1674: (*W)->m = V->m;
1675: (*W)->k = V->m;
1676: PetscCall(BVDuplicate_Private(V,*W));
1677: PetscFunctionReturn(PETSC_SUCCESS);
1678: }
1680: /*@
1681: BVDuplicateResize - Creates a new basis vectors object of the same type and
1682: dimensions as an existing one, but with possibly different number of columns.
1684: Collective
1686: Input Parameters:
1687: + V - basis vectors context
1688: - m - the new number of columns
1690: Output Parameter:
1691: . W - location to put the new BV
1693: Note:
1694: This is equivalent of a call to `BVDuplicate()` followed by `BVResize()`. The
1695: contents of `V` are not copied to `W`.
1697: Level: intermediate
1699: .seealso: [](sec:bv), `BVDuplicate()`, `BVResize()`
1700: @*/
1701: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1702: {
1703: PetscFunctionBegin;
1706: BVCheckSizes(V,1);
1708: PetscAssertPointer(W,3);
1709: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1710: (*W)->N = V->N;
1711: (*W)->n = V->n;
1712: (*W)->m = m;
1713: (*W)->k = m;
1714: PetscCall(BVDuplicate_Private(V,*W));
1715: PetscFunctionReturn(PETSC_SUCCESS);
1716: }
1718: /*@
1719: BVGetCachedBV - Returns a `BV` object stored internally that holds the
1720: result of $BX$ after a call to `BVApplyMatrixBV()`.
1722: Collective
1724: Input Parameter:
1725: . bv - the basis vectors context
1727: Output Parameter:
1728: . cached - the cached `BV`
1730: Note:
1731: The cached `BV` is created if not available previously.
1733: Level: developer
1735: .seealso: [](sec:bv), `BVApplyMatrixBV()`
1736: @*/
1737: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1738: {
1739: PetscFunctionBegin;
1741: PetscAssertPointer(cached,2);
1742: BVCheckSizes(bv,1);
1743: if (!bv->cached) {
1744: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1745: bv->cached->N = bv->N;
1746: bv->cached->n = bv->n;
1747: bv->cached->m = bv->m;
1748: bv->cached->k = bv->m;
1749: PetscCall(BVDuplicate_Private(bv,bv->cached));
1750: }
1751: *cached = bv->cached;
1752: PetscFunctionReturn(PETSC_SUCCESS);
1753: }
1755: /*@
1756: BVCopy - Copies a basis vector object into another one, $W \leftarrow V$.
1758: Logically Collective
1760: Input Parameter:
1761: . V - basis vectors context
1763: Output Parameter:
1764: . W - the copy
1766: Note:
1767: Both `V` and `W` must be distributed in the same manner; local copies are
1768: done. Only active columns (excluding the leading ones) are copied.
1769: In the destination `W`, columns are overwritten starting from the leading ones.
1770: Constraints are not copied.
1772: Level: beginner
1774: .seealso: [](sec:bv), `BVCopyVec()`, `BVCopyColumn()`, `BVDuplicate()`, `BVSetActiveColumns()`
1775: @*/
1776: PetscErrorCode BVCopy(BV V,BV W)
1777: {
1778: PetscScalar *womega;
1779: const PetscScalar *vomega;
1781: PetscFunctionBegin;
1784: BVCheckSizes(V,1);
1785: BVCheckOp(V,1,copy);
1788: BVCheckSizes(W,2);
1789: PetscCheckSameTypeAndComm(V,1,W,2);
1790: PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
1791: PetscCheck(V->k-V->l<=W->m-W->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"W has %" PetscInt_FMT " non-leading columns, not enough to store %" PetscInt_FMT " columns",W->m-W->l,V->k-V->l);
1792: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1794: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1795: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1796: /* copy signature */
1797: PetscCall(BV_AllocateSignature(W));
1798: PetscCall(VecGetArrayRead(V->omega,&vomega));
1799: PetscCall(VecGetArray(W->omega,&womega));
1800: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1801: PetscCall(VecRestoreArray(W->omega,&womega));
1802: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1803: }
1804: PetscUseTypeMethod(V,copy,W);
1805: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1806: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1807: PetscFunctionReturn(PETSC_SUCCESS);
1808: }
1810: /*@
1811: BVCopyVec - Copies one of the columns of a basis vectors object into a `Vec`.
1813: Logically Collective
1815: Input Parameters:
1816: + V - basis vectors context
1817: - j - the index of the column to be copied
1819: Output Parameter:
1820: . w - the copied column
1822: Note:
1823: Both `V` and `w` must be distributed in the same manner; local copies are done.
1825: Level: beginner
1827: .seealso: [](sec:bv), `BVCopy()`, `BVCopyColumn()`, `BVInsertVec()`
1828: @*/
1829: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1830: {
1831: PetscInt n,N;
1832: Vec z;
1834: PetscFunctionBegin;
1837: BVCheckSizes(V,1);
1840: PetscCheckSameComm(V,1,w,3);
1842: PetscCall(VecGetSize(w,&N));
1843: PetscCall(VecGetLocalSize(w,&n));
1844: PetscCheck(N==V->N && n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ") do not match BV sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ")",N,n,V->N,V->n);
1846: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1847: PetscCall(BVGetColumn(V,j,&z));
1848: PetscCall(VecCopy(z,w));
1849: PetscCall(BVRestoreColumn(V,j,&z));
1850: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1851: PetscFunctionReturn(PETSC_SUCCESS);
1852: }
1854: /*@
1855: BVCopyColumn - Copies the values from one of the columns to another one.
1857: Logically Collective
1859: Input Parameters:
1860: + V - basis vectors context
1861: . j - the number of the source column
1862: - i - the number of the destination column
1864: Level: beginner
1866: .seealso: [](sec:bv), `BVCopy()`, `BVCopyVec()`
1867: @*/
1868: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1869: {
1870: PetscScalar *omega;
1872: PetscFunctionBegin;
1875: BVCheckSizes(V,1);
1878: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1880: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1881: if (V->omega) {
1882: PetscCall(VecGetArray(V->omega,&omega));
1883: omega[i] = omega[j];
1884: PetscCall(VecRestoreArray(V->omega,&omega));
1885: }
1886: PetscUseTypeMethod(V,copycolumn,j,i);
1887: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1888: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1889: PetscFunctionReturn(PETSC_SUCCESS);
1890: }
1892: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1893: {
1894: PetscInt ncols;
1896: PetscFunctionBegin;
1897: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1898: if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1899: if (!*split) {
1900: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1901: (*split)->issplit = left? 1: 2;
1902: (*split)->splitparent = bv;
1903: (*split)->N = bv->N;
1904: (*split)->n = bv->n;
1905: (*split)->m = bv->m;
1906: (*split)->k = bv->m;
1907: PetscCall(BVDuplicate_Private(bv,*split));
1908: }
1909: (*split)->l = 0;
1910: (*split)->k = left? bv->l: bv->k-bv->l;
1911: (*split)->nc = left? bv->nc: 0;
1912: (*split)->m = ncols-(*split)->nc;
1913: if ((*split)->nc) {
1914: (*split)->ci[0] = -(*split)->nc-1;
1915: (*split)->ci[1] = -(*split)->nc-1;
1916: }
1917: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1918: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1919: PetscFunctionReturn(PETSC_SUCCESS);
1920: }
1922: /*@
1923: BVGetSplit - Splits the `BV` object into two `BV` objects that share the
1924: internal data, one of them containing the leading columns and the other
1925: one containing the remaining columns.
1927: Collective
1929: Input Parameter:
1930: . bv - the basis vectors context
1932: Output Parameters:
1933: + L - left `BV` containing leading columns (can be `NULL`)
1934: - R - right `BV` containing remaining columns (can be `NULL`)
1936: Notes:
1937: The columns are split in two sets. The leading columns (including the
1938: constraints) are assigned to the left `BV` and the remaining columns
1939: are assigned to the right `BV`. The number of leading columns, as
1940: specified with `BVSetActiveColumns()`, must be between `1` and `m-1` (to
1941: guarantee that both `L` and `R` have at least one column).
1943: The returned `BV`s must be seen as references (not copies) of the input
1944: `BV`, that is, modifying them will change the entries of `bv` as well.
1945: The returned `BV`s must not be destroyed. `BVRestoreSplit()` must be called
1946: when they are no longer needed.
1948: Pass `NULL` for any of the output `BV`s that is not needed.
1950: Level: advanced
1952: .seealso: [](sec:bv), `BVRestoreSplit()`, `BVSetActiveColumns()`, `BVSetNumConstraints()`, `BVGetSplitRows()`
1953: @*/
1954: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1955: {
1956: PetscFunctionBegin;
1959: BVCheckSizes(bv,1);
1960: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1961: PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1962: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1963: bv->lsplit = bv->nc+bv->l;
1964: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1965: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1966: if (L) *L = bv->L;
1967: if (R) *R = bv->R;
1968: PetscFunctionReturn(PETSC_SUCCESS);
1969: }
1971: /*@
1972: BVRestoreSplit - Restore the BV objects obtained with `BVGetSplit()`.
1974: Logically Collective
1976: Input Parameters:
1977: + bv - the basis vectors context
1978: . L - left BV obtained with BVGetSplit()
1979: - R - right BV obtained with BVGetSplit()
1981: Note:
1982: The arguments must match the corresponding call to `BVGetSplit()`.
1984: Level: advanced
1986: .seealso: [](sec:bv), `BVGetSplit()`
1987: @*/
1988: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1989: {
1990: PetscFunctionBegin;
1993: BVCheckSizes(bv,1);
1994: PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1995: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1996: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
1997: PetscCheck(!L || ((*L)->ci[0]<=(*L)->nc-1 && (*L)->ci[1]<=(*L)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 2 has unrestored columns, use BVRestoreColumn()");
1998: PetscCheck(!R || ((*R)->ci[0]<=(*R)->nc-1 && (*R)->ci[1]<=(*R)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 3 has unrestored columns, use BVRestoreColumn()");
2000: PetscTryTypeMethod(bv,restoresplit,L,R);
2001: bv->lsplit = 0;
2002: if (L) *L = NULL;
2003: if (R) *R = NULL;
2004: PetscFunctionReturn(PETSC_SUCCESS);
2005: }
2007: /*
2008: Copy all user-provided attributes of V to another BV object W with different layout
2009: */
2010: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
2011: {
2012: PetscFunctionBegin;
2013: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2014: PetscCall(BVSetVecType(W,V->vtype));
2015: W->ld = V->ld;
2016: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2017: W->orthog_type = V->orthog_type;
2018: W->orthog_ref = V->orthog_ref;
2019: W->orthog_eta = V->orthog_eta;
2020: W->orthog_block = V->orthog_block;
2021: W->vmm = V->vmm;
2022: W->rrandom = V->rrandom;
2023: W->deftol = V->deftol;
2024: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
2025: W->rand = V->rand;
2026: W->sfocalled = V->sfocalled;
2027: PetscTryTypeMethod(V,duplicate,W);
2028: PetscCall(PetscObjectStateIncrease((PetscObject)W));
2029: PetscFunctionReturn(PETSC_SUCCESS);
2030: }
2032: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2033: {
2034: PetscInt rstart,rend,lstart;
2035: PetscInt N,n;
2036: PetscBool contig;
2038: PetscFunctionBegin;
2039: PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2040: PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2041: PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2042: if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2043: else bv->lsplit = -lstart;
2044: PetscCall(ISGetSize(is,&N));
2045: PetscCall(ISGetLocalSize(is,&n));
2046: if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2047: if (!*split) {
2048: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2049: (*split)->issplit = top? -1: -2;
2050: (*split)->splitparent = bv;
2051: (*split)->N = N;
2052: (*split)->n = n;
2053: (*split)->m = bv->m;
2054: PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2055: }
2056: (*split)->k = bv->k;
2057: (*split)->l = bv->l;
2058: (*split)->nc = bv->nc;
2059: if ((*split)->nc) {
2060: (*split)->ci[0] = -(*split)->nc-1;
2061: (*split)->ci[1] = -(*split)->nc-1;
2062: }
2063: if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2064: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2065: PetscFunctionReturn(PETSC_SUCCESS);
2066: }
2068: /*@
2069: BVGetSplitRows - Splits the `BV` object into two `BV` objects that share the
2070: internal data, using a disjoint horizontal splitting.
2072: Collective
2074: Input Parameters:
2075: + bv - the basis vectors context
2076: . isup - the index set that defines the upper part of the horizontal splitting
2077: - islo - the index set that defines the lower part of the horizontal splitting
2079: Output Parameters:
2080: + U - the resulting `BV` containing the upper rows
2081: - L - the resulting `BV` containing the lower rows
2083: Notes:
2084: The index sets must be such that every MPI process can extract the selected
2085: rows from its local part of the input `BV`, and this part must be contiguous.
2086: With one process, `isup` will list contiguous indices starting from 0, and `islo`
2087: will contain the remaining indices, hence we refer to upper and lower part.
2088: However, with several processes the indices will be interleaved because
2089: `isup` will refer to the upper part of the local array.
2091: The intended use of this function is with matrices of `MATNEST` type, where
2092: `MatNestGetISs()` will return the appropriate index sets.
2094: The returned `BV`s must be seen as references (not copies) of the input
2095: `BV`, that is, modifying them will change the entries of `bv` as well.
2096: The returned `BV`s must not be destroyed. `BVRestoreSplitRows()` must be called
2097: when they are no longer needed.
2099: Pass `NULL` for any of the output `BV`s that is not needed.
2101: Level: advanced
2103: .seealso: [](sec:bv), `BVRestoreSplitRows()`, `BVGetSplit()`
2104: @*/
2105: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2106: {
2107: PetscFunctionBegin;
2110: BVCheckSizes(bv,1);
2113: PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2114: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2115: if (U) {
2116: PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2117: *U = bv->R;
2118: }
2119: if (L) {
2120: PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2121: *L = bv->L;
2122: }
2123: PetscFunctionReturn(PETSC_SUCCESS);
2124: }
2126: /*@
2127: BVRestoreSplitRows - Restore the BV objects obtained with `BVGetSplitRows()`.
2129: Logically Collective
2131: Input Parameters:
2132: + bv - the basis vectors context
2133: . isup - the index set that defines the upper part of the horizontal splitting
2134: . islo - the index set that defines the lower part of the horizontal splitting
2135: . U - the `BV` containing the upper rows
2136: - L - the `BV` containing the lower rows
2138: Note:
2139: The arguments must match the corresponding call to `BVGetSplitRows()`.
2141: Level: advanced
2143: .seealso: [](sec:bv), `BVGetSplitRows()`
2144: @*/
2145: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2146: {
2147: PetscFunctionBegin;
2150: BVCheckSizes(bv,1);
2151: PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2152: PetscCheck(!U || ((*U)->ci[0]<=(*U)->nc-1 && (*U)->ci[1]<=(*U)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The upper BV has unrestored columns, use BVRestoreColumn()");
2153: PetscCheck(!L || ((*L)->ci[0]<=(*L)->nc-1 && (*L)->ci[1]<=(*L)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The lower BV has unrestored columns, use BVRestoreColumn()");
2154: PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2155: bv->lsplit = 0;
2156: PetscFunctionReturn(PETSC_SUCCESS);
2157: }
2159: /*@
2160: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2161: definite inner product.
2163: Logically Collective
2165: Input Parameters:
2166: + bv - basis vectors
2167: - deftol - the tolerance
2169: Options Database Key:
2170: . -bv_definite_tol <deftol> - the tolerance
2172: Notes:
2173: When using a non-standard inner product, see `BVSetMatrix()`, the solver needs
2174: to compute $\sqrt{z^*B z}$ for various vectors $z$. If the inner product has not
2175: been declared indefinite, the value $z^*B z$ must be positive, but due to
2176: rounding error a tiny value may become negative. A tolerance is used to
2177: detect this situation. Likewise, in complex arithmetic $z^*B z$ should be
2178: real, and we use the same tolerance to check whether a nonzero imaginary part
2179: can be considered negligible.
2181: This function sets this tolerance, which defaults to `10*PETSC_MACHINE_EPSILON`.
2182: The default value should be good for most applications.
2184: Level: advanced
2186: .seealso: [](sec:bv), `BVSetMatrix()`
2187: @*/
2188: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2189: {
2190: PetscFunctionBegin;
2193: if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2194: else {
2195: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2196: bv->deftol = deftol;
2197: }
2198: PetscFunctionReturn(PETSC_SUCCESS);
2199: }
2201: /*@
2202: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2203: inner product.
2205: Not Collective
2207: Input Parameter:
2208: . bv - the basis vectors
2210: Output Parameter:
2211: . deftol - the tolerance
2213: Level: advanced
2215: .seealso: [](sec:bv), `BVSetDefiniteTolerance()`
2216: @*/
2217: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2218: {
2219: PetscFunctionBegin;
2221: PetscAssertPointer(deftol,2);
2222: *deftol = bv->deftol;
2223: PetscFunctionReturn(PETSC_SUCCESS);
2224: }
2226: /*@
2227: BVSetLeadingDimension - Set the leading dimension to be used for storing the `BV` data.
2229: Not Collective
2231: Input Parameters:
2232: + bv - basis vectors
2233: - ld - the leading dimension
2235: Notes:
2236: This parameter is relevant for `BVMAT`, though it might be employed in other types
2237: as well.
2239: When the internal data of the `BV` is stored as a dense matrix, the leading dimension
2240: has the same meaning as in `MatDenseSetLDA()`, i.e., the distance in number of
2241: elements from one entry of the matrix to the one in the next column at the same
2242: row. The leading dimension refers to the local array, and hence can be different
2243: in different processes.
2245: The user does not need to change this parameter. The default value is equal to the
2246: number of local rows, but this value may be increased a little to guarantee alignment
2247: (especially in the case of GPU storage).
2249: Level: advanced
2251: .seealso: [](sec:bv), `BVGetLeadingDimension()`
2252: @*/
2253: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2254: {
2255: PetscFunctionBegin;
2258: PetscCheck((bv->n<0 && bv->N<0) || !bv->ops->create,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Must call BVSetLeadingDimension() before setting the BV type and sizes");
2259: if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
2260: else {
2261: PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2262: bv->ld = ld;
2263: }
2264: PetscFunctionReturn(PETSC_SUCCESS);
2265: }
2267: /*@
2268: BVGetLeadingDimension - Returns the leading dimension of the `BV`.
2270: Not Collective
2272: Input Parameter:
2273: . bv - the basis vectors
2275: Output Parameter:
2276: . ld - the leading dimension
2278: Level: advanced
2280: Notes:
2281: The returned value may be different in different processes.
2283: The leading dimension must be used when accessing the internal array via
2284: `BVGetArray()` or `BVGetArrayRead()`.
2286: .seealso: [](sec:bv), `BVSetLeadingDimension()`, `BVGetArray()`, `BVGetArrayRead()`
2287: @*/
2288: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2289: {
2290: PetscFunctionBegin;
2292: PetscAssertPointer(ld,2);
2293: *ld = bv->ld;
2294: PetscFunctionReturn(PETSC_SUCCESS);
2295: }