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: Note:
78: `type` should not be retained for later use as it will be an invalid pointer
79: if the `BVType` of `bv` is changed.
81: Level: intermediate
83: .seealso: [](sec:bv), `BVSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
84: @*/
85: PetscErrorCode BVGetType(BV bv,BVType *type)
86: {
87: PetscFunctionBegin;
89: PetscAssertPointer(type,2);
90: *type = ((PetscObject)bv)->type_name;
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*@
95: BVSetSizes - Sets the local and global sizes, and the number of columns.
97: Collective
99: Input Parameters:
100: + bv - the basis vectors
101: . n - the local size (or `PETSC_DECIDE` to have it set)
102: . N - the global size (or `PETSC_DECIDE`)
103: - m - the number of columns
105: Notes:
106: `n` and `N` cannot be both `PETSC_DECIDE`.
107: If one process calls this with `N` equal to `PETSC_DECIDE` then all processes must,
108: otherwise the program will hang.
110: Level: beginner
112: .seealso: [](sec:bv), `BVSetSizesFromVec()`, `BVGetSizes()`, `BVResize()`
113: @*/
114: PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)
115: {
116: PetscInt ma;
117: PetscMPIInt size;
119: PetscFunctionBegin;
123: 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);
124: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
125: 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);
126: 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);
127: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
128: bv->n = n;
129: bv->N = N;
130: bv->m = m;
131: bv->k = m;
132: /* create layout and get actual dimensions */
133: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)bv),&bv->map));
134: PetscCall(PetscLayoutSetSize(bv->map,bv->N));
135: PetscCall(PetscLayoutSetLocalSize(bv->map,bv->n));
136: PetscCall(PetscLayoutSetUp(bv->map));
137: PetscCall(PetscLayoutGetSize(bv->map,&bv->N));
138: PetscCall(PetscLayoutGetLocalSize(bv->map,&bv->n));
139: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
140: PetscCall(BVSetVecType(bv,(size==1)?VECSEQ:VECMPI));
141: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
142: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
143: 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);
144: }
145: if (bv->ops->create) {
146: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
147: PetscUseTypeMethod(bv,create);
148: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
149: bv->ops->create = NULL;
150: bv->defersfo = PETSC_FALSE;
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*@
156: BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
157: Local and global sizes are specified indirectly by passing a template vector.
159: Collective
161: Input Parameters:
162: + bv - the basis vectors
163: . t - the template vector
164: - m - the number of columns
166: Notes:
167: Apart from the local and global sizes, the template vector `t` is also used
168: to get the vector type, which will be set as the `BV` vector type with
169: `BVSetVecType()`. This is relevant in cases where computation must be done
170: on the GPU.
172: If `t` is a nested vector (type `VECNEST`), then the vector type will be
173: obtained from the first subvector. The rationale is that `BV` will operate
174: internally with regular vectors, even though the user problem is defined via
175: nested matrices and vectors.
177: Level: beginner
179: .seealso: [](sec:bv), `BVSetSizes()`, `BVGetSizes()`, `BVResize()`, `BVSetVecType()`
180: @*/
181: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
182: {
183: PetscInt ma;
184: PetscLayout map;
185: VecType vtype;
186: PetscBool nest;
187: Vec v;
189: PetscFunctionBegin;
192: PetscCheckSameComm(bv,1,t,2);
194: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
195: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
196: PetscCall(PetscObjectTypeCompare((PetscObject)t,VECNEST,&nest));
197: if (nest) {
198: PetscCall(VecNestGetSubVec(t,0,&v));
199: PetscCall(VecGetType(v,&vtype));
200: } else PetscCall(VecGetType(t,&vtype));
201: PetscCall(BVSetVecType(bv,vtype));
202: PetscCall(VecGetLayout(t,&map));
203: PetscCall(PetscLayoutReference(map,&bv->map));
204: PetscCall(VecGetSize(t,&bv->N));
205: PetscCall(VecGetLocalSize(t,&bv->n));
206: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
207: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
208: 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);
209: }
210: bv->m = m;
211: bv->k = m;
212: if (bv->ops->create) {
213: PetscUseTypeMethod(bv,create);
214: bv->ops->create = NULL;
215: bv->defersfo = PETSC_FALSE;
216: }
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@
221: BVGetSizes - Returns the local and global sizes, and the number of columns.
223: Not Collective
225: Input Parameter:
226: . bv - the basis vectors
228: Output Parameters:
229: + n - the local size
230: . N - the global size
231: - m - the number of columns
233: Note:
234: Normal usage requires that `bv` has already been given its sizes, otherwise
235: the call fails. However, this function can also be used to determine if
236: a `BV` object has been initialized completely (sizes and type). For this,
237: call with `n=NULL` and `N=NULL`, then a return value of `m=0` indicates that
238: the `BV` object is not ready for use yet.
240: Level: beginner
242: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
243: @*/
244: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
245: {
246: PetscFunctionBegin;
247: if (!bv) {
248: if (m && !n && !N) *m = 0;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: if (n || N) BVCheckSizes(bv,1);
253: if (n) *n = bv->n;
254: if (N) *N = bv->N;
255: if (m) *m = bv->m;
256: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@
261: BVSetNumConstraints - Set the number of constraints.
263: Logically Collective
265: Input Parameters:
266: + V - basis vectors
267: - nc - number of constraints
269: Notes:
270: This function sets the number of constraints to `nc` and marks all remaining
271: columns as regular. Normal usage would be to call `BVInsertConstraints()` instead.
273: If `nc` is smaller than the previously set value, then some of the constraints
274: are discarded. In particular, using `nc=0` removes all constraints preserving
275: the content of regular columns.
277: Level: developer
279: .seealso: [](sec:bv), `BVInsertConstraints()`
280: @*/
281: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
282: {
283: PetscInt total,diff,i;
284: Vec x,y;
286: PetscFunctionBegin;
289: PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
291: BVCheckSizes(V,1);
292: PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
294: diff = nc-V->nc;
295: if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
296: total = V->nc+V->m;
297: PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
298: if (diff<0) { /* lessen constraints, shift contents of BV */
299: for (i=0;i<V->m;i++) {
300: PetscCall(BVGetColumn(V,i,&x));
301: PetscCall(BVGetColumn(V,i+diff,&y));
302: PetscCall(VecCopy(x,y));
303: PetscCall(BVRestoreColumn(V,i,&x));
304: PetscCall(BVRestoreColumn(V,i+diff,&y));
305: }
306: }
307: V->nc = nc;
308: V->ci[0] = -V->nc-1;
309: V->ci[1] = -V->nc-1;
310: V->m = total-nc;
311: V->l = PetscMin(V->l,V->m);
312: V->k = PetscMin(V->k,V->m);
313: PetscCall(PetscObjectStateIncrease((PetscObject)V));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@
318: BVGetNumConstraints - Returns the number of constraints.
320: Not Collective
322: Input Parameter:
323: . bv - the basis vectors
325: Output Parameter:
326: . nc - the number of constraints
328: Level: advanced
330: .seealso: [](sec:bv), `BVGetSizes()`, `BVInsertConstraints()`
331: @*/
332: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
333: {
334: PetscFunctionBegin;
336: PetscAssertPointer(nc,2);
337: *nc = bv->nc;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: BVResize - Change the number of columns.
344: Collective
346: Input Parameters:
347: + bv - the basis vectors
348: . m - the new number of columns
349: - copy - a flag indicating whether current values should be kept
351: Note:
352: Internal storage is reallocated. If the `copy` flag is set to true, then
353: the contents are copied to the leading part of the new space.
355: Level: advanced
357: .seealso: [](sec:bv), `BVSetSizes()`, `BVSetSizesFromVec()`
358: @*/
359: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
360: {
361: PetscScalar *array;
362: const PetscScalar *omega;
363: Vec v;
365: PetscFunctionBegin;
370: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
371: PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
372: if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
373: BVCheckOp(bv,1,resize);
375: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
376: PetscUseTypeMethod(bv,resize,m,copy);
377: PetscCall(VecDestroy(&bv->buffer));
378: PetscCall(BVDestroy(&bv->cached));
379: PetscCall(PetscFree2(bv->h,bv->c));
380: if (bv->omega) {
381: if (bv->cuda) {
382: #if defined(PETSC_HAVE_CUDA)
383: PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
384: #else
385: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
386: #endif
387: } else if (bv->hip) {
388: #if defined(PETSC_HAVE_HIP)
389: PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
390: #else
391: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
392: #endif
393: } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
394: if (copy) {
395: PetscCall(VecGetArray(v,&array));
396: PetscCall(VecGetArrayRead(bv->omega,&omega));
397: PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
398: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
399: PetscCall(VecRestoreArray(v,&array));
400: } else PetscCall(VecSet(v,1.0));
401: PetscCall(VecDestroy(&bv->omega));
402: bv->omega = v;
403: }
404: bv->m = m;
405: bv->k = PetscMin(bv->k,m);
406: bv->l = PetscMin(bv->l,m);
407: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
408: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*@
413: BVSetActiveColumns - Specify the columns that will be involved in operations.
415: Logically Collective
417: Input Parameters:
418: + bv - the basis vectors context
419: . l - number of leading columns
420: - k - number of active columns
422: Notes:
423: In operations such as `BVMult()` or `BVDot()`, only the first `k` columns are
424: considered. This is useful when the `BV` is filled from left to right, so
425: the last `m-k` columns do not have relevant information.
427: Also in operations such as `BVMult()` or `BVDot()`, the first `l` columns are
428: normally not included in the computation. See the manual page of each
429: operation.
431: In orthogonalization operations, the first `l` columns are treated
432: differently, they participate in the orthogonalization but the computed
433: coefficients are not stored.
435: Use `PETSC_CURRENT` to leave any of the values unchanged. Use `PETSC_DETERMINE`
436: to set `l` to the minimum value (`0`) and `k` to the maximum (`m`).
438: Level: intermediate
440: .seealso: [](sec:bv), `BVGetActiveColumns()`, `BVSetSizes()`
441: @*/
442: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
443: {
444: PetscFunctionBegin;
448: BVCheckSizes(bv,1);
449: if (PetscUnlikely(k == PETSC_DETERMINE)) {
450: bv->k = bv->m;
451: } else if (k != PETSC_CURRENT) {
452: 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);
453: bv->k = k;
454: }
455: if (PetscUnlikely(l == PETSC_DETERMINE)) {
456: bv->l = 0;
457: } else if (l != PETSC_CURRENT) {
458: 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);
459: bv->l = l;
460: }
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: BVGetActiveColumns - Returns the current active dimensions.
467: Not Collective
469: Input Parameter:
470: . bv - the basis vectors context
472: Output Parameters:
473: + l - number of leading columns
474: - k - number of active columns
476: Level: intermediate
478: .seealso: [](sec:bv), `BVSetActiveColumns()`
479: @*/
480: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
481: {
482: PetscFunctionBegin;
484: if (l) *l = bv->l;
485: if (k) *k = bv->k;
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /*@
490: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
492: Collective
494: Input Parameters:
495: + bv - the basis vectors context
496: . B - a Hermitian matrix (may be `NULL`)
497: - indef - a flag indicating if the matrix is indefinite
499: Notes:
500: This is used to specify a non-standard inner product, whose matrix
501: representation is given by `B`. Then, all inner products required during
502: orthogonalization are computed as $(x,y)_B=y^*Bx$ rather than the
503: standard form $(x,y)=y^*x$.
505: Matrix `B` must be real symmetric (or complex Hermitian). A genuine inner
506: product requires that `B` is also positive (semi-)definite. However, we
507: also allow for an indefinite `B` (setting `indef=PETSC_TRUE`), in which
508: case the orthogonalization uses an indefinite inner product.
510: This affects operations `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, and variants.
512: Setting `B=NULL` has the same effect as if the identity matrix was passed.
514: Level: advanced
516: .seealso: [](sec:bv), `BVGetMatrix()`, `BVDot()`, `BVNorm()`, `BVOrthogonalize()`, `BVSetDefiniteTolerance()`
517: @*/
518: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
519: {
520: PetscInt m,n;
522: PetscFunctionBegin;
525: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
526: if (B) {
528: PetscCall(MatGetLocalSize(B,&m,&n));
529: PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
530: 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);
531: }
532: PetscCall(PetscObjectReference((PetscObject)B));
533: PetscCall(MatDestroy(&bv->matrix));
534: bv->matrix = B;
535: bv->indef = indef;
536: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
537: if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
538: if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
539: }
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: /*@
544: BVGetMatrix - Retrieves the matrix representation of the inner product.
546: Not Collective
548: Input Parameter:
549: . bv - the basis vectors context
551: Output Parameters:
552: + B - the matrix of the inner product (may be `NULL`)
553: - indef - the flag indicating if the matrix is indefinite
555: Level: advanced
557: .seealso: [](sec:bv), `BVSetMatrix()`
558: @*/
559: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
560: {
561: PetscFunctionBegin;
563: if (B) *B = bv->matrix;
564: if (indef) *indef = bv->indef;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@
569: BVApplyMatrix - Multiplies a vector by the matrix representation of the
570: inner product.
572: Neighbor-wise Collective
574: Input Parameters:
575: + bv - the basis vectors context
576: - x - the vector
578: Output Parameter:
579: . y - the result
581: Note:
582: If no matrix was specified this function copies the vector.
584: Level: advanced
586: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrixBV()`
587: @*/
588: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
589: {
590: PetscFunctionBegin;
594: if (bv->matrix) {
595: PetscCall(BV_IPMatMult(bv,x));
596: PetscCall(VecCopy(bv->Bx,y));
597: } else PetscCall(VecCopy(x,y));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*@
602: BVApplyMatrixBV - Multiplies the `BV` vectors by the matrix representation
603: of the inner product.
605: Neighbor-wise Collective
607: Input Parameter:
608: . X - the basis vectors context
610: Output Parameter:
611: . Y - the basis vectors to store the result (optional)
613: Note:
614: This function computes $Y = B X$, where $B$ is the matrix given with
615: `BVSetMatrix()`. This operation is computed as in `BVMatMult()`.
616: If no matrix was specified, then it just copies $Y = X$.
618: If no `Y` is given, the result is stored internally in the cached `BV`.
620: Level: developer
622: .seealso: [](sec:bv), `BVSetMatrix()`, `BVApplyMatrix()`, `BVMatMult()`, `BVGetCachedBV()`
623: @*/
624: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
625: {
626: PetscFunctionBegin;
628: if (Y) {
630: if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
631: else PetscCall(BVCopy(X,Y));
632: } else PetscCall(BV_IPMatMultBV(X));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@
637: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
639: Logically Collective
641: Input Parameters:
642: + bv - the basis vectors context
643: - omega - a vector representing the diagonal of the signature matrix
645: Note:
646: The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.
648: Level: developer
650: .seealso: [](sec:bv), `BVSetMatrix()`, `BVGetSignature()`
651: @*/
652: PetscErrorCode BVSetSignature(BV bv,Vec omega)
653: {
654: PetscInt i,n;
655: const PetscScalar *pomega;
656: PetscScalar *intern;
658: PetscFunctionBegin;
660: BVCheckSizes(bv,1);
664: PetscCall(VecGetSize(omega,&n));
665: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
666: PetscCall(BV_AllocateSignature(bv));
667: if (bv->indef) {
668: PetscCall(VecGetArrayRead(omega,&pomega));
669: PetscCall(VecGetArray(bv->omega,&intern));
670: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
671: PetscCall(VecRestoreArray(bv->omega,&intern));
672: PetscCall(VecRestoreArrayRead(omega,&pomega));
673: } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
674: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: /*@
679: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
681: Not Collective
683: Input Parameter:
684: . bv - the basis vectors context
686: Output Parameter:
687: . omega - a vector representing the diagonal of the signature matrix
689: Note:
690: The signature matrix $\Omega = V^*B V$ is relevant only for an indefinite $B$.
692: Level: developer
694: .seealso: [](sec:bv), `BVSetMatrix()`, `BVSetSignature()`
695: @*/
696: PetscErrorCode BVGetSignature(BV bv,Vec omega)
697: {
698: PetscInt i,n;
699: PetscScalar *pomega;
700: const PetscScalar *intern;
702: PetscFunctionBegin;
704: BVCheckSizes(bv,1);
708: PetscCall(VecGetSize(omega,&n));
709: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
710: if (bv->indef && bv->omega) {
711: PetscCall(VecGetArray(omega,&pomega));
712: PetscCall(VecGetArrayRead(bv->omega,&intern));
713: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
714: PetscCall(VecRestoreArrayRead(bv->omega,&intern));
715: PetscCall(VecRestoreArray(omega,&pomega));
716: } else PetscCall(VecSet(omega,1.0));
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: BVSetBufferVec - Attach a vector object to be used as buffer space for
722: several operations.
724: Collective
726: Input Parameters:
727: + bv - the basis vectors context)
728: - buffer - the vector
730: Notes:
731: Use `BVGetBufferVec()` to retrieve the vector (for example, to free it
732: at the end of the computations).
734: The vector must be sequential of length $(n_c+m)m$, where $m$ is the number
735: of columns of `bv` and $n_c$ is the number of constraints.
737: Level: developer
739: .seealso: [](sec:bv), `BVGetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`
740: @*/
741: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
742: {
743: PetscInt ld,n;
744: PetscMPIInt size;
746: PetscFunctionBegin;
749: BVCheckSizes(bv,1);
750: PetscCall(VecGetSize(buffer,&n));
751: ld = bv->m+bv->nc;
752: PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
753: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
754: PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
756: PetscCall(PetscObjectReference((PetscObject)buffer));
757: PetscCall(VecDestroy(&bv->buffer));
758: bv->buffer = buffer;
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: /*@
763: BVGetBufferVec - Obtain the buffer vector associated with the `BV` object.
765: Collective
767: Input Parameter:
768: . bv - the basis vectors context
770: Output Parameter:
771: . buffer - vector
773: Notes:
774: The vector is created if not available previously. It is a sequential vector
775: of length $(n_c+m) m$, where $m$ is the number of columns of `bv` and $n_c$ is the number
776: of constraints.
778: Developer Notes:
779: The buffer vector is viewed as a column-major matrix with leading dimension equal
780: to $n_c+m$, and $m$ columns at most. In the most common usage, it has the structure
781: $$
782: \left[\begin{array}{c|c}
783: s & \begin{array}{c} C \\\hline H \end{array}
784: \end{array}\right]
785: $$
786: where $H$ is an upper Hessenberg matrix of order $m (m-1)$, $C$ contains coefficients
787: related to orthogonalization against constraints (first $n_c$ rows), and $s$ is the
788: first column that contains scratch values computed during Gram-Schmidt
789: orthogonalization. In particular, `BVDotColumn()` and `BVMultColumn()` use $s$ to
790: store the coefficients.
792: Level: developer
794: .seealso: [](sec:bv), `BVSetBufferVec()`, `BVSetSizes()`, `BVGetNumConstraints()`, `BVDotColumn()`, `BVMultColumn()`
795: @*/
796: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
797: {
798: PetscInt ld;
800: PetscFunctionBegin;
802: PetscAssertPointer(buffer,2);
803: BVCheckSizes(bv,1);
804: if (!bv->buffer) {
805: ld = bv->m+bv->nc;
806: PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
807: PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
808: PetscCall(VecSetType(bv->buffer,bv->vtype));
809: }
810: *buffer = bv->buffer;
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: /*@
815: BVSetRandomContext - Sets the `PetscRandom` object associated with the `BV`,
816: to be used in operations that need random numbers.
818: Collective
820: Input Parameters:
821: + bv - the basis vectors context
822: - rand - the random number generator context
824: Level: advanced
826: .seealso: [](sec:bv), `BVGetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`, `PetscRandomCreate()`
827: @*/
828: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
829: {
830: PetscFunctionBegin;
833: PetscCheckSameComm(bv,1,rand,2);
834: PetscCall(PetscObjectReference((PetscObject)rand));
835: PetscCall(PetscRandomDestroy(&bv->rand));
836: bv->rand = rand;
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: /*@
841: BVGetRandomContext - Gets the `PetscRandom` object associated with the `BV`.
843: Collective
845: Input Parameter:
846: . bv - the basis vectors context
848: Output Parameter:
849: . rand - the random number generator context
851: Level: advanced
853: .seealso: [](sec:bv), `BVSetRandomContext()`, `BVSetRandom()`, `BVSetRandomNormal()`, `BVSetRandomColumn()`, `BVSetRandomCond()`
854: @*/
855: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
856: {
857: PetscFunctionBegin;
859: PetscAssertPointer(rand,2);
860: if (!bv->rand) {
861: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
862: if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
863: if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
864: if (bv->rrandom) {
865: PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
866: PetscCall(PetscRandomSeed(bv->rand));
867: }
868: }
869: *rand = bv->rand;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: BVSetFromOptions - Sets `BV` options from the options database.
876: Collective
878: Input Parameter:
879: . bv - the basis vectors context
881: Note:
882: To see all options, run your program with the `-help` option.
884: Level: beginner
886: .seealso: [](sec:bv), `BVSetOptionsPrefix()`
887: @*/
888: PetscErrorCode BVSetFromOptions(BV bv)
889: {
890: char type[256];
891: PetscBool flg1,flg2,flg3,flg4;
892: PetscReal r;
893: BVOrthogType otype;
894: BVOrthogRefineType orefine;
895: BVOrthogBlockType oblock;
897: PetscFunctionBegin;
899: PetscCall(BVRegisterAll());
900: PetscObjectOptionsBegin((PetscObject)bv);
901: PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
902: if (flg1) PetscCall(BVSetType(bv,type));
903: else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));
905: otype = bv->orthog_type;
906: PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
907: orefine = bv->orthog_ref;
908: PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
909: r = bv->orthog_eta;
910: PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
911: oblock = bv->orthog_block;
912: PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
913: if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
915: PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
917: PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
918: if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
920: /* undocumented option to generate random vectors that are independent of the number of processes */
921: PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
923: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
924: else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
925: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
926: PetscOptionsEnd();
927: bv->sfocalled = PETSC_TRUE;
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: /*@
932: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
933: vectors (classical or modified Gram-Schmidt with or without refinement), and
934: for the block-orthogonalization (simultaneous orthogonalization of a set of
935: vectors).
937: Logically Collective
939: Input Parameters:
940: + bv - the basis vectors context
941: . type - the method of vector orthogonalization
942: . refine - type of refinement
943: . eta - parameter for selective refinement
944: - block - the method of block orthogonalization
946: Options Database Keys:
947: + -bv_orthog_type (cgs|mgs) - the vector orthogonalization
948: . -bv_orthog_refine (never|ifneeded|always) - the refinement type, default is `ifneeded`
949: . -bv_orthog_eta eta - the value of `eta`
950: - -bv_orthog_block (gs|chol|tsqr|tsqrchol|svqb) - the block-orthogonalization method
952: Notes:
953: The default settings work well for most problems.
955: The parameter `eta` should be a real value between 0 and 1, that is used only when
956: the refinement type is `ifneeded`. Use `PETSC_DETERMINE` to set a reasonable
957: default value. Use `PETSC_CURRENT` to leave the current value unchanged.
959: When using several processes, `BV_ORTHOG_MGS` is likely to result in bad scalability.
961: If the method set for block orthogonalization is `BV_ORTHOG_BLOCK_GS`, then the computation
962: is done column by column with the vector orthogonalization.
964: Level: advanced
966: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVGetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
967: @*/
968: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
969: {
970: PetscFunctionBegin;
976: switch (type) {
977: case BV_ORTHOG_CGS:
978: case BV_ORTHOG_MGS:
979: bv->orthog_type = type;
980: break;
981: default:
982: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
983: }
984: switch (refine) {
985: case BV_ORTHOG_REFINE_NEVER:
986: case BV_ORTHOG_REFINE_IFNEEDED:
987: case BV_ORTHOG_REFINE_ALWAYS:
988: bv->orthog_ref = refine;
989: break;
990: default:
991: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
992: }
993: if (eta == (PetscReal)PETSC_DETERMINE) {
994: bv->orthog_eta = 0.7071;
995: } else if (eta != (PetscReal)PETSC_CURRENT) {
996: PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
997: bv->orthog_eta = eta;
998: }
999: switch (block) {
1000: case BV_ORTHOG_BLOCK_GS:
1001: case BV_ORTHOG_BLOCK_CHOL:
1002: case BV_ORTHOG_BLOCK_TSQR:
1003: case BV_ORTHOG_BLOCK_TSQRCHOL:
1004: case BV_ORTHOG_BLOCK_SVQB:
1005: bv->orthog_block = block;
1006: break;
1007: default:
1008: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
1009: }
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /*@
1014: BVGetOrthogonalization - Gets the orthogonalization settings from the `BV` object.
1016: Not Collective
1018: Input Parameter:
1019: . bv - basis vectors context
1021: Output Parameters:
1022: + type - the method of vector orthogonalization
1023: . refine - type of refinement
1024: . eta - parameter for selective refinement
1025: - block - the method of block orthogonalization
1027: Level: advanced
1029: .seealso: [](sec:bv), `BVOrthogonalizeColumn()`, `BVSetOrthogonalization()`, `BVOrthogType`, `BVOrthogRefineType`, `BVOrthogBlockType`
1030: @*/
1031: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1032: {
1033: PetscFunctionBegin;
1035: if (type) *type = bv->orthog_type;
1036: if (refine) *refine = bv->orthog_ref;
1037: if (eta) *eta = bv->orthog_eta;
1038: if (block) *block = bv->orthog_block;
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: /*@
1043: BVSetMatMultMethod - Specifies the method used for the `BVMatMult()` operation.
1045: Logically Collective
1047: Input Parameters:
1048: + bv - the basis vectors context
1049: - method - the method for the `BVMatMult()` operation
1051: Options Database Key:
1052: . -bv_matmult (vecs|mat) - choose the method
1054: Note:
1055: The default is `BV_MATMULT_MAT` except in the case of `BVVECS`.
1057: Level: advanced
1059: .seealso: [](sec:bv), `BVMatMult()`, `BVGetMatMultMethod()`, `BVMatMultType`
1060: @*/
1061: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1062: {
1063: PetscFunctionBegin;
1066: switch (method) {
1067: case BV_MATMULT_VECS:
1068: case BV_MATMULT_MAT:
1069: bv->vmm = method;
1070: break;
1071: case BV_MATMULT_MAT_SAVE:
1072: PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1073: bv->vmm = BV_MATMULT_MAT;
1074: break;
1075: default:
1076: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1077: }
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }
1081: /*@
1082: BVGetMatMultMethod - Gets the method used for the `BVMatMult()` operation.
1084: Not Collective
1086: Input Parameter:
1087: . bv - basis vectors context
1089: Output Parameter:
1090: . method - the method for the `BVMatMult()` operation
1092: Level: advanced
1094: .seealso: [](sec:bv), `BVMatMult()`, `BVSetMatMultMethod()`, `BVMatMultType`
1095: @*/
1096: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1097: {
1098: PetscFunctionBegin;
1100: PetscAssertPointer(method,2);
1101: *method = bv->vmm;
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: /*@
1106: BVGetColumn - Returns a `Vec` object that contains the entries of the
1107: requested column of the basis vectors object.
1109: Logically Collective
1111: Input Parameters:
1112: + bv - the basis vectors context
1113: - j - the index of the requested column
1115: Output Parameter:
1116: . v - vector containing the jth column
1118: Notes:
1119: The returned `Vec` must be seen as a reference (not a copy) of the `BV`
1120: column, that is, modifying the `Vec` will change the `BV` entries as well.
1122: The returned `Vec` must not be destroyed. `BVRestoreColumn()` must be
1123: called when it is no longer needed. At most, two columns can be fetched,
1124: that is, this function can only be called twice before the corresponding
1125: `BVRestoreColumn()` is invoked.
1127: A negative index `j` selects the `i`-th constraint, where `i=-j`. Constraints
1128: should not be modified.
1130: Level: beginner
1132: .seealso: [](sec:bv), `BVRestoreColumn()`, `BVInsertConstraints()`
1133: @*/
1134: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1135: {
1136: PetscInt l;
1138: PetscFunctionBegin;
1141: BVCheckSizes(bv,1);
1142: BVCheckOp(bv,1,getcolumn);
1144: 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);
1145: 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);
1146: 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);
1147: l = BVAvailableVec;
1148: PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1149: PetscUseTypeMethod(bv,getcolumn,j,v);
1150: bv->ci[l] = j;
1151: PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
1152: PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1153: *v = bv->cv[l];
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: /*@
1158: BVRestoreColumn - Restore a column obtained with `BVGetColumn()`.
1160: Logically Collective
1162: Input Parameters:
1163: + bv - the basis vectors context
1164: . j - the index of the column
1165: - v - vector obtained with `BVGetColumn()`
1167: Note:
1168: The arguments must match the corresponding call to `BVGetColumn()`.
1170: Level: beginner
1172: .seealso: [](sec:bv), `BVGetColumn()`
1173: @*/
1174: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1175: {
1176: PetscObjectId id;
1177: PetscObjectState st;
1178: PetscInt l;
1180: PetscFunctionBegin;
1183: BVCheckSizes(bv,1);
1185: PetscAssertPointer(v,3);
1187: 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);
1188: 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);
1189: 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);
1190: l = (j==bv->ci[0])? 0: 1;
1191: PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1192: PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1193: PetscCall(VecGetState(*v,&st));
1194: if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1195: PetscUseTypeMethod(bv,restorecolumn,j,v);
1196: bv->ci[l] = -bv->nc-1;
1197: bv->st[l] = -1;
1198: bv->id[l] = 0;
1199: *v = NULL;
1200: PetscFunctionReturn(PETSC_SUCCESS);
1201: }
1203: /*@C
1204: BVGetArray - Returns a pointer to a contiguous array that contains this
1205: process' portion of the `BV` data.
1207: Logically Collective
1209: Input Parameter:
1210: . bv - the basis vectors context
1212: Output Parameter:
1213: . a - location to put pointer to the array
1215: Notes:
1216: `BVRestoreArray()` must be called when access to the array is no longer needed.
1217: This operation may imply a data copy, for `BV` types that do not store
1218: data contiguously in memory.
1220: The pointer will normally point to the first entry of the first column,
1221: but if the `BV` has constraints then these go before the regular columns.
1223: Note that for manipulating the pointer to the `BV` array, one must take into
1224: account the leading dimension, which might be different from the local
1225: number of rows, see `BVGetLeadingDimension()`.
1227: Use `BVGetArrayRead()` for read-only access.
1229: Level: advanced
1231: .seealso: [](sec:bv), `BVRestoreArray()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArrayRead()`, `BVType`
1232: @*/
1233: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1234: {
1235: PetscFunctionBegin;
1238: BVCheckSizes(bv,1);
1239: BVCheckOp(bv,1,getarray);
1240: PetscUseTypeMethod(bv,getarray,a);
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: /*@C
1245: BVRestoreArray - Restore the `BV` object after `BVGetArray()` has been called.
1247: Logically Collective
1249: Input Parameters:
1250: + bv - the basis vectors context
1251: - a - location of pointer to array obtained from `BVGetArray()`
1253: Note:
1254: This operation may imply a data copy, for `BV` types that do not store
1255: data contiguously in memory.
1257: Level: advanced
1259: .seealso: [](sec:bv), `BVGetArray()`
1260: @*/
1261: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1262: {
1263: PetscFunctionBegin;
1266: BVCheckSizes(bv,1);
1267: PetscTryTypeMethod(bv,restorearray,a);
1268: if (a) *a = NULL;
1269: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }
1273: /*@C
1274: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1275: contains this processor's portion of the `BV` data.
1277: Not Collective
1279: Input Parameter:
1280: . bv - the basis vectors context
1282: Output Parameter:
1283: . a - location to put pointer to the array
1285: Notes:
1286: `BVRestoreArrayRead()` must be called when access to the array is no
1287: longer needed. This operation may imply a data copy, for `BV` types that
1288: do not store data contiguously in memory.
1290: The pointer will normally point to the first entry of the first column,
1291: but if the `BV` has constraints then these go before the regular columns.
1293: Level: advanced
1295: .seealso: [](sec:bv), `BVRestoreArrayRead()`, `BVInsertConstraints()`, `BVGetLeadingDimension()`, `BVGetArray()`, `BVType`
1296: @*/
1297: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1298: {
1299: PetscFunctionBegin;
1302: BVCheckSizes(bv,1);
1303: BVCheckOp(bv,1,getarrayread);
1304: PetscUseTypeMethod(bv,getarrayread,a);
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: /*@C
1309: BVRestoreArrayRead - Restore the `BV` object after `BVGetArrayRead()` has
1310: been called.
1312: Not Collective
1314: Input Parameters:
1315: + bv - the basis vectors context
1316: - a - location of pointer to array obtained from `BVGetArrayRead()`
1318: Level: advanced
1320: .seealso: [](sec:bv), `BVGetArrayRead()`
1321: @*/
1322: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1323: {
1324: PetscFunctionBegin;
1327: BVCheckSizes(bv,1);
1328: PetscTryTypeMethod(bv,restorearrayread,a);
1329: if (a) *a = NULL;
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }
1333: /*@
1334: BVCreateVec - Creates a new `Vec` object with the same type and dimensions
1335: as the columns of the basis vectors object.
1337: Collective
1339: Input Parameter:
1340: . bv - the basis vectors context
1342: Output Parameter:
1343: . v - the new vector
1345: Note:
1346: The user is responsible for destroying the returned vector.
1348: Level: beginner
1350: .seealso: [](sec:bv), `BVCreateMat()`, `BVCreateVecEmpty()`
1351: @*/
1352: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1353: {
1354: PetscFunctionBegin;
1356: BVCheckSizes(bv,1);
1357: PetscAssertPointer(v,2);
1358: PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1359: PetscCall(VecSetLayout(*v,bv->map));
1360: PetscCall(VecSetType(*v,bv->vtype));
1361: PetscCall(VecSetUp(*v));
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }
1365: /*@
1366: BVCreateVecEmpty - Creates a new `Vec` object with the same type and dimensions
1367: as the columns of the basis vectors object, but without internal array.
1369: Collective
1371: Input Parameter:
1372: . bv - the basis vectors context
1374: Output Parameter:
1375: . v - the new vector
1377: Note:
1378: This works as `BVCreateVec()`, but the new vector does not have the array allocated,
1379: so the intended usage is with `VecPlaceArray()`.
1381: Level: developer
1383: .seealso: [](sec:bv), `BVCreateVec()`
1384: @*/
1385: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1386: {
1387: PetscBool standard,cuda,hip,mpi;
1388: PetscInt N,nloc,bs;
1390: PetscFunctionBegin;
1392: BVCheckSizes(bv,1);
1393: PetscAssertPointer(v,2);
1395: PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1396: PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1397: PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
1398: if (standard || cuda || hip) {
1399: PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
1400: PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1401: PetscCall(PetscLayoutGetSize(bv->map,&N));
1402: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1403: if (cuda) {
1404: #if defined(PETSC_HAVE_CUDA)
1405: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1406: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1407: #endif
1408: } else if (hip) {
1409: #if defined(PETSC_HAVE_HIP)
1410: if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1411: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1412: #endif
1413: } else {
1414: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1415: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1416: }
1417: } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: /*@
1422: BVSetVecType - Set the vector type to be used when creating vectors via `BVCreateVec()`.
1424: Collective
1426: Input Parameters:
1427: + bv - the basis vectors context
1428: - vtype - the vector type
1430: Level: advanced
1432: Note:
1433: This is not needed if the `BV` object is set up with `BVSetSizesFromVec()`, but may be
1434: required in the case of `BVSetSizes()` if one wants to work with non-standard vectors.
1436: .seealso: [](sec:bv), `BVCreateVec()`, `BVGetVecType()`, `BVSetSizesFromVec()`, `BVSetSizes()`
1437: @*/
1438: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1439: {
1440: PetscBool std;
1441: PetscMPIInt size;
1443: PetscFunctionBegin;
1445: PetscCall(PetscFree(bv->vtype));
1446: PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1447: if (std) {
1448: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1449: PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1450: } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1451: PetscFunctionReturn(PETSC_SUCCESS);
1452: }
1454: /*@
1455: BVGetVecType - Get the vector type to be used when creating vectors via `BVCreateVec()`.
1457: Not Collective
1459: Input Parameter:
1460: . bv - the basis vectors context
1462: Output Parameter:
1463: . vtype - the vector type
1465: Level: advanced
1467: .seealso: [](sec:bv), `BVSetVecType()`, `BVCreateVec()`
1468: @*/
1469: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1470: {
1471: PetscFunctionBegin;
1473: PetscAssertPointer(vtype,2);
1474: *vtype = bv->vtype;
1475: PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1478: /*@
1479: BVCreateMat - Creates a new `Mat` object of dense type and copies the contents
1480: of the `BV` object.
1482: Collective
1484: Input Parameter:
1485: . bv - the basis vectors context
1487: Output Parameter:
1488: . A - the new matrix
1490: Notes:
1491: The user is responsible for destroying the returned matrix.
1493: The matrix contains all columns of the `BV`, not just the active columns.
1495: Level: intermediate
1497: .seealso: [](sec:bv), `BVCreateFromMat()`, `BVCreateVec()`, `BVGetMat()`
1498: @*/
1499: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1500: {
1501: PetscInt ksave,lsave;
1502: Mat B;
1504: PetscFunctionBegin;
1506: BVCheckSizes(bv,1);
1507: PetscAssertPointer(A,2);
1509: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1510: lsave = bv->l;
1511: ksave = bv->k;
1512: bv->l = 0;
1513: bv->k = bv->m;
1514: PetscCall(BVGetMat(bv,&B));
1515: PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1516: PetscCall(BVRestoreMat(bv,&B));
1517: bv->l = lsave;
1518: bv->k = ksave;
1519: PetscFunctionReturn(PETSC_SUCCESS);
1520: }
1522: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1523: {
1524: PetscScalar *vv,*aa;
1525: PetscBool create=PETSC_FALSE;
1526: PetscInt m,cols;
1528: PetscFunctionBegin;
1529: m = bv->k-bv->l;
1530: if (!bv->Aget) create=PETSC_TRUE;
1531: else {
1532: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1533: PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1534: PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1535: if (cols!=m) {
1536: PetscCall(MatDestroy(&bv->Aget));
1537: create=PETSC_TRUE;
1538: }
1539: }
1540: PetscCall(BVGetArray(bv,&vv));
1541: if (create) {
1542: 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 */
1543: PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1544: }
1545: PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld)); /* set the actual pointer */
1546: *A = bv->Aget;
1547: PetscFunctionReturn(PETSC_SUCCESS);
1548: }
1550: /*@
1551: BVGetMat - Returns a `Mat` object of dense type that shares the memory of
1552: the `BV` object.
1554: Collective
1556: Input Parameter:
1557: . bv - the basis vectors context
1559: Output Parameter:
1560: . A - the matrix
1562: Notes:
1563: The returned matrix contains only the active columns. If the content of
1564: the `Mat` is modified, these changes are also done in the `BV` object. The
1565: user must call `BVRestoreMat()` when no longer needed.
1567: This operation implies a call to `BVGetArray()`, which may result in data
1568: copies.
1570: Level: advanced
1572: .seealso: [](sec:bv), `BVRestoreMat()`, `BVCreateMat()`, `BVGetArray()`
1573: @*/
1574: PetscErrorCode BVGetMat(BV bv,Mat *A)
1575: {
1576: char name[64];
1578: PetscFunctionBegin;
1580: BVCheckSizes(bv,1);
1581: PetscAssertPointer(A,2);
1582: PetscUseTypeMethod(bv,getmat,A);
1583: if (((PetscObject)bv)->name) { /* set A's name based on BV name */
1584: PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1585: PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1586: PetscCall(PetscObjectSetName((PetscObject)*A,name));
1587: }
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1592: {
1593: PetscScalar *vv,*aa;
1595: PetscFunctionBegin;
1596: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1597: vv = aa-(bv->nc+bv->l)*bv->ld;
1598: PetscCall(MatDenseResetArray(bv->Aget));
1599: PetscCall(BVRestoreArray(bv,&vv));
1600: *A = NULL;
1601: PetscFunctionReturn(PETSC_SUCCESS);
1602: }
1604: /*@
1605: BVRestoreMat - Restores the Mat obtained with `BVGetMat()`.
1607: Logically Collective
1609: Input Parameters:
1610: + bv - the basis vectors context
1611: - A - the fetched matrix
1613: Note:
1614: A call to this function must match a previous call of `BVGetMat()`.
1615: The effect is that the contents of the `Mat` are copied back to the
1616: `BV` internal data structures.
1618: Level: advanced
1620: .seealso: [](sec:bv), `BVGetMat()`, `BVRestoreArray()`
1621: @*/
1622: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1623: {
1624: PetscFunctionBegin;
1626: BVCheckSizes(bv,1);
1627: PetscAssertPointer(A,2);
1628: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1629: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1630: PetscUseTypeMethod(bv,restoremat,A);
1631: PetscFunctionReturn(PETSC_SUCCESS);
1632: }
1634: /*
1635: Copy all user-provided attributes of V to another BV object W
1636: */
1637: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1638: {
1639: PetscFunctionBegin;
1640: PetscCall(PetscLayoutReference(V->map,&W->map));
1641: PetscCall(BVSetVecType(W,V->vtype));
1642: W->ld = V->ld;
1643: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1644: W->orthog_type = V->orthog_type;
1645: W->orthog_ref = V->orthog_ref;
1646: W->orthog_eta = V->orthog_eta;
1647: W->orthog_block = V->orthog_block;
1648: PetscCall(PetscObjectReference((PetscObject)V->matrix));
1649: W->matrix = V->matrix;
1650: W->indef = V->indef;
1651: W->vmm = V->vmm;
1652: W->rrandom = V->rrandom;
1653: W->deftol = V->deftol;
1654: PetscCall(PetscObjectReference((PetscObject)V->rand));
1655: W->rand = V->rand;
1656: W->sfocalled = V->sfocalled;
1657: PetscTryTypeMethod(V,duplicate,W);
1658: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1659: PetscFunctionReturn(PETSC_SUCCESS);
1660: }
1662: /*@
1663: BVDuplicate - Creates a new basis vectors object of the same type and
1664: dimensions as an existing one.
1666: Collective
1668: Input Parameter:
1669: . V - basis vectors context
1671: Output Parameter:
1672: . W - location to put the new `BV`
1674: Notes:
1675: The new `BV` has the same type and dimensions as `V`. Also, the inner
1676: product matrix and orthogonalization options are copied.
1678: `BVDuplicate()` DOES NOT COPY the entries, but rather allocates storage
1679: for the new basis vectors. Use `BVCopy()` to copy the content.
1681: Level: intermediate
1683: .seealso: [](sec:bv), `BVDuplicateResize()`, `BVCreate()`, `BVSetSizesFromVec()`, `BVCopy()`
1684: @*/
1685: PetscErrorCode BVDuplicate(BV V,BV *W)
1686: {
1687: PetscFunctionBegin;
1690: BVCheckSizes(V,1);
1691: PetscAssertPointer(W,2);
1692: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1693: (*W)->N = V->N;
1694: (*W)->n = V->n;
1695: (*W)->m = V->m;
1696: (*W)->k = V->m;
1697: PetscCall(BVDuplicate_Private(V,*W));
1698: PetscFunctionReturn(PETSC_SUCCESS);
1699: }
1701: /*@
1702: BVDuplicateResize - Creates a new basis vectors object of the same type and
1703: dimensions as an existing one, but with possibly different number of columns.
1705: Collective
1707: Input Parameters:
1708: + V - basis vectors context
1709: - m - the new number of columns
1711: Output Parameter:
1712: . W - location to put the new BV
1714: Note:
1715: This is equivalent to a call to `BVDuplicate()` followed by `BVResize()`. The
1716: contents of `V` are not copied to `W`.
1718: Level: intermediate
1720: .seealso: [](sec:bv), `BVDuplicate()`, `BVResize()`
1721: @*/
1722: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1723: {
1724: PetscFunctionBegin;
1727: BVCheckSizes(V,1);
1729: PetscAssertPointer(W,3);
1730: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1731: (*W)->N = V->N;
1732: (*W)->n = V->n;
1733: (*W)->m = m;
1734: (*W)->k = m;
1735: PetscCall(BVDuplicate_Private(V,*W));
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: /*@
1740: BVGetCachedBV - Returns a `BV` object stored internally that holds the
1741: result of $BX$ after a call to `BVApplyMatrixBV()`.
1743: Collective
1745: Input Parameter:
1746: . bv - the basis vectors context
1748: Output Parameter:
1749: . cached - the cached `BV`
1751: Note:
1752: The cached `BV` is created if not available previously.
1754: Level: developer
1756: .seealso: [](sec:bv), `BVApplyMatrixBV()`
1757: @*/
1758: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1759: {
1760: PetscFunctionBegin;
1762: PetscAssertPointer(cached,2);
1763: BVCheckSizes(bv,1);
1764: if (!bv->cached) {
1765: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1766: bv->cached->N = bv->N;
1767: bv->cached->n = bv->n;
1768: bv->cached->m = bv->m;
1769: bv->cached->k = bv->m;
1770: PetscCall(BVDuplicate_Private(bv,bv->cached));
1771: }
1772: *cached = bv->cached;
1773: PetscFunctionReturn(PETSC_SUCCESS);
1774: }
1776: /*@
1777: BVCopy - Copies a basis vector object into another one, $W \leftarrow V$.
1779: Logically Collective
1781: Input Parameter:
1782: . V - basis vectors context
1784: Output Parameter:
1785: . W - the copy
1787: Note:
1788: Both `V` and `W` must be distributed in the same manner; local copies are
1789: done. Only active columns (excluding the leading ones) are copied.
1790: In the destination `W`, columns are overwritten starting from the leading ones.
1791: Constraints are not copied.
1793: Level: beginner
1795: .seealso: [](sec:bv), `BVCopyVec()`, `BVCopyColumn()`, `BVDuplicate()`, `BVSetActiveColumns()`
1796: @*/
1797: PetscErrorCode BVCopy(BV V,BV W)
1798: {
1799: PetscScalar *womega;
1800: const PetscScalar *vomega;
1802: PetscFunctionBegin;
1805: BVCheckSizes(V,1);
1806: BVCheckOp(V,1,copy);
1809: BVCheckSizes(W,2);
1810: PetscCheckSameTypeAndComm(V,1,W,2);
1811: 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);
1812: 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);
1813: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1815: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1816: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1817: /* copy signature */
1818: PetscCall(BV_AllocateSignature(W));
1819: PetscCall(VecGetArrayRead(V->omega,&vomega));
1820: PetscCall(VecGetArray(W->omega,&womega));
1821: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1822: PetscCall(VecRestoreArray(W->omega,&womega));
1823: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1824: }
1825: PetscUseTypeMethod(V,copy,W);
1826: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1827: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1828: PetscFunctionReturn(PETSC_SUCCESS);
1829: }
1831: /*@
1832: BVCopyVec - Copies one of the columns of a basis vectors object into a `Vec`.
1834: Logically Collective
1836: Input Parameters:
1837: + V - basis vectors context
1838: - j - the index of the column to be copied
1840: Output Parameter:
1841: . w - the copied column
1843: Note:
1844: Both `V` and `w` must be distributed in the same manner; local copies are done.
1846: Level: beginner
1848: .seealso: [](sec:bv), `BVCopy()`, `BVCopyColumn()`, `BVInsertVec()`
1849: @*/
1850: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1851: {
1852: PetscInt n,N;
1853: Vec z;
1855: PetscFunctionBegin;
1858: BVCheckSizes(V,1);
1861: PetscCheckSameComm(V,1,w,3);
1863: PetscCall(VecGetSize(w,&N));
1864: PetscCall(VecGetLocalSize(w,&n));
1865: 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);
1867: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1868: PetscCall(BVGetColumn(V,j,&z));
1869: PetscCall(VecCopy(z,w));
1870: PetscCall(BVRestoreColumn(V,j,&z));
1871: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1872: PetscFunctionReturn(PETSC_SUCCESS);
1873: }
1875: /*@
1876: BVCopyColumn - Copies the values from one of the columns to another one.
1878: Logically Collective
1880: Input Parameters:
1881: + V - basis vectors context
1882: . j - the index of the source column
1883: - i - the index of the destination column
1885: Level: beginner
1887: .seealso: [](sec:bv), `BVCopy()`, `BVCopyVec()`
1888: @*/
1889: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1890: {
1891: PetscScalar *omega;
1893: PetscFunctionBegin;
1896: BVCheckSizes(V,1);
1899: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1901: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1902: if (V->omega) {
1903: PetscCall(VecGetArray(V->omega,&omega));
1904: omega[i] = omega[j];
1905: PetscCall(VecRestoreArray(V->omega,&omega));
1906: }
1907: PetscUseTypeMethod(V,copycolumn,j,i);
1908: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1909: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1910: PetscFunctionReturn(PETSC_SUCCESS);
1911: }
1913: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1914: {
1915: PetscInt ncols;
1917: PetscFunctionBegin;
1918: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1919: if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1920: if (!*split) {
1921: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1922: (*split)->issplit = left? 1: 2;
1923: (*split)->splitparent = bv;
1924: (*split)->N = bv->N;
1925: (*split)->n = bv->n;
1926: (*split)->m = bv->m;
1927: (*split)->k = bv->m;
1928: PetscCall(BVDuplicate_Private(bv,*split));
1929: }
1930: (*split)->l = 0;
1931: (*split)->k = left? bv->l: bv->k-bv->l;
1932: (*split)->nc = left? bv->nc: 0;
1933: (*split)->m = ncols-(*split)->nc;
1934: if ((*split)->nc) {
1935: (*split)->ci[0] = -(*split)->nc-1;
1936: (*split)->ci[1] = -(*split)->nc-1;
1937: }
1938: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1939: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1940: PetscFunctionReturn(PETSC_SUCCESS);
1941: }
1943: /*@
1944: BVGetSplit - Splits the `BV` object into two `BV` objects that share the
1945: internal data, one of them containing the leading columns and the other
1946: one containing the remaining columns.
1948: Collective
1950: Input Parameter:
1951: . bv - the basis vectors context
1953: Output Parameters:
1954: + L - left `BV` containing leading columns (can be `NULL`)
1955: - R - right `BV` containing remaining columns (can be `NULL`)
1957: Notes:
1958: The columns are split in two sets. The leading columns (including the
1959: constraints) are assigned to the left `BV` and the remaining columns
1960: are assigned to the right `BV`. The number of leading columns, as
1961: specified with `BVSetActiveColumns()`, must be between `1` and `m-1` (to
1962: guarantee that both `L` and `R` have at least one column).
1964: The returned `BV`s must be seen as references (not copies) of the input
1965: `BV`, that is, modifying them will change the entries of `bv` as well.
1966: The returned `BV`s must not be destroyed. `BVRestoreSplit()` must be called
1967: when they are no longer needed.
1969: Pass `NULL` for any of the output `BV`s that is not needed.
1971: Level: advanced
1973: .seealso: [](sec:bv), `BVRestoreSplit()`, `BVSetActiveColumns()`, `BVSetNumConstraints()`, `BVGetSplitRows()`
1974: @*/
1975: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1976: {
1977: PetscFunctionBegin;
1980: BVCheckSizes(bv,1);
1981: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1982: PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1983: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1984: bv->lsplit = bv->nc+bv->l;
1985: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1986: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1987: if (L) *L = bv->L;
1988: if (R) *R = bv->R;
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: /*@
1993: BVRestoreSplit - Restore the BV objects obtained with `BVGetSplit()`.
1995: Logically Collective
1997: Input Parameters:
1998: + bv - the basis vectors context
1999: . L - left BV obtained with BVGetSplit()
2000: - R - right BV obtained with BVGetSplit()
2002: Note:
2003: The arguments must match the corresponding call to `BVGetSplit()`.
2005: Level: advanced
2007: .seealso: [](sec:bv), `BVGetSplit()`
2008: @*/
2009: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
2010: {
2011: PetscFunctionBegin;
2014: BVCheckSizes(bv,1);
2015: PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
2016: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
2017: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
2018: 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()");
2019: 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()");
2021: PetscTryTypeMethod(bv,restoresplit,L,R);
2022: bv->lsplit = 0;
2023: if (L) *L = NULL;
2024: if (R) *R = NULL;
2025: PetscFunctionReturn(PETSC_SUCCESS);
2026: }
2028: /*
2029: Copy all user-provided attributes of V to another BV object W with different layout
2030: */
2031: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
2032: {
2033: PetscFunctionBegin;
2034: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2035: PetscCall(BVSetVecType(W,V->vtype));
2036: W->ld = V->ld;
2037: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2038: W->orthog_type = V->orthog_type;
2039: W->orthog_ref = V->orthog_ref;
2040: W->orthog_eta = V->orthog_eta;
2041: W->orthog_block = V->orthog_block;
2042: W->vmm = V->vmm;
2043: W->rrandom = V->rrandom;
2044: W->deftol = V->deftol;
2045: PetscCall(PetscObjectReference((PetscObject)V->rand));
2046: W->rand = V->rand;
2047: W->sfocalled = V->sfocalled;
2048: PetscTryTypeMethod(V,duplicate,W);
2049: PetscCall(PetscObjectStateIncrease((PetscObject)W));
2050: PetscFunctionReturn(PETSC_SUCCESS);
2051: }
2053: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2054: {
2055: PetscInt rstart,rend,lstart;
2056: PetscInt N,n;
2057: PetscBool contig;
2059: PetscFunctionBegin;
2060: PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2061: PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2062: PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2063: if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2064: else bv->lsplit = -lstart;
2065: PetscCall(ISGetSize(is,&N));
2066: PetscCall(ISGetLocalSize(is,&n));
2067: if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2068: if (!*split) {
2069: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2070: (*split)->issplit = top? -1: -2;
2071: (*split)->splitparent = bv;
2072: (*split)->N = N;
2073: (*split)->n = n;
2074: (*split)->m = bv->m;
2075: PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2076: }
2077: (*split)->k = bv->k;
2078: (*split)->l = bv->l;
2079: (*split)->nc = bv->nc;
2080: if ((*split)->nc) {
2081: (*split)->ci[0] = -(*split)->nc-1;
2082: (*split)->ci[1] = -(*split)->nc-1;
2083: }
2084: if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2085: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2086: PetscFunctionReturn(PETSC_SUCCESS);
2087: }
2089: /*@
2090: BVGetSplitRows - Splits the `BV` object into two `BV` objects that share the
2091: internal data, using a disjoint horizontal splitting.
2093: Collective
2095: Input Parameters:
2096: + bv - the basis vectors context
2097: . isup - the index set that defines the upper part of the horizontal splitting
2098: - islo - the index set that defines the lower part of the horizontal splitting
2100: Output Parameters:
2101: + U - the resulting `BV` containing the upper rows
2102: - L - the resulting `BV` containing the lower rows
2104: Notes:
2105: The index sets must be such that every MPI process can extract the selected
2106: rows from its local part of the input `BV`, and this part must be contiguous.
2107: With one process, `isup` will list contiguous indices starting from 0, and `islo`
2108: will contain the remaining indices, hence we refer to upper and lower part.
2109: However, with several processes the indices will be interleaved because
2110: `isup` will refer to the upper part of the local array.
2112: The intended use of this function is with matrices of `MATNEST` type, where
2113: `MatNestGetISs()` will return the appropriate index sets.
2115: The returned `BV`s must be seen as references (not copies) of the input
2116: `BV`, that is, modifying them will change the entries of `bv` as well.
2117: The returned `BV`s must not be destroyed. `BVRestoreSplitRows()` must be called
2118: when they are no longer needed.
2120: Pass `NULL` for any of the output `BV`s that is not needed.
2122: Level: advanced
2124: .seealso: [](sec:bv), `BVRestoreSplitRows()`, `BVGetSplit()`
2125: @*/
2126: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2127: {
2128: PetscFunctionBegin;
2131: BVCheckSizes(bv,1);
2134: PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2135: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2136: if (U) {
2137: PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2138: *U = bv->R;
2139: }
2140: if (L) {
2141: PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2142: *L = bv->L;
2143: }
2144: PetscFunctionReturn(PETSC_SUCCESS);
2145: }
2147: /*@
2148: BVRestoreSplitRows - Restore the BV objects obtained with `BVGetSplitRows()`.
2150: Logically Collective
2152: Input Parameters:
2153: + bv - the basis vectors context
2154: . isup - the index set that defines the upper part of the horizontal splitting
2155: . islo - the index set that defines the lower part of the horizontal splitting
2156: . U - the `BV` containing the upper rows
2157: - L - the `BV` containing the lower rows
2159: Note:
2160: The arguments must match the corresponding call to `BVGetSplitRows()`.
2162: Level: advanced
2164: .seealso: [](sec:bv), `BVGetSplitRows()`
2165: @*/
2166: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2167: {
2168: PetscFunctionBegin;
2171: BVCheckSizes(bv,1);
2172: PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2173: 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()");
2174: 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()");
2175: PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2176: bv->lsplit = 0;
2177: PetscFunctionReturn(PETSC_SUCCESS);
2178: }
2180: /*@
2181: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2182: definite inner product.
2184: Logically Collective
2186: Input Parameters:
2187: + bv - basis vectors
2188: - deftol - the tolerance
2190: Options Database Key:
2191: . -bv_definite_tol deftol - the tolerance
2193: Notes:
2194: When using a non-standard inner product, see `BVSetMatrix()`, the solver needs
2195: to compute $\sqrt{z^*B z}$ for various vectors $z$. If the inner product has not
2196: been declared indefinite, the value $z^*B z$ must be positive, but due to
2197: rounding error a tiny value may become negative. A tolerance is used to
2198: detect this situation. Likewise, in complex arithmetic $z^*B z$ should be
2199: real, and we use the same tolerance to check whether a nonzero imaginary part
2200: can be considered negligible.
2202: This function sets this tolerance, which defaults to `10*PETSC_MACHINE_EPSILON`.
2203: The default value should be good for most applications.
2205: Level: advanced
2207: .seealso: [](sec:bv), `BVSetMatrix()`
2208: @*/
2209: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2210: {
2211: PetscFunctionBegin;
2214: if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2215: else {
2216: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2217: bv->deftol = deftol;
2218: }
2219: PetscFunctionReturn(PETSC_SUCCESS);
2220: }
2222: /*@
2223: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2224: inner product.
2226: Not Collective
2228: Input Parameter:
2229: . bv - the basis vectors
2231: Output Parameter:
2232: . deftol - the tolerance
2234: Level: advanced
2236: .seealso: [](sec:bv), `BVSetDefiniteTolerance()`
2237: @*/
2238: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2239: {
2240: PetscFunctionBegin;
2242: PetscAssertPointer(deftol,2);
2243: *deftol = bv->deftol;
2244: PetscFunctionReturn(PETSC_SUCCESS);
2245: }
2247: /*@
2248: BVSetLeadingDimension - Set the leading dimension to be used for storing the `BV` data.
2250: Not Collective
2252: Input Parameters:
2253: + bv - basis vectors
2254: - ld - the leading dimension
2256: Notes:
2257: This parameter is relevant for `BVMAT`, though it might be employed in other types
2258: as well.
2260: When the internal data of the `BV` is stored as a dense matrix, the leading dimension
2261: has the same meaning as in `MatDenseSetLDA()`, i.e., the distance in number of
2262: elements from one entry of the matrix to the one in the next column at the same
2263: row. The leading dimension refers to the local array, and hence can be different
2264: in different processes.
2266: The user does not need to change this parameter. The default value is equal to the
2267: number of local rows, but this value may be increased a little to guarantee alignment
2268: (especially in the case of GPU storage).
2270: Level: advanced
2272: .seealso: [](sec:bv), `BVGetLeadingDimension()`
2273: @*/
2274: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2275: {
2276: PetscFunctionBegin;
2279: 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");
2280: if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
2281: else {
2282: PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2283: bv->ld = ld;
2284: }
2285: PetscFunctionReturn(PETSC_SUCCESS);
2286: }
2288: /*@
2289: BVGetLeadingDimension - Returns the leading dimension of the `BV`.
2291: Not Collective
2293: Input Parameter:
2294: . bv - the basis vectors
2296: Output Parameter:
2297: . ld - the leading dimension
2299: Level: advanced
2301: Notes:
2302: The returned value may be different in different processes.
2304: The leading dimension must be used when accessing the internal array via
2305: `BVGetArray()` or `BVGetArrayRead()`.
2307: .seealso: [](sec:bv), `BVSetLeadingDimension()`, `BVGetArray()`, `BVGetArrayRead()`
2308: @*/
2309: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2310: {
2311: PetscFunctionBegin;
2313: PetscAssertPointer(ld,2);
2314: *ld = bv->ld;
2315: PetscFunctionReturn(PETSC_SUCCESS);
2316: }