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: Notes:
163: Apart from the local and global sizes, the template vector `t` is also used
164: to get the vector type, which will be set as the `BV` vector type with
165: `BVSetVecType()`. This is relevant in cases where computation must be done
166: on the GPU.
168: If `t` is a nested vector (type `VECNEST`), then the vector type will be
169: obtained from the first subvector. The rationale is that `BV` will operate
170: internally with regular vectors, even though the user problem is defined via
171: nested matrices and vectors.
173: Level: beginner
175: .seealso: [](sec:bv), `BVSetSizes()`, `BVGetSizes()`, `BVResize()`, `BVSetVecType()`
176: @*/
177: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
178: {
179: PetscInt ma;
180: PetscLayout map;
181: VecType vtype;
182: PetscBool nest;
183: Vec v;
185: PetscFunctionBegin;
188: PetscCheckSameComm(bv,1,t,2);
190: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
191: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
192: PetscCall(PetscObjectTypeCompare((PetscObject)t,VECNEST,&nest));
193: if (nest) {
194: PetscCall(VecNestGetSubVec(t,0,&v));
195: PetscCall(VecGetType(v,&vtype));
196: } else PetscCall(VecGetType(t,&vtype));
197: PetscCall(BVSetVecType(bv,vtype));
198: PetscCall(VecGetLayout(t,&map));
199: PetscCall(PetscLayoutReference(map,&bv->map));
200: PetscCall(VecGetSize(t,&bv->N));
201: PetscCall(VecGetLocalSize(t,&bv->n));
202: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
203: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
204: 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);
205: }
206: bv->m = m;
207: bv->k = m;
208: if (bv->ops->create) {
209: PetscUseTypeMethod(bv,create);
210: bv->ops->create = NULL;
211: bv->defersfo = PETSC_FALSE;
212: }
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*@
217: BVGetSizes - Returns the local and global sizes, and the number of columns.
219: Not Collective
221: Input Parameter:
222: . bv - the basis vectors
224: Output Parameters:
225: + n - the local size
226: . N - the global size
227: - m - the number of columns
229: Note:
230: Normal usage requires that `bv` has already been given its sizes, otherwise
231: the call fails. However, this function can also be used to determine if
232: a `BV` object has been initialized completely (sizes and type). For this,
233: call with `n=NULL` and `N=NULL`, then a return value of `m=0` indicates that
234: the `BV` object is not ready for use yet.
236: Level: beginner
238: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
239: @*/
240: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
241: {
242: PetscFunctionBegin;
243: if (!bv) {
244: if (m && !n && !N) *m = 0;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: if (n || N) BVCheckSizes(bv,1);
249: if (n) *n = bv->n;
250: if (N) *N = bv->N;
251: if (m) *m = bv->m;
252: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: BVSetNumConstraints - Set the number of constraints.
259: Logically Collective
261: Input Parameters:
262: + V - basis vectors
263: - nc - number of constraints
265: Notes:
266: This function sets the number of constraints to `nc` and marks all remaining
267: columns as regular. Normal usage would be to call `BVInsertConstraints()` instead.
269: If `nc` is smaller than the previously set value, then some of the constraints
270: are discarded. In particular, using `nc=0` removes all constraints preserving
271: the content of regular columns.
273: Level: developer
275: .seealso: [](sec:bv), `BVInsertConstraints()`
276: @*/
277: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
278: {
279: PetscInt total,diff,i;
280: Vec x,y;
282: PetscFunctionBegin;
285: PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
287: BVCheckSizes(V,1);
288: PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
290: diff = nc-V->nc;
291: if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
292: total = V->nc+V->m;
293: PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
294: if (diff<0) { /* lessen constraints, shift contents of BV */
295: for (i=0;i<V->m;i++) {
296: PetscCall(BVGetColumn(V,i,&x));
297: PetscCall(BVGetColumn(V,i+diff,&y));
298: PetscCall(VecCopy(x,y));
299: PetscCall(BVRestoreColumn(V,i,&x));
300: PetscCall(BVRestoreColumn(V,i+diff,&y));
301: }
302: }
303: V->nc = nc;
304: V->ci[0] = -V->nc-1;
305: V->ci[1] = -V->nc-1;
306: V->m = total-nc;
307: V->l = PetscMin(V->l,V->m);
308: V->k = PetscMin(V->k,V->m);
309: PetscCall(PetscObjectStateIncrease((PetscObject)V));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: BVGetNumConstraints - Returns the number of constraints.
316: Not Collective
318: Input Parameter:
319: . bv - the basis vectors
321: Output Parameter:
322: . nc - the number of constraints
324: Level: advanced
326: .seealso: [](sec:bv), `BVGetSizes()`, `BVInsertConstraints()`
327: @*/
328: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
329: {
330: PetscFunctionBegin;
332: PetscAssertPointer(nc,2);
333: *nc = bv->nc;
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: /*@
338: BVResize - Change the number of columns.
340: Collective
342: Input Parameters:
343: + bv - the basis vectors
344: . m - the new number of columns
345: - copy - a flag indicating whether current values should be kept
347: Note:
348: Internal storage is reallocated. If the `copy` flag is set to true, then
349: the contents are copied to the leading part of the new space.
351: Level: advanced
353: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
354: @*/
355: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
356: {
357: PetscScalar *array;
358: const PetscScalar *omega;
359: Vec v;
361: PetscFunctionBegin;
366: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
367: PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
368: if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
369: BVCheckOp(bv,1,resize);
371: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
372: PetscUseTypeMethod(bv,resize,m,copy);
373: PetscCall(VecDestroy(&bv->buffer));
374: PetscCall(BVDestroy(&bv->cached));
375: PetscCall(PetscFree2(bv->h,bv->c));
376: if (bv->omega) {
377: if (bv->cuda) {
378: #if defined(PETSC_HAVE_CUDA)
379: PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
380: #else
381: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
382: #endif
383: } else if (bv->hip) {
384: #if defined(PETSC_HAVE_HIP)
385: PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
386: #else
387: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
388: #endif
389: } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
390: if (copy) {
391: PetscCall(VecGetArray(v,&array));
392: PetscCall(VecGetArrayRead(bv->omega,&omega));
393: PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
394: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
395: PetscCall(VecRestoreArray(v,&array));
396: } else PetscCall(VecSet(v,1.0));
397: PetscCall(VecDestroy(&bv->omega));
398: bv->omega = v;
399: }
400: bv->m = m;
401: bv->k = PetscMin(bv->k,m);
402: bv->l = PetscMin(bv->l,m);
403: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
404: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: BVSetActiveColumns - Specify the columns that will be involved in operations.
411: Logically Collective
413: Input Parameters:
414: + bv - the basis vectors context
415: . l - number of leading columns
416: - k - number of active columns
418: Notes:
419: In operations such as `BVMult()` or `BVDot()`, only the first `k` columns are
420: considered. This is useful when the `BV` is filled from left to right, so
421: the last `m-k` columns do not have relevant information.
423: Also in operations such as `BVMult()` or `BVDot()`, the first `l` columns are
424: normally not included in the computation. See the manual page of each
425: operation.
427: In orthogonalization operations, the first `l` columns are treated
428: differently, they participate in the orthogonalization but the computed
429: coefficients are not stored.
431: Use `PETSC_CURRENT` to leave any of the values unchanged. Use `PETSC_DETERMINE`
432: to set `l` to the minimum value (`0`) and `k` to the maximum (`m`).
434: Level: intermediate
436: .seealso: [](sec:bv), `BVGetActiveColumns()`, `BVSetSizes()`
437: @*/
438: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
439: {
440: PetscFunctionBegin;
444: BVCheckSizes(bv,1);
445: if (PetscUnlikely(k == PETSC_DETERMINE)) {
446: bv->k = bv->m;
447: } else if (k != PETSC_CURRENT) {
448: 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);
449: bv->k = k;
450: }
451: if (PetscUnlikely(l == PETSC_DETERMINE)) {
452: bv->l = 0;
453: } else if (l != PETSC_CURRENT) {
454: 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);
455: bv->l = l;
456: }
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@
461: BVGetActiveColumns - Returns the current active dimensions.
463: Not Collective
465: Input Parameter:
466: . bv - the basis vectors context
468: Output Parameters:
469: + l - number of leading columns
470: - k - number of active columns
472: Level: intermediate
474: .seealso: [](sec:bv), `BVSetActiveColumns()`
475: @*/
476: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
477: {
478: PetscFunctionBegin;
480: if (l) *l = bv->l;
481: if (k) *k = bv->k;
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: /*@
486: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
488: Collective
490: Input Parameters:
491: + bv - the basis vectors context
492: . B - a Hermitian matrix (may be `NULL`)
493: - indef - a flag indicating if the matrix is indefinite
495: Notes:
496: This is used to specify a non-standard inner product, whose matrix
497: representation is given by `B`. Then, all inner products required during
498: orthogonalization are computed as $(x,y)_B=y^*Bx$ rather than the
499: standard form $(x,y)=y^*x$.
501: Matrix `B` must be real symmetric (or complex Hermitian). A genuine inner
502: product requires that `B` is also positive (semi-)definite. However, we
503: also allow for an indefinite `B` (setting `indef=PETSC_TRUE`), in which
504: case the orthogonalization uses an indefinite inner product.
506: This affects operations `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, and variants.
508: Setting `B=NULL` has the same effect as if the identity matrix was passed.
510: Level: advanced
512: .seealso: [](sec:bv), `BVGetMatrix()`, `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, `BVSetDefiniteTolerance()`
513: @*/
514: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
515: {
516: PetscInt m,n;
518: PetscFunctionBegin;
521: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
522: if (B) {
524: PetscCall(MatGetLocalSize(B,&m,&n));
525: PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
526: 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);
527: }
528: if (B) PetscCall(PetscObjectReference((PetscObject)B));
529: PetscCall(MatDestroy(&bv->matrix));
530: bv->matrix = B;
531: bv->indef = indef;
532: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
533: if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
534: if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
535: }
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@
540: BVGetMatrix - Retrieves the matrix representation of the inner product.
542: Not Collective
544: Input Parameter:
545: . bv - the basis vectors context
547: Output Parameters:
548: + B - the matrix of the inner product (may be `NULL`)
549: - indef - the flag indicating if the matrix is indefinite
551: Level: advanced
553: .seealso: [](sec:bv), `BVSetMatrix()`
554: @*/
555: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
556: {
557: PetscFunctionBegin;
559: if (B) *B = bv->matrix;
560: if (indef) *indef = bv->indef;
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: /*@
565: BVApplyMatrix - Multiplies a vector by the matrix representation of the
566: inner product.
568: Neighbor-wise Collective
570: Input Parameters:
571: + bv - the basis vectors context
572: - x - the vector
574: Output Parameter:
575: . y - the result
577: Note:
578: If no matrix was specified this function copies the vector.
580: Level: advanced
582: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrixBV()`
583: @*/
584: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
585: {
586: PetscFunctionBegin;
590: if (bv->matrix) {
591: PetscCall(BV_IPMatMult(bv,x));
592: PetscCall(VecCopy(bv->Bx,y));
593: } else PetscCall(VecCopy(x,y));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: /*@
598: BVApplyMatrixBV - Multiplies the `BV` vectors by the matrix representation
599: of the inner product.
601: Neighbor-wise Collective
603: Input Parameter:
604: . X - the basis vectors context
606: Output Parameter:
607: . Y - the basis vectors to store the result (optional)
609: Note:
610: This function computes $Y = B X$, where $B$ is the matrix given with
611: `BVSetMatrix()`. This operation is computed as in `BVMatMult()`.
612: If no matrix was specified, then it just copies $Y = X$.
614: If no `Y` is given, the result is stored internally in the cached `BV`.
616: Level: developer
618: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrix()`, `BVMatMult()`, `BVGetCachedBV()`
619: @*/
620: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
621: {
622: PetscFunctionBegin;
624: if (Y) {
626: if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
627: else PetscCall(BVCopy(X,Y));
628: } else PetscCall(BV_IPMatMultBV(X));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /*@
633: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
635: Logically Collective
637: Input Parameters:
638: + bv - the basis vectors context
639: - omega - a vector representing the diagonal of the signature matrix
641: Note:
642: The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.
644: Level: developer
646: .seealso: [](sec:bv), `BVSetMatrix()`, `BVGetSignature()`
647: @*/
648: PetscErrorCode BVSetSignature(BV bv,Vec omega)
649: {
650: PetscInt i,n;
651: const PetscScalar *pomega;
652: PetscScalar *intern;
654: PetscFunctionBegin;
656: BVCheckSizes(bv,1);
660: PetscCall(VecGetSize(omega,&n));
661: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
662: PetscCall(BV_AllocateSignature(bv));
663: if (bv->indef) {
664: PetscCall(VecGetArrayRead(omega,&pomega));
665: PetscCall(VecGetArray(bv->omega,&intern));
666: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
667: PetscCall(VecRestoreArray(bv->omega,&intern));
668: PetscCall(VecRestoreArrayRead(omega,&pomega));
669: } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
670: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /*@
675: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
677: Not Collective
679: Input Parameter:
680: . bv - the basis vectors context
682: Output Parameter:
683: . omega - a vector representing the diagonal of the signature matrix
685: Note:
686: The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.
688: Level: developer
690: .seealso: [](sec:bv), `BVSetMatrix()`, `BVSetSignature()`
691: @*/
692: PetscErrorCode BVGetSignature(BV bv,Vec omega)
693: {
694: PetscInt i,n;
695: PetscScalar *pomega;
696: const PetscScalar *intern;
698: PetscFunctionBegin;
700: BVCheckSizes(bv,1);
704: PetscCall(VecGetSize(omega,&n));
705: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
706: if (bv->indef && bv->omega) {
707: PetscCall(VecGetArray(omega,&pomega));
708: PetscCall(VecGetArrayRead(bv->omega,&intern));
709: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
710: PetscCall(VecRestoreArrayRead(bv->omega,&intern));
711: PetscCall(VecRestoreArray(omega,&pomega));
712: } else PetscCall(VecSet(omega,1.0));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /*@
717: BVSetBufferVec - Attach a vector object to be used as buffer space for
718: several operations.
720: Collective
722: Input Parameters:
723: + bv - the basis vectors context)
724: - buffer - the vector
726: Notes:
727: Use `BVGetBufferVec()` to retrieve the vector (for example, to free it
728: at the end of the computations).
730: The vector must be sequential of length $(n_c+m)m$, where $m$ is the number
731: of columns of `bv` and $n_c$ is the number of constraints.
733: Level: developer
735: .seealso: [](sec:bv), `BVGetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`
736: @*/
737: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
738: {
739: PetscInt ld,n;
740: PetscMPIInt size;
742: PetscFunctionBegin;
745: BVCheckSizes(bv,1);
746: PetscCall(VecGetSize(buffer,&n));
747: ld = bv->m+bv->nc;
748: PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
749: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
750: PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
752: PetscCall(PetscObjectReference((PetscObject)buffer));
753: PetscCall(VecDestroy(&bv->buffer));
754: bv->buffer = buffer;
755: PetscFunctionReturn(PETSC_SUCCESS);
756: }
758: /*@
759: BVGetBufferVec - Obtain the buffer vector associated with the `BV` object.
761: Collective
763: Input Parameter:
764: . bv - the basis vectors context
766: Output Parameter:
767: . buffer - vector
769: Notes:
770: The vector is created if not available previously. It is a sequential vector
771: of length $(n_c+m) m$, where $m$ is the number of columns of `bv` and $n_c$ is the number
772: of constraints.
774: Developer Notes:
775: The buffer vector is viewed as a column-major matrix with leading dimension equal
776: to $n_c+m$, and $m$ columns at most. In the most common usage, it has the structure
777: $$
778: \left[\begin{array}{c|c}
779: s & \begin{array}{c} C \\\hline H \end{array}
780: \end{array}\right]
781: $$
782: where $H$ is an upper Hessenberg matrix of order $m (m-1)$, $C$ contains coefficients
783: related to orthogonalization against constraints (first $n_c$ rows), and $s$ is the
784: first column that contains scratch values computed during Gram-Schmidt
785: orthogonalization. In particular, `BVDotColumn()` and `BVMultColumn()` use $s$ to
786: store the coefficients.
788: Level: developer
790: .seealso: [](sec:bv), `BVSetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`, `BVDotColumn()`, `BVMultColumn()`
791: @*/
792: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
793: {
794: PetscInt ld;
796: PetscFunctionBegin;
798: PetscAssertPointer(buffer,2);
799: BVCheckSizes(bv,1);
800: if (!bv->buffer) {
801: ld = bv->m+bv->nc;
802: PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
803: PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
804: PetscCall(VecSetType(bv->buffer,bv->vtype));
805: }
806: *buffer = bv->buffer;
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /*@
811: BVSetRandomContext - Sets the `PetscRandom` object associated with the `BV`,
812: to be used in operations that need random numbers.
814: Collective
816: Input Parameters:
817: + bv - the basis vectors context
818: - rand - the random number generator context
820: Level: advanced
822: .seealso: [](sec:bv), `BVGetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`, `PetscRandomCreate()`
823: @*/
824: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
825: {
826: PetscFunctionBegin;
829: PetscCheckSameComm(bv,1,rand,2);
830: PetscCall(PetscObjectReference((PetscObject)rand));
831: PetscCall(PetscRandomDestroy(&bv->rand));
832: bv->rand = rand;
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: /*@
837: BVGetRandomContext - Gets the `PetscRandom` object associated with the `BV`.
839: Collective
841: Input Parameter:
842: . bv - the basis vectors context
844: Output Parameter:
845: . rand - the random number generator context
847: Level: advanced
849: .seealso: [](sec:bv), `BVSetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`
850: @*/
851: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
852: {
853: PetscFunctionBegin;
855: PetscAssertPointer(rand,2);
856: if (!bv->rand) {
857: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
858: if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
859: if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
860: if (bv->rrandom) {
861: PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
862: PetscCall(PetscRandomSeed(bv->rand));
863: }
864: }
865: *rand = bv->rand;
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*@
870: BVSetFromOptions - Sets `BV` options from the options database.
872: Collective
874: Input Parameter:
875: . bv - the basis vectors context
877: Note:
878: To see all options, run your program with the `-help` option.
880: Level: beginner
882: .seealso: [](sec:bv), `BVSetOptionsPrefix()`
883: @*/
884: PetscErrorCode BVSetFromOptions(BV bv)
885: {
886: char type[256];
887: PetscBool flg1,flg2,flg3,flg4;
888: PetscReal r;
889: BVOrthogType otype;
890: BVOrthogRefineType orefine;
891: BVOrthogBlockType oblock;
893: PetscFunctionBegin;
895: PetscCall(BVRegisterAll());
896: PetscObjectOptionsBegin((PetscObject)bv);
897: PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
898: if (flg1) PetscCall(BVSetType(bv,type));
899: else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));
901: otype = bv->orthog_type;
902: PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
903: orefine = bv->orthog_ref;
904: PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
905: r = bv->orthog_eta;
906: PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
907: oblock = bv->orthog_block;
908: PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
909: if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
911: PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
913: PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
914: if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
916: /* undocumented option to generate random vectors that are independent of the number of processes */
917: PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
919: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
920: else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
921: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
922: PetscOptionsEnd();
923: bv->sfocalled = PETSC_TRUE;
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: /*@
928: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
929: vectors (classical or modified Gram-Schmidt with or without refinement), and
930: for the block-orthogonalization (simultaneous orthogonalization of a set of
931: vectors).
933: Logically Collective
935: Input Parameters:
936: + bv - the basis vectors context
937: . type - the method of vector orthogonalization
938: . refine - type of refinement
939: . eta - parameter for selective refinement
940: - block - the method of block orthogonalization
942: Options Database Keys:
943: + -bv_orthog_type \<type\> - the vector orthogonalization, either `cgs` or `mgs`
944: . -bv_orthog_refine \<refine\> - the refinement type, `never`, `ifneeded` (default) or `always`
945: . -bv_orthog_eta \<eta\> - the value of `eta`
946: - -bv_orthog_block \<block\> - the block-orthogonalization method
948: Notes:
949: The default settings work well for most problems.
951: The parameter `eta` should be a real value between 0 and 1, that is used only when
952: the refinement type is `ifneeded`. Use `PETSC_DETERMINE` to set a reasonable
953: default value. Use `PETSC_CURRENT` to leave the current value unchanged.
955: When using several processes, `BV_ORTHOG_MGS` is likely to result in bad scalability.
957: If the method set for block orthogonalization is `BV_ORTHOG_BLOCK_GS`, then the computation
958: is done column by column with the vector orthogonalization.
960: Level: advanced
962: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVGetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
963: @*/
964: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
965: {
966: PetscFunctionBegin;
972: switch (type) {
973: case BV_ORTHOG_CGS:
974: case BV_ORTHOG_MGS:
975: bv->orthog_type = type;
976: break;
977: default:
978: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
979: }
980: switch (refine) {
981: case BV_ORTHOG_REFINE_NEVER:
982: case BV_ORTHOG_REFINE_IFNEEDED:
983: case BV_ORTHOG_REFINE_ALWAYS:
984: bv->orthog_ref = refine;
985: break;
986: default:
987: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
988: }
989: if (eta == (PetscReal)PETSC_DETERMINE) {
990: bv->orthog_eta = 0.7071;
991: } else if (eta != (PetscReal)PETSC_CURRENT) {
992: PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
993: bv->orthog_eta = eta;
994: }
995: switch (block) {
996: case BV_ORTHOG_BLOCK_GS:
997: case BV_ORTHOG_BLOCK_CHOL:
998: case BV_ORTHOG_BLOCK_TSQR:
999: case BV_ORTHOG_BLOCK_TSQRCHOL:
1000: case BV_ORTHOG_BLOCK_SVQB:
1001: bv->orthog_block = block;
1002: break;
1003: default:
1004: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
1005: }
1006: PetscFunctionReturn(PETSC_SUCCESS);
1007: }
1009: /*@
1010: BVGetOrthogonalization - Gets the orthogonalization settings from the `BV` object.
1012: Not Collective
1014: Input Parameter:
1015: . bv - basis vectors context
1017: Output Parameters:
1018: + type - the method of vector orthogonalization
1019: . refine - type of refinement
1020: . eta - parameter for selective refinement
1021: - block - the method of block orthogonalization
1023: Level: advanced
1025: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVSetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
1026: @*/
1027: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1028: {
1029: PetscFunctionBegin;
1031: if (type) *type = bv->orthog_type;
1032: if (refine) *refine = bv->orthog_ref;
1033: if (eta) *eta = bv->orthog_eta;
1034: if (block) *block = bv->orthog_block;
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: /*@
1039: BVSetMatMultMethod - Specifies the method used for the `BVMatMult()` operation.
1041: Logically Collective
1043: Input Parameters:
1044: + bv - the basis vectors context
1045: - method - the method for the `BVMatMult()` operation
1047: Options Database Key:
1048: . -bv_matmult \<method\> - choose one of the methods: `vecs`, `mat`
1050: Note:
1051: The default is `BV_MATMULT_MAT` except in the case of `BVVECS`.
1053: Level: advanced
1055: .seealso: [](sec:bv), `BVMatMult()`, `BVGetMatMultMethod()`, `BVMatMultType`
1056: @*/
1057: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1058: {
1059: PetscFunctionBegin;
1062: switch (method) {
1063: case BV_MATMULT_VECS:
1064: case BV_MATMULT_MAT:
1065: bv->vmm = method;
1066: break;
1067: case BV_MATMULT_MAT_SAVE:
1068: PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1069: bv->vmm = BV_MATMULT_MAT;
1070: break;
1071: default:
1072: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1073: }
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@
1078: BVGetMatMultMethod - Gets the method used for the `BVMatMult()` operation.
1080: Not Collective
1082: Input Parameter:
1083: . bv - basis vectors context
1085: Output Parameter:
1086: . method - the method for the `BVMatMult()` operation
1088: Level: advanced
1090: .seealso: [](sec:bv), `BVMatMult()`, `BVSetMatMultMethod()`, `BVMatMultType`
1091: @*/
1092: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1093: {
1094: PetscFunctionBegin;
1096: PetscAssertPointer(method,2);
1097: *method = bv->vmm;
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: /*@
1102: BVGetColumn - Returns a `Vec` object that contains the entries of the
1103: requested column of the basis vectors object.
1105: Logically Collective
1107: Input Parameters:
1108: + bv - the basis vectors context
1109: - j - the index of the requested column
1111: Output Parameter:
1112: . v - vector containing the jth column
1114: Notes:
1115: The returned `Vec` must be seen as a reference (not a copy) of the `BV`
1116: column, that is, modifying the `Vec` will change the `BV` entries as well.
1118: The returned `Vec` must not be destroyed. `BVRestoreColumn()` must be
1119: called when it is no longer needed. At most, two columns can be fetched,
1120: that is, this function can only be called twice before the corresponding
1121: `BVRestoreColumn()` is invoked.
1123: A negative index `j` selects the `i`-th constraint, where `i=-j`. Constraints
1124: should not be modified.
1126: Level: beginner
1128: .seealso: [](sec:bv), `BVRestoreColumn()`, `BVInsertConstraints()`
1129: @*/
1130: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1131: {
1132: PetscInt l;
1134: PetscFunctionBegin;
1137: BVCheckSizes(bv,1);
1138: BVCheckOp(bv,1,getcolumn);
1140: 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);
1141: 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);
1142: 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);
1143: l = BVAvailableVec;
1144: PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1145: PetscUseTypeMethod(bv,getcolumn,j,v);
1146: bv->ci[l] = j;
1147: PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
1148: PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1149: *v = bv->cv[l];
1150: PetscFunctionReturn(PETSC_SUCCESS);
1151: }
1153: /*@
1154: BVRestoreColumn - Restore a column obtained with `BVGetColumn()`.
1156: Logically Collective
1158: Input Parameters:
1159: + bv - the basis vectors context
1160: . j - the index of the column
1161: - v - vector obtained with `BVGetColumn()`
1163: Note:
1164: The arguments must match the corresponding call to `BVGetColumn()`.
1166: Level: beginner
1168: .seealso: [](sec:bv), `BVGetColumn()`
1169: @*/
1170: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1171: {
1172: PetscObjectId id;
1173: PetscObjectState st;
1174: PetscInt l;
1176: PetscFunctionBegin;
1179: BVCheckSizes(bv,1);
1181: PetscAssertPointer(v,3);
1183: 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);
1184: 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);
1185: 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);
1186: l = (j==bv->ci[0])? 0: 1;
1187: PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1188: PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1189: PetscCall(VecGetState(*v,&st));
1190: if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1191: PetscUseTypeMethod(bv,restorecolumn,j,v);
1192: bv->ci[l] = -bv->nc-1;
1193: bv->st[l] = -1;
1194: bv->id[l] = 0;
1195: *v = NULL;
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@C
1200: BVGetArray - Returns a pointer to a contiguous array that contains this
1201: process' portion of the `BV` data.
1203: Logically Collective
1205: Input Parameter:
1206: . bv - the basis vectors context
1208: Output Parameter:
1209: . a - location to put pointer to the array
1211: Notes:
1212: `BVRestoreArray()` must be called when access to the array is no longer needed.
1213: This operation may imply a data copy, for `BV` types that do not store
1214: data contiguously in memory.
1216: The pointer will normally point to the first entry of the first column,
1217: but if the `BV` has constraints then these go before the regular columns.
1219: Note that for manipulating the pointer to the `BV` array, one must take into
1220: account the leading dimension, which might be different from the local
1221: number of rows, see `BVGetLeadingDimension()`.
1223: Use `BVGetArrayRead()` for read-only access.
1225: Level: advanced
1227: .seealso: [](sec:bv), `BVRestoreArray()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArrayRead()`, `BVType`
1228: @*/
1229: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1230: {
1231: PetscFunctionBegin;
1234: BVCheckSizes(bv,1);
1235: BVCheckOp(bv,1,getarray);
1236: PetscUseTypeMethod(bv,getarray,a);
1237: PetscFunctionReturn(PETSC_SUCCESS);
1238: }
1240: /*@C
1241: BVRestoreArray - Restore the `BV` object after `BVGetArray()` has been called.
1243: Logically Collective
1245: Input Parameters:
1246: + bv - the basis vectors context
1247: - a - location of pointer to array obtained from `BVGetArray()`
1249: Note:
1250: This operation may imply a data copy, for `BV` types that do not store
1251: data contiguously in memory.
1253: Level: advanced
1255: .seealso: [](sec:bv), `BVGetArray()`
1256: @*/
1257: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1258: {
1259: PetscFunctionBegin;
1262: BVCheckSizes(bv,1);
1263: PetscTryTypeMethod(bv,restorearray,a);
1264: if (a) *a = NULL;
1265: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: /*@C
1270: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1271: contains this processor's portion of the `BV` data.
1273: Not Collective
1275: Input Parameter:
1276: . bv - the basis vectors context
1278: Output Parameter:
1279: . a - location to put pointer to the array
1281: Notes:
1282: `BVRestoreArrayRead()` must be called when access to the array is no
1283: longer needed. This operation may imply a data copy, for `BV` types that
1284: do not store data contiguously in memory.
1286: The pointer will normally point to the first entry of the first column,
1287: but if the `BV` has constraints then these go before the regular columns.
1289: Level: advanced
1291: .seealso: [](sec:bv), `BVRestoreArrayRead()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArray()`, `BVType`
1292: @*/
1293: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1294: {
1295: PetscFunctionBegin;
1298: BVCheckSizes(bv,1);
1299: BVCheckOp(bv,1,getarrayread);
1300: PetscUseTypeMethod(bv,getarrayread,a);
1301: PetscFunctionReturn(PETSC_SUCCESS);
1302: }
1304: /*@C
1305: BVRestoreArrayRead - Restore the `BV` object after `BVGetArrayRead()` has
1306: been called.
1308: Not Collective
1310: Input Parameters:
1311: + bv - the basis vectors context
1312: - a - location of pointer to array obtained from `BVGetArrayRead()`
1314: Level: advanced
1316: .seealso: [](sec:bv), `BVGetArrayRead()`
1317: @*/
1318: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1319: {
1320: PetscFunctionBegin;
1323: BVCheckSizes(bv,1);
1324: PetscTryTypeMethod(bv,restorearrayread,a);
1325: if (a) *a = NULL;
1326: PetscFunctionReturn(PETSC_SUCCESS);
1327: }
1329: /*@
1330: BVCreateVec - Creates a new `Vec` object with the same type and dimensions
1331: as the columns of the basis vectors object.
1333: Collective
1335: Input Parameter:
1336: . bv - the basis vectors context
1338: Output Parameter:
1339: . v - the new vector
1341: Note:
1342: The user is responsible for destroying the returned vector.
1344: Level: beginner
1346: .seealso: [](sec:bv), `BVCreateMat()`, `BVCreateVecEmpty()`
1347: @*/
1348: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1349: {
1350: PetscFunctionBegin;
1352: BVCheckSizes(bv,1);
1353: PetscAssertPointer(v,2);
1354: PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1355: PetscCall(VecSetLayout(*v,bv->map));
1356: PetscCall(VecSetType(*v,bv->vtype));
1357: PetscCall(VecSetUp(*v));
1358: PetscFunctionReturn(PETSC_SUCCESS);
1359: }
1361: /*@
1362: BVCreateVecEmpty - Creates a new `Vec` object with the same type and dimensions
1363: as the columns of the basis vectors object, but without internal array.
1365: Collective
1367: Input Parameter:
1368: . bv - the basis vectors context
1370: Output Parameter:
1371: . v - the new vector
1373: Note:
1374: This works as `BVCreateVec()`, but the new vector does not have the array allocated,
1375: so the intended usage is with `VecPlaceArray()`.
1377: Level: developer
1379: .seealso: [](sec:bv), `BVCreateVec()`
1380: @*/
1381: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1382: {
1383: PetscBool standard,cuda,hip,mpi;
1384: PetscInt N,nloc,bs;
1386: PetscFunctionBegin;
1388: BVCheckSizes(bv,1);
1389: PetscAssertPointer(v,2);
1391: PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1392: PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1393: PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
1394: if (standard || cuda || hip) {
1395: PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
1396: PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1397: PetscCall(PetscLayoutGetSize(bv->map,&N));
1398: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1399: if (cuda) {
1400: #if defined(PETSC_HAVE_CUDA)
1401: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1402: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1403: #endif
1404: } else if (hip) {
1405: #if defined(PETSC_HAVE_HIP)
1406: if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1407: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1408: #endif
1409: } else {
1410: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1411: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1412: }
1413: } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@
1418: BVSetVecType - Set the vector type to be used when creating vectors via `BVCreateVec()`.
1420: Collective
1422: Input Parameters:
1423: + bv - the basis vectors context
1424: - vtype - the vector type
1426: Level: advanced
1428: Note:
1429: This is not needed if the `BV` object is set up with `BVSetSizesFromVec()`, but may be
1430: required in the case of `BVSetSizes()` if one wants to work with non-standard vectors.
1432: .seealso: [](sec:bv), `BVCreateVec()`, `BVGetVecType()`, `BVSetSizesFromVec()`, `BVSetSizes()`
1433: @*/
1434: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1435: {
1436: PetscBool std;
1437: PetscMPIInt size;
1439: PetscFunctionBegin;
1441: PetscCall(PetscFree(bv->vtype));
1442: PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1443: if (std) {
1444: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1445: PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1446: } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: /*@
1451: BVGetVecType - Get the vector type to be used when creating vectors via `BVCreateVec()`.
1453: Not Collective
1455: Input Parameter:
1456: . bv - the basis vectors context
1458: Output Parameter:
1459: . vtype - the vector type
1461: Level: advanced
1463: .seealso: [](sec:bv), `BVSetVecType()`, `BVCreateVec()`
1464: @*/
1465: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1466: {
1467: PetscFunctionBegin;
1469: PetscAssertPointer(vtype,2);
1470: *vtype = bv->vtype;
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: /*@
1475: BVCreateMat - Creates a new `Mat` object of dense type and copies the contents
1476: of the `BV` object.
1478: Collective
1480: Input Parameter:
1481: . bv - the basis vectors context
1483: Output Parameter:
1484: . A - the new matrix
1486: Notes:
1487: The user is responsible for destroying the returned matrix.
1489: The matrix contains all columns of the `BV`, not just the active columns.
1491: Level: intermediate
1493: .seealso: [](sec:bv), `BVCreateFromMat()`, `BVCreateVec()`, `BVGetMat()`
1494: @*/
1495: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1496: {
1497: PetscInt ksave,lsave;
1498: Mat B;
1500: PetscFunctionBegin;
1502: BVCheckSizes(bv,1);
1503: PetscAssertPointer(A,2);
1505: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1506: lsave = bv->l;
1507: ksave = bv->k;
1508: bv->l = 0;
1509: bv->k = bv->m;
1510: PetscCall(BVGetMat(bv,&B));
1511: PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1512: PetscCall(BVRestoreMat(bv,&B));
1513: bv->l = lsave;
1514: bv->k = ksave;
1515: PetscFunctionReturn(PETSC_SUCCESS);
1516: }
1518: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1519: {
1520: PetscScalar *vv,*aa;
1521: PetscBool create=PETSC_FALSE;
1522: PetscInt m,cols;
1524: PetscFunctionBegin;
1525: m = bv->k-bv->l;
1526: if (!bv->Aget) create=PETSC_TRUE;
1527: else {
1528: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1529: PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1530: PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1531: if (cols!=m) {
1532: PetscCall(MatDestroy(&bv->Aget));
1533: create=PETSC_TRUE;
1534: }
1535: }
1536: PetscCall(BVGetArray(bv,&vv));
1537: if (create) {
1538: 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 */
1539: PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1540: }
1541: PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld)); /* set the actual pointer */
1542: *A = bv->Aget;
1543: PetscFunctionReturn(PETSC_SUCCESS);
1544: }
1546: /*@
1547: BVGetMat - Returns a `Mat` object of dense type that shares the memory of
1548: the `BV` object.
1550: Collective
1552: Input Parameter:
1553: . bv - the basis vectors context
1555: Output Parameter:
1556: . A - the matrix
1558: Notes:
1559: The returned matrix contains only the active columns. If the content of
1560: the `Mat` is modified, these changes are also done in the `BV` object. The
1561: user must call `BVRestoreMat()` when no longer needed.
1563: This operation implies a call to `BVGetArray()`, which may result in data
1564: copies.
1566: Level: advanced
1568: .seealso: [](sec:bv), `BVRestoreMat()`, `BVCreateMat()`, `BVGetArray()`
1569: @*/
1570: PetscErrorCode BVGetMat(BV bv,Mat *A)
1571: {
1572: char name[64];
1574: PetscFunctionBegin;
1576: BVCheckSizes(bv,1);
1577: PetscAssertPointer(A,2);
1578: PetscUseTypeMethod(bv,getmat,A);
1579: if (((PetscObject)bv)->name) { /* set A's name based on BV name */
1580: PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1581: PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1582: PetscCall(PetscObjectSetName((PetscObject)*A,name));
1583: }
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1588: {
1589: PetscScalar *vv,*aa;
1591: PetscFunctionBegin;
1592: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1593: vv = aa-(bv->nc+bv->l)*bv->ld;
1594: PetscCall(MatDenseResetArray(bv->Aget));
1595: PetscCall(BVRestoreArray(bv,&vv));
1596: *A = NULL;
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /*@
1601: BVRestoreMat - Restores the Mat obtained with `BVGetMat()`.
1603: Logically Collective
1605: Input Parameters:
1606: + bv - the basis vectors context
1607: - A - the fetched matrix
1609: Note:
1610: A call to this function must match a previous call of `BVGetMat()`.
1611: The effect is that the contents of the `Mat` are copied back to the
1612: `BV` internal data structures.
1614: Level: advanced
1616: .seealso: [](sec:bv), `BVGetMat()`, `BVRestoreArray()`
1617: @*/
1618: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1619: {
1620: PetscFunctionBegin;
1622: BVCheckSizes(bv,1);
1623: PetscAssertPointer(A,2);
1624: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1625: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1626: PetscUseTypeMethod(bv,restoremat,A);
1627: PetscFunctionReturn(PETSC_SUCCESS);
1628: }
1630: /*
1631: Copy all user-provided attributes of V to another BV object W
1632: */
1633: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1634: {
1635: PetscFunctionBegin;
1636: PetscCall(PetscLayoutReference(V->map,&W->map));
1637: PetscCall(BVSetVecType(W,V->vtype));
1638: W->ld = V->ld;
1639: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1640: W->orthog_type = V->orthog_type;
1641: W->orthog_ref = V->orthog_ref;
1642: W->orthog_eta = V->orthog_eta;
1643: W->orthog_block = V->orthog_block;
1644: if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1645: W->matrix = V->matrix;
1646: W->indef = V->indef;
1647: W->vmm = V->vmm;
1648: W->rrandom = V->rrandom;
1649: W->deftol = V->deftol;
1650: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1651: W->rand = V->rand;
1652: W->sfocalled = V->sfocalled;
1653: PetscTryTypeMethod(V,duplicate,W);
1654: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1655: PetscFunctionReturn(PETSC_SUCCESS);
1656: }
1658: /*@
1659: BVDuplicate - Creates a new basis vectors object of the same type and
1660: dimensions as an existing one.
1662: Collective
1664: Input Parameter:
1665: . V - basis vectors context
1667: Output Parameter:
1668: . W - location to put the new `BV`
1670: Notes:
1671: The new `BV` has the same type and dimensions as `V`. Also, the inner
1672: product matrix and orthogonalization options are copied.
1674: `BVDuplicate()` DOES NOT COPY the entries, but rather allocates storage
1675: for the new basis vectors. Use `BVCopy()` to copy the content.
1677: Level: intermediate
1679: .seealso: [](sec:bv), `BVDuplicateResize()`, `BVCreate()`, `BVSetSizesFromVec()`, `BVCopy()`
1680: @*/
1681: PetscErrorCode BVDuplicate(BV V,BV *W)
1682: {
1683: PetscFunctionBegin;
1686: BVCheckSizes(V,1);
1687: PetscAssertPointer(W,2);
1688: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1689: (*W)->N = V->N;
1690: (*W)->n = V->n;
1691: (*W)->m = V->m;
1692: (*W)->k = V->m;
1693: PetscCall(BVDuplicate_Private(V,*W));
1694: PetscFunctionReturn(PETSC_SUCCESS);
1695: }
1697: /*@
1698: BVDuplicateResize - Creates a new basis vectors object of the same type and
1699: dimensions as an existing one, but with possibly different number of columns.
1701: Collective
1703: Input Parameters:
1704: + V - basis vectors context
1705: - m - the new number of columns
1707: Output Parameter:
1708: . W - location to put the new BV
1710: Note:
1711: This is equivalent to a call to `BVDuplicate()` followed by `BVResize()`. The
1712: contents of `V` are not copied to `W`.
1714: Level: intermediate
1716: .seealso: [](sec:bv), `BVDuplicate()`, `BVResize()`
1717: @*/
1718: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1719: {
1720: PetscFunctionBegin;
1723: BVCheckSizes(V,1);
1725: PetscAssertPointer(W,3);
1726: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1727: (*W)->N = V->N;
1728: (*W)->n = V->n;
1729: (*W)->m = m;
1730: (*W)->k = m;
1731: PetscCall(BVDuplicate_Private(V,*W));
1732: PetscFunctionReturn(PETSC_SUCCESS);
1733: }
1735: /*@
1736: BVGetCachedBV - Returns a `BV` object stored internally that holds the
1737: result of $BX$ after a call to `BVApplyMatrixBV()`.
1739: Collective
1741: Input Parameter:
1742: . bv - the basis vectors context
1744: Output Parameter:
1745: . cached - the cached `BV`
1747: Note:
1748: The cached `BV` is created if not available previously.
1750: Level: developer
1752: .seealso: [](sec:bv), `BVApplyMatrixBV()`
1753: @*/
1754: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1755: {
1756: PetscFunctionBegin;
1758: PetscAssertPointer(cached,2);
1759: BVCheckSizes(bv,1);
1760: if (!bv->cached) {
1761: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1762: bv->cached->N = bv->N;
1763: bv->cached->n = bv->n;
1764: bv->cached->m = bv->m;
1765: bv->cached->k = bv->m;
1766: PetscCall(BVDuplicate_Private(bv,bv->cached));
1767: }
1768: *cached = bv->cached;
1769: PetscFunctionReturn(PETSC_SUCCESS);
1770: }
1772: /*@
1773: BVCopy - Copies a basis vector object into another one, $W \leftarrow V$.
1775: Logically Collective
1777: Input Parameter:
1778: . V - basis vectors context
1780: Output Parameter:
1781: . W - the copy
1783: Note:
1784: Both `V` and `W` must be distributed in the same manner; local copies are
1785: done. Only active columns (excluding the leading ones) are copied.
1786: In the destination `W`, columns are overwritten starting from the leading ones.
1787: Constraints are not copied.
1789: Level: beginner
1791: .seealso: [](sec:bv), `BVCopyVec()`, `BVCopyColumn()`, `BVDuplicate()`, `BVSetActiveColumns()`
1792: @*/
1793: PetscErrorCode BVCopy(BV V,BV W)
1794: {
1795: PetscScalar *womega;
1796: const PetscScalar *vomega;
1798: PetscFunctionBegin;
1801: BVCheckSizes(V,1);
1802: BVCheckOp(V,1,copy);
1805: BVCheckSizes(W,2);
1806: PetscCheckSameTypeAndComm(V,1,W,2);
1807: 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);
1808: 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);
1809: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1811: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1812: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1813: /* copy signature */
1814: PetscCall(BV_AllocateSignature(W));
1815: PetscCall(VecGetArrayRead(V->omega,&vomega));
1816: PetscCall(VecGetArray(W->omega,&womega));
1817: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1818: PetscCall(VecRestoreArray(W->omega,&womega));
1819: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1820: }
1821: PetscUseTypeMethod(V,copy,W);
1822: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1823: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1824: PetscFunctionReturn(PETSC_SUCCESS);
1825: }
1827: /*@
1828: BVCopyVec - Copies one of the columns of a basis vectors object into a `Vec`.
1830: Logically Collective
1832: Input Parameters:
1833: + V - basis vectors context
1834: - j - the index of the column to be copied
1836: Output Parameter:
1837: . w - the copied column
1839: Note:
1840: Both `V` and `w` must be distributed in the same manner; local copies are done.
1842: Level: beginner
1844: .seealso: [](sec:bv), `BVCopy()`, `BVCopyColumn()`, `BVInsertVec()`
1845: @*/
1846: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1847: {
1848: PetscInt n,N;
1849: Vec z;
1851: PetscFunctionBegin;
1854: BVCheckSizes(V,1);
1857: PetscCheckSameComm(V,1,w,3);
1859: PetscCall(VecGetSize(w,&N));
1860: PetscCall(VecGetLocalSize(w,&n));
1861: 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);
1863: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1864: PetscCall(BVGetColumn(V,j,&z));
1865: PetscCall(VecCopy(z,w));
1866: PetscCall(BVRestoreColumn(V,j,&z));
1867: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1868: PetscFunctionReturn(PETSC_SUCCESS);
1869: }
1871: /*@
1872: BVCopyColumn - Copies the values from one of the columns to another one.
1874: Logically Collective
1876: Input Parameters:
1877: + V - basis vectors context
1878: . j - the index of the source column
1879: - i - the index of the destination column
1881: Level: beginner
1883: .seealso: [](sec:bv), `BVCopy()`, `BVCopyVec()`
1884: @*/
1885: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1886: {
1887: PetscScalar *omega;
1889: PetscFunctionBegin;
1892: BVCheckSizes(V,1);
1895: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1897: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1898: if (V->omega) {
1899: PetscCall(VecGetArray(V->omega,&omega));
1900: omega[i] = omega[j];
1901: PetscCall(VecRestoreArray(V->omega,&omega));
1902: }
1903: PetscUseTypeMethod(V,copycolumn,j,i);
1904: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1905: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1906: PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1909: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1910: {
1911: PetscInt ncols;
1913: PetscFunctionBegin;
1914: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1915: if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1916: if (!*split) {
1917: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1918: (*split)->issplit = left? 1: 2;
1919: (*split)->splitparent = bv;
1920: (*split)->N = bv->N;
1921: (*split)->n = bv->n;
1922: (*split)->m = bv->m;
1923: (*split)->k = bv->m;
1924: PetscCall(BVDuplicate_Private(bv,*split));
1925: }
1926: (*split)->l = 0;
1927: (*split)->k = left? bv->l: bv->k-bv->l;
1928: (*split)->nc = left? bv->nc: 0;
1929: (*split)->m = ncols-(*split)->nc;
1930: if ((*split)->nc) {
1931: (*split)->ci[0] = -(*split)->nc-1;
1932: (*split)->ci[1] = -(*split)->nc-1;
1933: }
1934: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1935: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1936: PetscFunctionReturn(PETSC_SUCCESS);
1937: }
1939: /*@
1940: BVGetSplit - Splits the `BV` object into two `BV` objects that share the
1941: internal data, one of them containing the leading columns and the other
1942: one containing the remaining columns.
1944: Collective
1946: Input Parameter:
1947: . bv - the basis vectors context
1949: Output Parameters:
1950: + L - left `BV` containing leading columns (can be `NULL`)
1951: - R - right `BV` containing remaining columns (can be `NULL`)
1953: Notes:
1954: The columns are split in two sets. The leading columns (including the
1955: constraints) are assigned to the left `BV` and the remaining columns
1956: are assigned to the right `BV`. The number of leading columns, as
1957: specified with `BVSetActiveColumns()`, must be between `1` and `m-1` (to
1958: guarantee that both `L` and `R` have at least one column).
1960: The returned `BV`s must be seen as references (not copies) of the input
1961: `BV`, that is, modifying them will change the entries of `bv` as well.
1962: The returned `BV`s must not be destroyed. `BVRestoreSplit()` must be called
1963: when they are no longer needed.
1965: Pass `NULL` for any of the output `BV`s that is not needed.
1967: Level: advanced
1969: .seealso: [](sec:bv), `BVRestoreSplit()`, `BVSetActiveColumns()`, `BVSetNumConstraints()`, `BVGetSplitRows()`
1970: @*/
1971: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1972: {
1973: PetscFunctionBegin;
1976: BVCheckSizes(bv,1);
1977: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1978: PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1979: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1980: bv->lsplit = bv->nc+bv->l;
1981: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1982: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1983: if (L) *L = bv->L;
1984: if (R) *R = bv->R;
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: /*@
1989: BVRestoreSplit - Restore the BV objects obtained with `BVGetSplit()`.
1991: Logically Collective
1993: Input Parameters:
1994: + bv - the basis vectors context
1995: . L - left BV obtained with BVGetSplit()
1996: - R - right BV obtained with BVGetSplit()
1998: Note:
1999: The arguments must match the corresponding call to `BVGetSplit()`.
2001: Level: advanced
2003: .seealso: [](sec:bv), `BVGetSplit()`
2004: @*/
2005: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
2006: {
2007: PetscFunctionBegin;
2010: BVCheckSizes(bv,1);
2011: PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
2012: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
2013: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
2014: 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()");
2015: 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()");
2017: PetscTryTypeMethod(bv,restoresplit,L,R);
2018: bv->lsplit = 0;
2019: if (L) *L = NULL;
2020: if (R) *R = NULL;
2021: PetscFunctionReturn(PETSC_SUCCESS);
2022: }
2024: /*
2025: Copy all user-provided attributes of V to another BV object W with different layout
2026: */
2027: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
2028: {
2029: PetscFunctionBegin;
2030: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2031: PetscCall(BVSetVecType(W,V->vtype));
2032: W->ld = V->ld;
2033: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2034: W->orthog_type = V->orthog_type;
2035: W->orthog_ref = V->orthog_ref;
2036: W->orthog_eta = V->orthog_eta;
2037: W->orthog_block = V->orthog_block;
2038: W->vmm = V->vmm;
2039: W->rrandom = V->rrandom;
2040: W->deftol = V->deftol;
2041: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
2042: W->rand = V->rand;
2043: W->sfocalled = V->sfocalled;
2044: PetscTryTypeMethod(V,duplicate,W);
2045: PetscCall(PetscObjectStateIncrease((PetscObject)W));
2046: PetscFunctionReturn(PETSC_SUCCESS);
2047: }
2049: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2050: {
2051: PetscInt rstart,rend,lstart;
2052: PetscInt N,n;
2053: PetscBool contig;
2055: PetscFunctionBegin;
2056: PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2057: PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2058: PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2059: if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2060: else bv->lsplit = -lstart;
2061: PetscCall(ISGetSize(is,&N));
2062: PetscCall(ISGetLocalSize(is,&n));
2063: if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2064: if (!*split) {
2065: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2066: (*split)->issplit = top? -1: -2;
2067: (*split)->splitparent = bv;
2068: (*split)->N = N;
2069: (*split)->n = n;
2070: (*split)->m = bv->m;
2071: PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2072: }
2073: (*split)->k = bv->k;
2074: (*split)->l = bv->l;
2075: (*split)->nc = bv->nc;
2076: if ((*split)->nc) {
2077: (*split)->ci[0] = -(*split)->nc-1;
2078: (*split)->ci[1] = -(*split)->nc-1;
2079: }
2080: if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2081: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2082: PetscFunctionReturn(PETSC_SUCCESS);
2083: }
2085: /*@
2086: BVGetSplitRows - Splits the `BV` object into two `BV` objects that share the
2087: internal data, using a disjoint horizontal splitting.
2089: Collective
2091: Input Parameters:
2092: + bv - the basis vectors context
2093: . isup - the index set that defines the upper part of the horizontal splitting
2094: - islo - the index set that defines the lower part of the horizontal splitting
2096: Output Parameters:
2097: + U - the resulting `BV` containing the upper rows
2098: - L - the resulting `BV` containing the lower rows
2100: Notes:
2101: The index sets must be such that every MPI process can extract the selected
2102: rows from its local part of the input `BV`, and this part must be contiguous.
2103: With one process, `isup` will list contiguous indices starting from 0, and `islo`
2104: will contain the remaining indices, hence we refer to upper and lower part.
2105: However, with several processes the indices will be interleaved because
2106: `isup` will refer to the upper part of the local array.
2108: The intended use of this function is with matrices of `MATNEST` type, where
2109: `MatNestGetISs()` will return the appropriate index sets.
2111: The returned `BV`s must be seen as references (not copies) of the input
2112: `BV`, that is, modifying them will change the entries of `bv` as well.
2113: The returned `BV`s must not be destroyed. `BVRestoreSplitRows()` must be called
2114: when they are no longer needed.
2116: Pass `NULL` for any of the output `BV`s that is not needed.
2118: Level: advanced
2120: .seealso: [](sec:bv), `BVRestoreSplitRows()`, `BVGetSplit()`
2121: @*/
2122: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2123: {
2124: PetscFunctionBegin;
2127: BVCheckSizes(bv,1);
2130: PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2131: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2132: if (U) {
2133: PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2134: *U = bv->R;
2135: }
2136: if (L) {
2137: PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2138: *L = bv->L;
2139: }
2140: PetscFunctionReturn(PETSC_SUCCESS);
2141: }
2143: /*@
2144: BVRestoreSplitRows - Restore the BV objects obtained with `BVGetSplitRows()`.
2146: Logically Collective
2148: Input Parameters:
2149: + bv - the basis vectors context
2150: . isup - the index set that defines the upper part of the horizontal splitting
2151: . islo - the index set that defines the lower part of the horizontal splitting
2152: . U - the `BV` containing the upper rows
2153: - L - the `BV` containing the lower rows
2155: Note:
2156: The arguments must match the corresponding call to `BVGetSplitRows()`.
2158: Level: advanced
2160: .seealso: [](sec:bv), `BVGetSplitRows()`
2161: @*/
2162: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2163: {
2164: PetscFunctionBegin;
2167: BVCheckSizes(bv,1);
2168: PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2169: 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()");
2170: 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()");
2171: PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2172: bv->lsplit = 0;
2173: PetscFunctionReturn(PETSC_SUCCESS);
2174: }
2176: /*@
2177: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2178: definite inner product.
2180: Logically Collective
2182: Input Parameters:
2183: + bv - basis vectors
2184: - deftol - the tolerance
2186: Options Database Key:
2187: . -bv_definite_tol \<deftol\> - the tolerance
2189: Notes:
2190: When using a non-standard inner product, see `BVSetMatrix()`, the solver needs
2191: to compute $\sqrt{z^*B z}$ for various vectors $z$. If the inner product has not
2192: been declared indefinite, the value $z^*B z$ must be positive, but due to
2193: rounding error a tiny value may become negative. A tolerance is used to
2194: detect this situation. Likewise, in complex arithmetic $z^*B z$ should be
2195: real, and we use the same tolerance to check whether a nonzero imaginary part
2196: can be considered negligible.
2198: This function sets this tolerance, which defaults to `10*PETSC_MACHINE_EPSILON`.
2199: The default value should be good for most applications.
2201: Level: advanced
2203: .seealso: [](sec:bv), `BVSetMatrix()`
2204: @*/
2205: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2206: {
2207: PetscFunctionBegin;
2210: if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2211: else {
2212: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2213: bv->deftol = deftol;
2214: }
2215: PetscFunctionReturn(PETSC_SUCCESS);
2216: }
2218: /*@
2219: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2220: inner product.
2222: Not Collective
2224: Input Parameter:
2225: . bv - the basis vectors
2227: Output Parameter:
2228: . deftol - the tolerance
2230: Level: advanced
2232: .seealso: [](sec:bv), `BVSetDefiniteTolerance()`
2233: @*/
2234: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2235: {
2236: PetscFunctionBegin;
2238: PetscAssertPointer(deftol,2);
2239: *deftol = bv->deftol;
2240: PetscFunctionReturn(PETSC_SUCCESS);
2241: }
2243: /*@
2244: BVSetLeadingDimension - Set the leading dimension to be used for storing the `BV` data.
2246: Not Collective
2248: Input Parameters:
2249: + bv - basis vectors
2250: - ld - the leading dimension
2252: Notes:
2253: This parameter is relevant for `BVMAT`, though it might be employed in other types
2254: as well.
2256: When the internal data of the `BV` is stored as a dense matrix, the leading dimension
2257: has the same meaning as in `MatDenseSetLDA()`, i.e., the distance in number of
2258: elements from one entry of the matrix to the one in the next column at the same
2259: row. The leading dimension refers to the local array, and hence can be different
2260: in different processes.
2262: The user does not need to change this parameter. The default value is equal to the
2263: number of local rows, but this value may be increased a little to guarantee alignment
2264: (especially in the case of GPU storage).
2266: Level: advanced
2268: .seealso: [](sec:bv), `BVGetLeadingDimension()`
2269: @*/
2270: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2271: {
2272: PetscFunctionBegin;
2275: 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");
2276: if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
2277: else {
2278: PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2279: bv->ld = ld;
2280: }
2281: PetscFunctionReturn(PETSC_SUCCESS);
2282: }
2284: /*@
2285: BVGetLeadingDimension - Returns the leading dimension of the `BV`.
2287: Not Collective
2289: Input Parameter:
2290: . bv - the basis vectors
2292: Output Parameter:
2293: . ld - the leading dimension
2295: Level: advanced
2297: Notes:
2298: The returned value may be different in different processes.
2300: The leading dimension must be used when accessing the internal array via
2301: `BVGetArray()` or `BVGetArrayRead()`.
2303: .seealso: [](sec:bv), `BVSetLeadingDimension()`, `BVGetArray()`, `BVGetArrayRead()`
2304: @*/
2305: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2306: {
2307: PetscFunctionBegin;
2309: PetscAssertPointer(ld,2);
2310: *ld = bv->ld;
2311: PetscFunctionReturn(PETSC_SUCCESS);
2312: }