Actual source code: bvbasic.c
slepc-3.22.1 2024-10-28
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: 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: 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 processor calls this with N of PETSC_DECIDE then all processors must,
104: otherwise the program will hang.
106: Level: beginner
108: .seealso: BVSetSizesFromVec(), BVGetSizes(), BVResize()
109: @*/
110: PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)
111: {
112: PetscInt ma;
113: PetscMPIInt size;
115: PetscFunctionBegin;
119: PetscCheck(N<0 || n<=N,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT,n,N);
120: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
121: PetscCheck((bv->n<0 && bv->N<0) || (bv->n==n && bv->N==N),PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change/reset vector sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global",n,N,bv->n,bv->N);
122: PetscCheck(bv->m<=0 || bv->m==m,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change the number of columns to %" PetscInt_FMT " after previously setting it to %" PetscInt_FMT "; use BVResize()",m,bv->m);
123: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
124: bv->n = n;
125: bv->N = N;
126: bv->m = m;
127: bv->k = m;
128: /* create layout and get actual dimensions */
129: PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)bv),&bv->map));
130: PetscCall(PetscLayoutSetSize(bv->map,bv->N));
131: PetscCall(PetscLayoutSetLocalSize(bv->map,bv->n));
132: PetscCall(PetscLayoutSetUp(bv->map));
133: PetscCall(PetscLayoutGetSize(bv->map,&bv->N));
134: PetscCall(PetscLayoutGetLocalSize(bv->map,&bv->n));
135: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
136: PetscCall(BVSetVecType(bv,(size==1)?VECSEQ:VECMPI));
137: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
138: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
139: PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
140: }
141: if (bv->ops->create) {
142: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
143: PetscUseTypeMethod(bv,create);
144: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
145: bv->ops->create = NULL;
146: bv->defersfo = PETSC_FALSE;
147: }
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
153: Local and global sizes are specified indirectly by passing a template vector.
155: Collective
157: Input Parameters:
158: + bv - the basis vectors
159: . t - the template vector
160: - m - the number of columns
162: Level: beginner
164: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
165: @*/
166: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
167: {
168: PetscInt ma;
169: PetscLayout map;
170: VecType vtype;
172: PetscFunctionBegin;
175: PetscCheckSameComm(bv,1,t,2);
177: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
178: PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
179: PetscCall(VecGetType(t,&vtype));
180: PetscCall(BVSetVecType(bv,vtype));
181: PetscCall(VecGetLayout(t,&map));
182: PetscCall(PetscLayoutReference(map,&bv->map));
183: PetscCall(VecGetSize(t,&bv->N));
184: PetscCall(VecGetLocalSize(t,&bv->n));
185: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
186: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
187: PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
188: }
189: bv->m = m;
190: bv->k = m;
191: if (bv->ops->create) {
192: PetscUseTypeMethod(bv,create);
193: bv->ops->create = NULL;
194: bv->defersfo = PETSC_FALSE;
195: }
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: BVGetSizes - Returns the local and global sizes, and the number of columns.
202: Not Collective
204: Input Parameter:
205: . bv - the basis vectors
207: Output Parameters:
208: + n - the local size
209: . N - the global size
210: - m - the number of columns
212: Note:
213: Normal usage requires that bv has already been given its sizes, otherwise
214: the call fails. However, this function can also be used to determine if
215: a BV object has been initialized completely (sizes and type). For this,
216: call with n=NULL and N=NULL, then a return value of m=0 indicates that
217: the BV object is not ready for use yet.
219: Level: beginner
221: .seealso: BVSetSizes(), BVSetSizesFromVec()
222: @*/
223: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
224: {
225: PetscFunctionBegin;
226: if (!bv) {
227: if (m && !n && !N) *m = 0;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: if (n || N) BVCheckSizes(bv,1);
232: if (n) *n = bv->n;
233: if (N) *N = bv->N;
234: if (m) *m = bv->m;
235: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: BVSetNumConstraints - Set the number of constraints.
242: Logically Collective
244: Input Parameters:
245: + V - basis vectors
246: - nc - number of constraints
248: Notes:
249: This function sets the number of constraints to nc and marks all remaining
250: columns as regular. Normal user would call BVInsertConstraints() instead.
252: If nc is smaller than the previously set value, then some of the constraints
253: are discarded. In particular, using nc=0 removes all constraints preserving
254: the content of regular columns.
256: Level: developer
258: .seealso: BVInsertConstraints()
259: @*/
260: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
261: {
262: PetscInt total,diff,i;
263: Vec x,y;
265: PetscFunctionBegin;
268: PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
270: BVCheckSizes(V,1);
271: PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
273: diff = nc-V->nc;
274: if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
275: total = V->nc+V->m;
276: PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
277: if (diff<0) { /* lessen constraints, shift contents of BV */
278: for (i=0;i<V->m;i++) {
279: PetscCall(BVGetColumn(V,i,&x));
280: PetscCall(BVGetColumn(V,i+diff,&y));
281: PetscCall(VecCopy(x,y));
282: PetscCall(BVRestoreColumn(V,i,&x));
283: PetscCall(BVRestoreColumn(V,i+diff,&y));
284: }
285: }
286: V->nc = nc;
287: V->ci[0] = -V->nc-1;
288: V->ci[1] = -V->nc-1;
289: V->m = total-nc;
290: V->l = PetscMin(V->l,V->m);
291: V->k = PetscMin(V->k,V->m);
292: PetscCall(PetscObjectStateIncrease((PetscObject)V));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: BVGetNumConstraints - Returns the number of constraints.
299: Not Collective
301: Input Parameter:
302: . bv - the basis vectors
304: Output Parameters:
305: . nc - the number of constraints
307: Level: advanced
309: .seealso: BVGetSizes(), BVInsertConstraints()
310: @*/
311: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
312: {
313: PetscFunctionBegin;
315: PetscAssertPointer(nc,2);
316: *nc = bv->nc;
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@
321: BVResize - Change the number of columns.
323: Collective
325: Input Parameters:
326: + bv - the basis vectors
327: . m - the new number of columns
328: - copy - a flag indicating whether current values should be kept
330: Note:
331: Internal storage is reallocated. If the copy flag is set to true, then
332: the contents are copied to the leading part of the new space.
334: Level: advanced
336: .seealso: BVSetSizes(), BVSetSizesFromVec()
337: @*/
338: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
339: {
340: PetscScalar *array;
341: const PetscScalar *omega;
342: Vec v;
344: PetscFunctionBegin;
349: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
350: PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
351: if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
352: BVCheckOp(bv,1,resize);
354: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
355: PetscUseTypeMethod(bv,resize,m,copy);
356: PetscCall(VecDestroy(&bv->buffer));
357: PetscCall(BVDestroy(&bv->cached));
358: PetscCall(PetscFree2(bv->h,bv->c));
359: if (bv->omega) {
360: if (bv->cuda) {
361: #if defined(PETSC_HAVE_CUDA)
362: PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
363: #else
364: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
365: #endif
366: } else if (bv->hip) {
367: #if defined(PETSC_HAVE_HIP)
368: PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
369: #else
370: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
371: #endif
372: } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
373: if (copy) {
374: PetscCall(VecGetArray(v,&array));
375: PetscCall(VecGetArrayRead(bv->omega,&omega));
376: PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
377: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
378: PetscCall(VecRestoreArray(v,&array));
379: } else PetscCall(VecSet(v,1.0));
380: PetscCall(VecDestroy(&bv->omega));
381: bv->omega = v;
382: }
383: bv->m = m;
384: bv->k = PetscMin(bv->k,m);
385: bv->l = PetscMin(bv->l,m);
386: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
387: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: BVSetActiveColumns - Specify the columns that will be involved in operations.
394: Logically Collective
396: Input Parameters:
397: + bv - the basis vectors context
398: . l - number of leading columns
399: - k - number of active columns
401: Notes:
402: In operations such as BVMult() or BVDot(), only the first k columns are
403: considered. This is useful when the BV is filled from left to right, so
404: the last m-k columns do not have relevant information.
406: Also in operations such as BVMult() or BVDot(), the first l columns are
407: normally not included in the computation. See the manpage of each
408: operation.
410: In orthogonalization operations, the first l columns are treated
411: differently, they participate in the orthogonalization but the computed
412: coefficients are not stored.
414: Use PETSC_CURRENT to leave any of the values unchanged. Use PETSC_DETERMINE
415: to set l to the minimum value (0) and k to the maximum (m).
417: Level: intermediate
419: .seealso: BVGetActiveColumns(), BVSetSizes()
420: @*/
421: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
422: {
423: PetscFunctionBegin;
427: BVCheckSizes(bv,1);
428: if (PetscUnlikely(k == PETSC_DETERMINE)) {
429: bv->k = bv->m;
430: } else if (k != PETSC_CURRENT) {
431: PetscCheck(k>=0 && k<=bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k (%" PetscInt_FMT "). Must be between 0 and m (%" PetscInt_FMT ")",k,bv->m);
432: bv->k = k;
433: }
434: if (PetscUnlikely(l == PETSC_DETERMINE)) {
435: bv->l = 0;
436: } else if (l != PETSC_CURRENT) {
437: PetscCheck(l>=0 && l<=bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l (%" PetscInt_FMT "). Must be between 0 and k (%" PetscInt_FMT ")",l,bv->k);
438: bv->l = l;
439: }
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@
444: BVGetActiveColumns - Returns the current active dimensions.
446: Not Collective
448: Input Parameter:
449: . bv - the basis vectors context
451: Output Parameters:
452: + l - number of leading columns
453: - k - number of active columns
455: Level: intermediate
457: .seealso: BVSetActiveColumns()
458: @*/
459: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
460: {
461: PetscFunctionBegin;
463: if (l) *l = bv->l;
464: if (k) *k = bv->k;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
471: Collective
473: Input Parameters:
474: + bv - the basis vectors context
475: . B - a symmetric matrix (may be NULL)
476: - indef - a flag indicating if the matrix is indefinite
478: Notes:
479: This is used to specify a non-standard inner product, whose matrix
480: representation is given by B. Then, all inner products required during
481: orthogonalization are computed as (x,y)_B=y^H*B*x rather than the
482: standard form (x,y)=y^H*x.
484: Matrix B must be real symmetric (or complex Hermitian). A genuine inner
485: product requires that B is also positive (semi-)definite. However, we
486: also allow for an indefinite B (setting indef=PETSC_TRUE), in which
487: case the orthogonalization uses an indefinite inner product.
489: This affects operations BVDot(), BVNorm(), BVOrthogonalize(), and variants.
491: Setting B=NULL has the same effect as if the identity matrix was passed.
493: Level: advanced
495: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
496: @*/
497: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
498: {
499: PetscInt m,n;
501: PetscFunctionBegin;
504: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
505: if (B) {
507: PetscCall(MatGetLocalSize(B,&m,&n));
508: PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
509: PetscCheck(!bv->m || bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension BV %" PetscInt_FMT ", Mat %" PetscInt_FMT,bv->n,n);
510: }
511: if (B) PetscCall(PetscObjectReference((PetscObject)B));
512: PetscCall(MatDestroy(&bv->matrix));
513: bv->matrix = B;
514: bv->indef = indef;
515: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
516: if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
517: if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
518: }
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: BVGetMatrix - Retrieves the matrix representation of the inner product.
525: Not Collective
527: Input Parameter:
528: . bv - the basis vectors context
530: Output Parameters:
531: + B - the matrix of the inner product (may be NULL)
532: - indef - the flag indicating if the matrix is indefinite
534: Level: advanced
536: .seealso: BVSetMatrix()
537: @*/
538: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
539: {
540: PetscFunctionBegin;
542: if (B) *B = bv->matrix;
543: if (indef) *indef = bv->indef;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: BVApplyMatrix - Multiplies a vector by the matrix representation of the
549: inner product.
551: Neighbor-wise Collective
553: Input Parameters:
554: + bv - the basis vectors context
555: - x - the vector
557: Output Parameter:
558: . y - the result
560: Note:
561: If no matrix was specified this function copies the vector.
563: Level: advanced
565: .seealso: BVSetMatrix(), BVApplyMatrixBV()
566: @*/
567: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
568: {
569: PetscFunctionBegin;
573: if (bv->matrix) {
574: PetscCall(BV_IPMatMult(bv,x));
575: PetscCall(VecCopy(bv->Bx,y));
576: } else PetscCall(VecCopy(x,y));
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
582: of the inner product.
584: Neighbor-wise Collective
586: Input Parameter:
587: . X - the basis vectors context
589: Output Parameter:
590: . Y - the basis vectors to store the result (optional)
592: Note:
593: This function computes Y = B*X, where B is the matrix given with
594: BVSetMatrix(). This operation is computed as in BVMatMult().
595: If no matrix was specified, then it just copies Y = X.
597: If no Y is given, the result is stored internally in the cached BV.
599: Level: developer
601: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
602: @*/
603: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
604: {
605: PetscFunctionBegin;
607: if (Y) {
609: if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
610: else PetscCall(BVCopy(X,Y));
611: } else PetscCall(BV_IPMatMultBV(X));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@
616: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
618: Logically Collective
620: Input Parameters:
621: + bv - the basis vectors context
622: - omega - a vector representing the diagonal of the signature matrix
624: Note:
625: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
627: Level: developer
629: .seealso: BVSetMatrix(), BVGetSignature()
630: @*/
631: PetscErrorCode BVSetSignature(BV bv,Vec omega)
632: {
633: PetscInt i,n;
634: const PetscScalar *pomega;
635: PetscScalar *intern;
637: PetscFunctionBegin;
639: BVCheckSizes(bv,1);
643: PetscCall(VecGetSize(omega,&n));
644: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
645: PetscCall(BV_AllocateSignature(bv));
646: if (bv->indef) {
647: PetscCall(VecGetArrayRead(omega,&pomega));
648: PetscCall(VecGetArray(bv->omega,&intern));
649: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
650: PetscCall(VecRestoreArray(bv->omega,&intern));
651: PetscCall(VecRestoreArrayRead(omega,&pomega));
652: } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
653: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: /*@
658: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
660: Not Collective
662: Input Parameter:
663: . bv - the basis vectors context
665: Output Parameter:
666: . omega - a vector representing the diagonal of the signature matrix
668: Note:
669: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
671: Level: developer
673: .seealso: BVSetMatrix(), BVSetSignature()
674: @*/
675: PetscErrorCode BVGetSignature(BV bv,Vec omega)
676: {
677: PetscInt i,n;
678: PetscScalar *pomega;
679: const PetscScalar *intern;
681: PetscFunctionBegin;
683: BVCheckSizes(bv,1);
687: PetscCall(VecGetSize(omega,&n));
688: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
689: if (bv->indef && bv->omega) {
690: PetscCall(VecGetArray(omega,&pomega));
691: PetscCall(VecGetArrayRead(bv->omega,&intern));
692: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
693: PetscCall(VecRestoreArrayRead(bv->omega,&intern));
694: PetscCall(VecRestoreArray(omega,&pomega));
695: } else PetscCall(VecSet(omega,1.0));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: BVSetBufferVec - Attach a vector object to be used as buffer space for
701: several operations.
703: Collective
705: Input Parameters:
706: + bv - the basis vectors context)
707: - buffer - the vector
709: Notes:
710: Use BVGetBufferVec() to retrieve the vector (for example, to free it
711: at the end of the computations).
713: The vector must be sequential of length (nc+m)*m, where m is the number
714: of columns of bv and nc is the number of constraints.
716: Level: developer
718: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
719: @*/
720: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
721: {
722: PetscInt ld,n;
723: PetscMPIInt size;
725: PetscFunctionBegin;
728: BVCheckSizes(bv,1);
729: PetscCall(VecGetSize(buffer,&n));
730: ld = bv->m+bv->nc;
731: PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
732: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
733: PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
735: PetscCall(PetscObjectReference((PetscObject)buffer));
736: PetscCall(VecDestroy(&bv->buffer));
737: bv->buffer = buffer;
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: /*@
742: BVGetBufferVec - Obtain the buffer vector associated with the BV object.
744: Collective
746: Input Parameters:
747: . bv - the basis vectors context
749: Output Parameter:
750: . buffer - vector
752: Notes:
753: The vector is created if not available previously. It is a sequential vector
754: of length (nc+m)*m, where m is the number of columns of bv and nc is the number
755: of constraints.
757: Developer Notes:
758: The buffer vector is viewed as a column-major matrix with leading dimension
759: ld=nc+m, and m columns at most. In the most common usage, it has the structure
760: .vb
761: | | C |
762: |s|---|
763: | | H |
764: .ve
765: where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
766: related to orthogonalization against constraints (first nc rows), and s is the
767: first column that contains scratch values computed during Gram-Schmidt
768: orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
769: store the coefficients.
771: Level: developer
773: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
774: @*/
775: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
776: {
777: PetscInt ld;
779: PetscFunctionBegin;
781: PetscAssertPointer(buffer,2);
782: BVCheckSizes(bv,1);
783: if (!bv->buffer) {
784: ld = bv->m+bv->nc;
785: PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
786: PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
787: PetscCall(VecSetType(bv->buffer,bv->vtype));
788: }
789: *buffer = bv->buffer;
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: BVSetRandomContext - Sets the PetscRandom object associated with the BV,
795: to be used in operations that need random numbers.
797: Collective
799: Input Parameters:
800: + bv - the basis vectors context
801: - rand - the random number generator context
803: Level: advanced
805: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
806: @*/
807: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
808: {
809: PetscFunctionBegin;
812: PetscCheckSameComm(bv,1,rand,2);
813: PetscCall(PetscObjectReference((PetscObject)rand));
814: PetscCall(PetscRandomDestroy(&bv->rand));
815: bv->rand = rand;
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: /*@
820: BVGetRandomContext - Gets the PetscRandom object associated with the BV.
822: Collective
824: Input Parameter:
825: . bv - the basis vectors context
827: Output Parameter:
828: . rand - the random number generator context
830: Level: advanced
832: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
833: @*/
834: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
835: {
836: PetscFunctionBegin;
838: PetscAssertPointer(rand,2);
839: if (!bv->rand) {
840: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
841: if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
842: if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
843: if (bv->rrandom) {
844: PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
845: PetscCall(PetscRandomSeed(bv->rand));
846: }
847: }
848: *rand = bv->rand;
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: /*@
853: BVSetFromOptions - Sets BV options from the options database.
855: Collective
857: Input Parameter:
858: . bv - the basis vectors context
860: Level: beginner
862: .seealso: BVSetOptionsPrefix()
863: @*/
864: PetscErrorCode BVSetFromOptions(BV bv)
865: {
866: char type[256];
867: PetscBool flg1,flg2,flg3,flg4;
868: PetscReal r;
869: BVOrthogType otype;
870: BVOrthogRefineType orefine;
871: BVOrthogBlockType oblock;
873: PetscFunctionBegin;
875: PetscCall(BVRegisterAll());
876: PetscObjectOptionsBegin((PetscObject)bv);
877: PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
878: if (flg1) PetscCall(BVSetType(bv,type));
879: else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));
881: otype = bv->orthog_type;
882: PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
883: orefine = bv->orthog_ref;
884: PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
885: r = bv->orthog_eta;
886: PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
887: oblock = bv->orthog_block;
888: PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
889: if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
891: PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
893: PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
894: if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
896: /* undocumented option to generate random vectors that are independent of the number of processes */
897: PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
899: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
900: else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
901: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
902: PetscOptionsEnd();
903: bv->sfocalled = PETSC_TRUE;
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: /*@
908: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
909: vectors (classical or modified Gram-Schmidt with or without refinement), and
910: for the block-orthogonalization (simultaneous orthogonalization of a set of
911: vectors).
913: Logically Collective
915: Input Parameters:
916: + bv - the basis vectors context
917: . type - the method of vector orthogonalization
918: . refine - type of refinement
919: . eta - parameter for selective refinement
920: - block - the method of block orthogonalization
922: Options Database Keys:
923: + -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
924: (default) or mgs for Modified Gram-Schmidt orthogonalization
925: . -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
926: . -bv_orthog_eta <eta> - For setting the value of eta
927: - -bv_orthog_block <block> - Where <block> is the block-orthogonalization method
929: Notes:
930: The default settings work well for most problems.
932: The parameter eta should be a real value between 0 and 1, that is used only when
933: the refinement type is "ifneeded". Use PETSC_DETERMINE to set a reasonable
934: default value. Use PETSC_CURRENT to leave the current value unchanged.
936: When using several processors, MGS is likely to result in bad scalability.
938: If the method set for block orthogonalization is GS, then the computation
939: is done column by column with the vector orthogonalization.
941: Level: advanced
943: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
944: @*/
945: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
946: {
947: PetscFunctionBegin;
953: switch (type) {
954: case BV_ORTHOG_CGS:
955: case BV_ORTHOG_MGS:
956: bv->orthog_type = type;
957: break;
958: default:
959: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
960: }
961: switch (refine) {
962: case BV_ORTHOG_REFINE_NEVER:
963: case BV_ORTHOG_REFINE_IFNEEDED:
964: case BV_ORTHOG_REFINE_ALWAYS:
965: bv->orthog_ref = refine;
966: break;
967: default:
968: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
969: }
970: if (eta == (PetscReal)PETSC_DETERMINE) {
971: bv->orthog_eta = 0.7071;
972: } else if (eta != (PetscReal)PETSC_CURRENT) {
973: PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
974: bv->orthog_eta = eta;
975: }
976: switch (block) {
977: case BV_ORTHOG_BLOCK_GS:
978: case BV_ORTHOG_BLOCK_CHOL:
979: case BV_ORTHOG_BLOCK_TSQR:
980: case BV_ORTHOG_BLOCK_TSQRCHOL:
981: case BV_ORTHOG_BLOCK_SVQB:
982: bv->orthog_block = block;
983: break;
984: default:
985: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
986: }
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: /*@
991: BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.
993: Not Collective
995: Input Parameter:
996: . bv - basis vectors context
998: Output Parameters:
999: + type - the method of vector orthogonalization
1000: . refine - type of refinement
1001: . eta - parameter for selective refinement
1002: - block - the method of block orthogonalization
1004: Level: advanced
1006: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
1007: @*/
1008: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1009: {
1010: PetscFunctionBegin;
1012: if (type) *type = bv->orthog_type;
1013: if (refine) *refine = bv->orthog_ref;
1014: if (eta) *eta = bv->orthog_eta;
1015: if (block) *block = bv->orthog_block;
1016: PetscFunctionReturn(PETSC_SUCCESS);
1017: }
1019: /*@
1020: BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.
1022: Logically Collective
1024: Input Parameters:
1025: + bv - the basis vectors context
1026: - method - the method for the BVMatMult() operation
1028: Options Database Keys:
1029: . -bv_matmult <meth> - choose one of the methods: vecs, mat
1031: Notes:
1032: Allowed values are
1033: + BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1034: . BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1035: - BV_MATMULT_MAT_SAVE - this case is deprecated
1037: The default is BV_MATMULT_MAT except in the case of BVVECS.
1039: Level: advanced
1041: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1042: @*/
1043: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1044: {
1045: PetscFunctionBegin;
1048: switch (method) {
1049: case BV_MATMULT_VECS:
1050: case BV_MATMULT_MAT:
1051: bv->vmm = method;
1052: break;
1053: case BV_MATMULT_MAT_SAVE:
1054: PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1055: bv->vmm = BV_MATMULT_MAT;
1056: break;
1057: default:
1058: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1059: }
1060: PetscFunctionReturn(PETSC_SUCCESS);
1061: }
1063: /*@
1064: BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.
1066: Not Collective
1068: Input Parameter:
1069: . bv - basis vectors context
1071: Output Parameter:
1072: . method - the method for the BVMatMult() operation
1074: Level: advanced
1076: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1077: @*/
1078: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1079: {
1080: PetscFunctionBegin;
1082: PetscAssertPointer(method,2);
1083: *method = bv->vmm;
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: /*@
1088: BVGetColumn - Returns a Vec object that contains the entries of the
1089: requested column of the basis vectors object.
1091: Logically Collective
1093: Input Parameters:
1094: + bv - the basis vectors context
1095: - j - the index of the requested column
1097: Output Parameter:
1098: . v - vector containing the jth column
1100: Notes:
1101: The returned Vec must be seen as a reference (not a copy) of the BV
1102: column, that is, modifying the Vec will change the BV entries as well.
1104: The returned Vec must not be destroyed. BVRestoreColumn() must be
1105: called when it is no longer needed. At most, two columns can be fetched,
1106: that is, this function can only be called twice before the corresponding
1107: BVRestoreColumn() is invoked.
1109: A negative index j selects the i-th constraint, where i=-j. Constraints
1110: should not be modified.
1112: Level: beginner
1114: .seealso: BVRestoreColumn(), BVInsertConstraints()
1115: @*/
1116: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1117: {
1118: PetscInt l;
1120: PetscFunctionBegin;
1123: BVCheckSizes(bv,1);
1124: BVCheckOp(bv,1,getcolumn);
1126: 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);
1127: 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);
1128: 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);
1129: l = BVAvailableVec;
1130: PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1131: PetscUseTypeMethod(bv,getcolumn,j,v);
1132: bv->ci[l] = j;
1133: PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
1134: PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1135: *v = bv->cv[l];
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: /*@
1140: BVRestoreColumn - Restore a column obtained with BVGetColumn().
1142: Logically Collective
1144: Input Parameters:
1145: + bv - the basis vectors context
1146: . j - the index of the column
1147: - v - vector obtained with BVGetColumn()
1149: Note:
1150: The arguments must match the corresponding call to BVGetColumn().
1152: Level: beginner
1154: .seealso: BVGetColumn()
1155: @*/
1156: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1157: {
1158: PetscObjectId id;
1159: PetscObjectState st;
1160: PetscInt l;
1162: PetscFunctionBegin;
1165: BVCheckSizes(bv,1);
1167: PetscAssertPointer(v,3);
1169: 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);
1170: 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);
1171: 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);
1172: l = (j==bv->ci[0])? 0: 1;
1173: PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1174: PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1175: PetscCall(VecGetState(*v,&st));
1176: if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1177: PetscUseTypeMethod(bv,restorecolumn,j,v);
1178: bv->ci[l] = -bv->nc-1;
1179: bv->st[l] = -1;
1180: bv->id[l] = 0;
1181: *v = NULL;
1182: PetscFunctionReturn(PETSC_SUCCESS);
1183: }
1185: /*@C
1186: BVGetArray - Returns a pointer to a contiguous array that contains this
1187: processor's portion of the BV data.
1189: Logically Collective
1191: Input Parameters:
1192: . bv - the basis vectors context
1194: Output Parameter:
1195: . a - location to put pointer to the array
1197: Notes:
1198: BVRestoreArray() must be called when access to the array is no longer needed.
1199: This operation may imply a data copy, for BV types that do not store
1200: data contiguously in memory.
1202: The pointer will normally point to the first entry of the first column,
1203: but if the BV has constraints then these go before the regular columns.
1205: Note that for manipulating the pointer to the BV array, one must take into
1206: account the leading dimension, which might be different from the local
1207: number of rows, see BVGetLeadingDimension().
1209: Use BVGetArrayRead() for read-only access.
1211: Level: advanced
1213: .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArrayRead()
1214: @*/
1215: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1216: {
1217: PetscFunctionBegin;
1220: BVCheckSizes(bv,1);
1221: BVCheckOp(bv,1,getarray);
1222: PetscUseTypeMethod(bv,getarray,a);
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: /*@C
1227: BVRestoreArray - Restore the BV object after BVGetArray() has been called.
1229: Logically Collective
1231: Input Parameters:
1232: + bv - the basis vectors context
1233: - a - location of pointer to array obtained from BVGetArray()
1235: Note:
1236: This operation may imply a data copy, for BV types that do not store
1237: data contiguously in memory.
1239: Level: advanced
1241: .seealso: BVGetColumn()
1242: @*/
1243: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1244: {
1245: PetscFunctionBegin;
1248: BVCheckSizes(bv,1);
1249: PetscTryTypeMethod(bv,restorearray,a);
1250: if (a) *a = NULL;
1251: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@C
1256: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1257: contains this processor's portion of the BV data.
1259: Not Collective
1261: Input Parameters:
1262: . bv - the basis vectors context
1264: Output Parameter:
1265: . a - location to put pointer to the array
1267: Notes:
1268: BVRestoreArrayRead() must be called when access to the array is no
1269: longer needed. This operation may imply a data copy, for BV types that
1270: do not store data contiguously in memory.
1272: The pointer will normally point to the first entry of the first column,
1273: but if the BV has constraints then these go before the regular columns.
1275: Level: advanced
1277: .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArray()
1278: @*/
1279: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1280: {
1281: PetscFunctionBegin;
1284: BVCheckSizes(bv,1);
1285: BVCheckOp(bv,1,getarrayread);
1286: PetscUseTypeMethod(bv,getarrayread,a);
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: /*@C
1291: BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1292: been called.
1294: Not Collective
1296: Input Parameters:
1297: + bv - the basis vectors context
1298: - a - location of pointer to array obtained from BVGetArrayRead()
1300: Level: advanced
1302: .seealso: BVGetColumn()
1303: @*/
1304: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1305: {
1306: PetscFunctionBegin;
1309: BVCheckSizes(bv,1);
1310: PetscTryTypeMethod(bv,restorearrayread,a);
1311: if (a) *a = NULL;
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: /*@
1316: BVCreateVec - Creates a new Vec object with the same type and dimensions
1317: as the columns of the basis vectors object.
1319: Collective
1321: Input Parameter:
1322: . bv - the basis vectors context
1324: Output Parameter:
1325: . v - the new vector
1327: Note:
1328: The user is responsible of destroying the returned vector.
1330: Level: beginner
1332: .seealso: BVCreateMat(), BVCreateVecEmpty()
1333: @*/
1334: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1335: {
1336: PetscFunctionBegin;
1338: BVCheckSizes(bv,1);
1339: PetscAssertPointer(v,2);
1340: PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1341: PetscCall(VecSetLayout(*v,bv->map));
1342: PetscCall(VecSetType(*v,bv->vtype));
1343: PetscCall(VecSetUp(*v));
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: /*@
1348: BVCreateVecEmpty - Creates a new Vec object with the same type and dimensions
1349: as the columns of the basis vectors object, but without internal array.
1351: Collective
1353: Input Parameter:
1354: . bv - the basis vectors context
1356: Output Parameter:
1357: . v - the new vector
1359: Note:
1360: This works as BVCreateVec(), but the new vector does not have the array allocated,
1361: so the intended usage is with VecPlaceArray().
1363: Level: developer
1365: .seealso: BVCreateVec()
1366: @*/
1367: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1368: {
1369: PetscBool standard,cuda,hip,mpi;
1370: PetscInt N,nloc,bs;
1372: PetscFunctionBegin;
1374: BVCheckSizes(bv,1);
1375: PetscAssertPointer(v,2);
1377: PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1378: PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1379: PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
1380: if (standard || cuda || hip) {
1381: PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
1382: PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1383: PetscCall(PetscLayoutGetSize(bv->map,&N));
1384: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1385: if (cuda) {
1386: #if defined(PETSC_HAVE_CUDA)
1387: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1388: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1389: #endif
1390: } else if (hip) {
1391: #if defined(PETSC_HAVE_HIP)
1392: if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1393: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1394: #endif
1395: } else {
1396: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1397: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1398: }
1399: } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1400: PetscFunctionReturn(PETSC_SUCCESS);
1401: }
1403: /*@
1404: BVSetVecType - Set the vector type to be used when creating vectors via BVCreateVec().
1406: Collective
1408: Input Parameters:
1409: + bv - the basis vectors context
1410: - vtype - the vector type
1412: Level: advanced
1414: Note:
1415: This is not needed if the BV object is set up with BVSetSizesFromVec(), but may be
1416: required in the case of BVSetSizes() if one wants to work with non-standard vectors.
1418: .seealso: BVGetVecType(), BVSetSizesFromVec(), BVSetSizes()
1419: @*/
1420: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1421: {
1422: PetscBool std;
1423: PetscMPIInt size;
1425: PetscFunctionBegin;
1427: PetscCall(PetscFree(bv->vtype));
1428: PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1429: if (std) {
1430: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1431: PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1432: } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@
1437: BVGetVecType - Get the vector type to be used when creating vectors via BVCreateVec().
1439: Not Collective
1441: Input Parameter:
1442: . bv - the basis vectors context
1444: Output Parameter:
1445: . vtype - the vector type
1447: Level: advanced
1449: .seealso: BVSetVecType()
1450: @*/
1451: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1452: {
1453: PetscFunctionBegin;
1455: PetscAssertPointer(vtype,2);
1456: *vtype = bv->vtype;
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@
1461: BVCreateMat - Creates a new Mat object of dense type and copies the contents
1462: of the BV object.
1464: Collective
1466: Input Parameter:
1467: . bv - the basis vectors context
1469: Output Parameter:
1470: . A - the new matrix
1472: Notes:
1473: The user is responsible of destroying the returned matrix.
1475: The matrix contains all columns of the BV, not just the active columns.
1477: Level: intermediate
1479: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1480: @*/
1481: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1482: {
1483: PetscInt ksave,lsave;
1484: Mat B;
1486: PetscFunctionBegin;
1488: BVCheckSizes(bv,1);
1489: PetscAssertPointer(A,2);
1491: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1492: lsave = bv->l;
1493: ksave = bv->k;
1494: bv->l = 0;
1495: bv->k = bv->m;
1496: PetscCall(BVGetMat(bv,&B));
1497: PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1498: PetscCall(BVRestoreMat(bv,&B));
1499: bv->l = lsave;
1500: bv->k = ksave;
1501: PetscFunctionReturn(PETSC_SUCCESS);
1502: }
1504: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1505: {
1506: PetscScalar *vv,*aa;
1507: PetscBool create=PETSC_FALSE;
1508: PetscInt m,cols;
1510: PetscFunctionBegin;
1511: m = bv->k-bv->l;
1512: if (!bv->Aget) create=PETSC_TRUE;
1513: else {
1514: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1515: PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1516: PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1517: if (cols!=m) {
1518: PetscCall(MatDestroy(&bv->Aget));
1519: create=PETSC_TRUE;
1520: }
1521: }
1522: PetscCall(BVGetArray(bv,&vv));
1523: if (create) {
1524: 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 */
1525: PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1526: }
1527: PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld)); /* set the actual pointer */
1528: *A = bv->Aget;
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: /*@
1533: BVGetMat - Returns a Mat object of dense type that shares the memory of
1534: the BV object.
1536: Collective
1538: Input Parameter:
1539: . bv - the basis vectors context
1541: Output Parameter:
1542: . A - the matrix
1544: Notes:
1545: The returned matrix contains only the active columns. If the content of
1546: the Mat is modified, these changes are also done in the BV object. The
1547: user must call BVRestoreMat() when no longer needed.
1549: This operation implies a call to BVGetArray(), which may result in data
1550: copies.
1552: Level: advanced
1554: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1555: @*/
1556: PetscErrorCode BVGetMat(BV bv,Mat *A)
1557: {
1558: char name[64];
1560: PetscFunctionBegin;
1562: BVCheckSizes(bv,1);
1563: PetscAssertPointer(A,2);
1564: PetscUseTypeMethod(bv,getmat,A);
1565: if (((PetscObject)bv)->name) { /* set A's name based on BV name */
1566: PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1567: PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1568: PetscCall(PetscObjectSetName((PetscObject)*A,name));
1569: }
1570: PetscFunctionReturn(PETSC_SUCCESS);
1571: }
1573: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1574: {
1575: PetscScalar *vv,*aa;
1577: PetscFunctionBegin;
1578: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1579: vv = aa-(bv->nc+bv->l)*bv->ld;
1580: PetscCall(MatDenseResetArray(bv->Aget));
1581: PetscCall(BVRestoreArray(bv,&vv));
1582: *A = NULL;
1583: PetscFunctionReturn(PETSC_SUCCESS);
1584: }
1586: /*@
1587: BVRestoreMat - Restores the Mat obtained with BVGetMat().
1589: Logically Collective
1591: Input Parameters:
1592: + bv - the basis vectors context
1593: - A - the fetched matrix
1595: Note:
1596: A call to this function must match a previous call of BVGetMat().
1597: The effect is that the contents of the Mat are copied back to the
1598: BV internal data structures.
1600: Level: advanced
1602: .seealso: BVGetMat(), BVRestoreArray()
1603: @*/
1604: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1605: {
1606: PetscFunctionBegin;
1608: BVCheckSizes(bv,1);
1609: PetscAssertPointer(A,2);
1610: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1611: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1612: PetscUseTypeMethod(bv,restoremat,A);
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: /*
1617: Copy all user-provided attributes of V to another BV object W
1618: */
1619: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1620: {
1621: PetscFunctionBegin;
1622: PetscCall(PetscLayoutReference(V->map,&W->map));
1623: PetscCall(BVSetVecType(W,V->vtype));
1624: W->ld = V->ld;
1625: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1626: W->orthog_type = V->orthog_type;
1627: W->orthog_ref = V->orthog_ref;
1628: W->orthog_eta = V->orthog_eta;
1629: W->orthog_block = V->orthog_block;
1630: if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1631: W->matrix = V->matrix;
1632: W->indef = V->indef;
1633: W->vmm = V->vmm;
1634: W->rrandom = V->rrandom;
1635: W->deftol = V->deftol;
1636: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1637: W->rand = V->rand;
1638: W->sfocalled = V->sfocalled;
1639: PetscTryTypeMethod(V,duplicate,W);
1640: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1641: PetscFunctionReturn(PETSC_SUCCESS);
1642: }
1644: /*@
1645: BVDuplicate - Creates a new basis vector object of the same type and
1646: dimensions as an existing one.
1648: Collective
1650: Input Parameter:
1651: . V - basis vectors context
1653: Output Parameter:
1654: . W - location to put the new BV
1656: Notes:
1657: The new BV has the same type and dimensions as V. Also, the inner
1658: product matrix and orthogonalization options are copied.
1660: BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1661: for the new basis vectors. Use BVCopy() to copy the contents.
1663: Level: intermediate
1665: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1666: @*/
1667: PetscErrorCode BVDuplicate(BV V,BV *W)
1668: {
1669: PetscFunctionBegin;
1672: BVCheckSizes(V,1);
1673: PetscAssertPointer(W,2);
1674: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1675: (*W)->N = V->N;
1676: (*W)->n = V->n;
1677: (*W)->m = V->m;
1678: (*W)->k = V->m;
1679: PetscCall(BVDuplicate_Private(V,*W));
1680: PetscFunctionReturn(PETSC_SUCCESS);
1681: }
1683: /*@
1684: BVDuplicateResize - Creates a new basis vector object of the same type and
1685: dimensions as an existing one, but with possibly different number of columns.
1687: Collective
1689: Input Parameters:
1690: + V - basis vectors context
1691: - m - the new number of columns
1693: Output Parameter:
1694: . W - location to put the new BV
1696: Note:
1697: This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1698: contents of V are not copied to W.
1700: Level: intermediate
1702: .seealso: BVDuplicate(), BVResize()
1703: @*/
1704: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1705: {
1706: PetscFunctionBegin;
1709: BVCheckSizes(V,1);
1711: PetscAssertPointer(W,3);
1712: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1713: (*W)->N = V->N;
1714: (*W)->n = V->n;
1715: (*W)->m = m;
1716: (*W)->k = m;
1717: PetscCall(BVDuplicate_Private(V,*W));
1718: PetscFunctionReturn(PETSC_SUCCESS);
1719: }
1721: /*@
1722: BVGetCachedBV - Returns a BV object stored internally that holds the
1723: result of B*X after a call to BVApplyMatrixBV().
1725: Collective
1727: Input Parameter:
1728: . bv - the basis vectors context
1730: Output Parameter:
1731: . cached - the cached BV
1733: Note:
1734: The cached BV is created if not available previously.
1736: Level: developer
1738: .seealso: BVApplyMatrixBV()
1739: @*/
1740: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1741: {
1742: PetscFunctionBegin;
1744: PetscAssertPointer(cached,2);
1745: BVCheckSizes(bv,1);
1746: if (!bv->cached) {
1747: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1748: bv->cached->N = bv->N;
1749: bv->cached->n = bv->n;
1750: bv->cached->m = bv->m;
1751: bv->cached->k = bv->m;
1752: PetscCall(BVDuplicate_Private(bv,bv->cached));
1753: }
1754: *cached = bv->cached;
1755: PetscFunctionReturn(PETSC_SUCCESS);
1756: }
1758: /*@
1759: BVCopy - Copies a basis vector object into another one, W <- V.
1761: Logically Collective
1763: Input Parameter:
1764: . V - basis vectors context
1766: Output Parameter:
1767: . W - the copy
1769: Note:
1770: Both V and W must be distributed in the same manner; local copies are
1771: done. Only active columns (excluding the leading ones) are copied.
1772: In the destination W, columns are overwritten starting from the leading ones.
1773: Constraints are not copied.
1775: Level: beginner
1777: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1778: @*/
1779: PetscErrorCode BVCopy(BV V,BV W)
1780: {
1781: PetscScalar *womega;
1782: const PetscScalar *vomega;
1784: PetscFunctionBegin;
1787: BVCheckSizes(V,1);
1788: BVCheckOp(V,1,copy);
1791: BVCheckSizes(W,2);
1792: PetscCheckSameTypeAndComm(V,1,W,2);
1793: 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);
1794: 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);
1795: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1797: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1798: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1799: /* copy signature */
1800: PetscCall(BV_AllocateSignature(W));
1801: PetscCall(VecGetArrayRead(V->omega,&vomega));
1802: PetscCall(VecGetArray(W->omega,&womega));
1803: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1804: PetscCall(VecRestoreArray(W->omega,&womega));
1805: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1806: }
1807: PetscUseTypeMethod(V,copy,W);
1808: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1809: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1810: PetscFunctionReturn(PETSC_SUCCESS);
1811: }
1813: /*@
1814: BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1816: Logically Collective
1818: Input Parameters:
1819: + V - basis vectors context
1820: - j - the column number to be copied
1822: Output Parameter:
1823: . w - the copied column
1825: Note:
1826: Both V and w must be distributed in the same manner; local copies are done.
1828: Level: beginner
1830: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1831: @*/
1832: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1833: {
1834: PetscInt n,N;
1835: Vec z;
1837: PetscFunctionBegin;
1840: BVCheckSizes(V,1);
1843: PetscCheckSameComm(V,1,w,3);
1845: PetscCall(VecGetSize(w,&N));
1846: PetscCall(VecGetLocalSize(w,&n));
1847: 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);
1849: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1850: PetscCall(BVGetColumn(V,j,&z));
1851: PetscCall(VecCopy(z,w));
1852: PetscCall(BVRestoreColumn(V,j,&z));
1853: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1854: PetscFunctionReturn(PETSC_SUCCESS);
1855: }
1857: /*@
1858: BVCopyColumn - Copies the values from one of the columns to another one.
1860: Logically Collective
1862: Input Parameters:
1863: + V - basis vectors context
1864: . j - the number of the source column
1865: - i - the number of the destination column
1867: Level: beginner
1869: .seealso: BVCopy(), BVCopyVec()
1870: @*/
1871: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1872: {
1873: PetscScalar *omega;
1875: PetscFunctionBegin;
1878: BVCheckSizes(V,1);
1881: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1883: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1884: if (V->omega) {
1885: PetscCall(VecGetArray(V->omega,&omega));
1886: omega[i] = omega[j];
1887: PetscCall(VecRestoreArray(V->omega,&omega));
1888: }
1889: PetscUseTypeMethod(V,copycolumn,j,i);
1890: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1891: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1892: PetscFunctionReturn(PETSC_SUCCESS);
1893: }
1895: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1896: {
1897: PetscInt ncols;
1899: PetscFunctionBegin;
1900: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1901: if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1902: if (!*split) {
1903: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1904: (*split)->issplit = left? 1: 2;
1905: (*split)->splitparent = bv;
1906: (*split)->N = bv->N;
1907: (*split)->n = bv->n;
1908: (*split)->m = bv->m;
1909: (*split)->k = bv->m;
1910: PetscCall(BVDuplicate_Private(bv,*split));
1911: }
1912: (*split)->l = 0;
1913: (*split)->k = left? bv->l: bv->k-bv->l;
1914: (*split)->nc = left? bv->nc: 0;
1915: (*split)->m = ncols-(*split)->nc;
1916: if ((*split)->nc) {
1917: (*split)->ci[0] = -(*split)->nc-1;
1918: (*split)->ci[1] = -(*split)->nc-1;
1919: }
1920: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1921: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1922: PetscFunctionReturn(PETSC_SUCCESS);
1923: }
1925: /*@
1926: BVGetSplit - Splits the BV object into two BV objects that share the
1927: internal data, one of them containing the leading columns and the other
1928: one containing the remaining columns.
1930: Collective
1932: Input Parameter:
1933: . bv - the basis vectors context
1935: Output Parameters:
1936: + L - left BV containing leading columns (can be NULL)
1937: - R - right BV containing remaining columns (can be NULL)
1939: Notes:
1940: The columns are split in two sets. The leading columns (including the
1941: constraints) are assigned to the left BV and the remaining columns
1942: are assigned to the right BV. The number of leading columns, as
1943: specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1944: guarantee that both L and R have at least one column).
1946: The returned BV's must be seen as references (not copies) of the input
1947: BV, that is, modifying them will change the entries of bv as well.
1948: The returned BV's must not be destroyed. BVRestoreSplit() must be called
1949: when they are no longer needed.
1951: Pass NULL for any of the output BV's that is not needed.
1953: Level: advanced
1955: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints(), BVGetSplitRows()
1956: @*/
1957: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1958: {
1959: PetscFunctionBegin;
1962: BVCheckSizes(bv,1);
1963: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1964: PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1965: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1966: bv->lsplit = bv->nc+bv->l;
1967: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1968: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1969: if (L) *L = bv->L;
1970: if (R) *R = bv->R;
1971: PetscFunctionReturn(PETSC_SUCCESS);
1972: }
1974: /*@
1975: BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1977: Logically Collective
1979: Input Parameters:
1980: + bv - the basis vectors context
1981: . L - left BV obtained with BVGetSplit()
1982: - R - right BV obtained with BVGetSplit()
1984: Note:
1985: The arguments must match the corresponding call to BVGetSplit().
1987: Level: advanced
1989: .seealso: BVGetSplit()
1990: @*/
1991: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1992: {
1993: PetscFunctionBegin;
1996: BVCheckSizes(bv,1);
1997: PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1998: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1999: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
2000: 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()");
2001: 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()");
2003: PetscTryTypeMethod(bv,restoresplit,L,R);
2004: bv->lsplit = 0;
2005: if (L) *L = NULL;
2006: if (R) *R = NULL;
2007: PetscFunctionReturn(PETSC_SUCCESS);
2008: }
2010: /*
2011: Copy all user-provided attributes of V to another BV object W with different layout
2012: */
2013: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
2014: {
2015: PetscFunctionBegin;
2016: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2017: PetscCall(BVSetVecType(W,V->vtype));
2018: W->ld = V->ld;
2019: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2020: W->orthog_type = V->orthog_type;
2021: W->orthog_ref = V->orthog_ref;
2022: W->orthog_eta = V->orthog_eta;
2023: W->orthog_block = V->orthog_block;
2024: W->vmm = V->vmm;
2025: W->rrandom = V->rrandom;
2026: W->deftol = V->deftol;
2027: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
2028: W->rand = V->rand;
2029: W->sfocalled = V->sfocalled;
2030: PetscTryTypeMethod(V,duplicate,W);
2031: PetscCall(PetscObjectStateIncrease((PetscObject)W));
2032: PetscFunctionReturn(PETSC_SUCCESS);
2033: }
2035: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2036: {
2037: PetscInt rstart,rend,lstart;
2038: PetscInt N,n;
2039: PetscBool contig;
2041: PetscFunctionBegin;
2042: PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2043: PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2044: PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2045: if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2046: else bv->lsplit = -lstart;
2047: PetscCall(ISGetSize(is,&N));
2048: PetscCall(ISGetLocalSize(is,&n));
2049: if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2050: if (!*split) {
2051: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2052: (*split)->issplit = top? -1: -2;
2053: (*split)->splitparent = bv;
2054: (*split)->N = N;
2055: (*split)->n = n;
2056: (*split)->m = bv->m;
2057: PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2058: }
2059: (*split)->k = bv->k;
2060: (*split)->l = bv->l;
2061: (*split)->nc = bv->nc;
2062: if ((*split)->nc) {
2063: (*split)->ci[0] = -(*split)->nc-1;
2064: (*split)->ci[1] = -(*split)->nc-1;
2065: }
2066: if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2067: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2068: PetscFunctionReturn(PETSC_SUCCESS);
2069: }
2071: /*@
2072: BVGetSplitRows - Splits the BV object into two BV objects that share the
2073: internal data, using a disjoint horizontal splitting.
2075: Collective
2077: Input Parameters:
2078: + bv - the basis vectors context
2079: . isup - the index set that defines the upper part of the horizontal splitting
2080: - islo - the index set that defines the lower part of the horizontal splitting
2082: Output Parameters:
2083: + U - the resulting BV containing the upper rows
2084: - L - the resulting BV containing the lower rows
2086: Notes:
2087: The index sets must be such that every MPI process can extract the selected
2088: rows from its local part of the input BV, and this part must be contiguous.
2089: With one process, isup will list contiguous indices starting from 0, and islo
2090: will contain the remaining indices, hence we refer to upper and lower part.
2091: However, with several processes the indices will be interleaved because
2092: isup will refer to the upper part of the local array.
2094: The intended use of this function is with matrices of MATNEST type, where
2095: MatNestGetISs() will return the appropriate index sets.
2097: The returned BV's must be seen as references (not copies) of the input
2098: BV, that is, modifying them will change the entries of bv as well.
2099: The returned BV's must not be destroyed. BVRestoreSplitRows() must be called
2100: when they are no longer needed.
2102: Pass NULL for any of the output BV's that is not needed.
2104: Level: advanced
2106: .seealso: BVRestoreSplitRows(), BVGetSplit()
2107: @*/
2108: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2109: {
2110: PetscFunctionBegin;
2113: BVCheckSizes(bv,1);
2116: PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2117: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2118: if (U) {
2119: PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2120: *U = bv->R;
2121: }
2122: if (L) {
2123: PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2124: *L = bv->L;
2125: }
2126: PetscFunctionReturn(PETSC_SUCCESS);
2127: }
2129: /*@
2130: BVRestoreSplitRows - Restore the BV objects obtained with BVGetSplitRows().
2132: Logically Collective
2134: Input Parameters:
2135: + bv - the basis vectors context
2136: . isup - the index set that defines the upper part of the horizontal splitting
2137: . islo - the index set that defines the lower part of the horizontal splitting
2138: . U - the BV containing the upper rows
2139: - L - the BV containing the lower rows
2141: Note:
2142: The arguments must match the corresponding call to BVGetSplitRows().
2144: Level: advanced
2146: .seealso: BVGetSplitRows()
2147: @*/
2148: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2149: {
2150: PetscFunctionBegin;
2153: BVCheckSizes(bv,1);
2154: PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2155: 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()");
2156: 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()");
2157: PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2158: bv->lsplit = 0;
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: /*@
2163: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2164: definite inner product.
2166: Logically Collective
2168: Input Parameters:
2169: + bv - basis vectors
2170: - deftol - the tolerance
2172: Options Database Key:
2173: . -bv_definite_tol <deftol> - the tolerance
2175: Notes:
2176: When using a non-standard inner product, see BVSetMatrix(), the solver needs
2177: to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
2178: been declared indefinite, the value z'*B*z must be positive, but due to
2179: rounding error a tiny value may become negative. A tolerance is used to
2180: detect this situation. Likewise, in complex arithmetic z'*B*z should be
2181: real, and we use the same tolerance to check whether a nonzero imaginary part
2182: can be considered negligible.
2184: This function sets this tolerance, which defaults to 10*eps, where eps is
2185: the machine epsilon. The default value should be good for most applications.
2187: Level: advanced
2189: .seealso: BVSetMatrix()
2190: @*/
2191: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2192: {
2193: PetscFunctionBegin;
2196: if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2197: else {
2198: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2199: bv->deftol = deftol;
2200: }
2201: PetscFunctionReturn(PETSC_SUCCESS);
2202: }
2204: /*@
2205: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2206: inner product.
2208: Not Collective
2210: Input Parameter:
2211: . bv - the basis vectors
2213: Output Parameter:
2214: . deftol - the tolerance
2216: Level: advanced
2218: .seealso: BVSetDefiniteTolerance()
2219: @*/
2220: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2221: {
2222: PetscFunctionBegin;
2224: PetscAssertPointer(deftol,2);
2225: *deftol = bv->deftol;
2226: PetscFunctionReturn(PETSC_SUCCESS);
2227: }
2229: /*@
2230: BVSetLeadingDimension - Set the leading dimension to be used for storing the BV data.
2232: Not Collective
2234: Input Parameters:
2235: + bv - basis vectors
2236: - ld - the leading dimension
2238: Notes:
2239: This parameter is relevant for BVMAT, though it might be employed in other types
2240: as well.
2242: When the internal data of the BV is stored as a dense matrix, the leading dimension
2243: has the same meaning as in MatDenseSetLDA(), i.e., the distance in number of
2244: elements from one entry of the matrix to the one in the next column at the same
2245: row. The leading dimension refers to the local array, and hence can be different
2246: in different processes.
2248: The user does not need to change this parameter. The default value is equal to the
2249: number of local rows, but this value may be increased a little to guarantee alignment
2250: (especially in the case of GPU storage).
2252: Level: advanced
2254: .seealso: BVGetLeadingDimension()
2255: @*/
2256: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2257: {
2258: PetscFunctionBegin;
2261: 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");
2262: if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
2263: else {
2264: PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2265: bv->ld = ld;
2266: }
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }
2270: /*@
2271: BVGetLeadingDimension - Returns the leading dimension of the BV.
2273: Not Collective
2275: Input Parameter:
2276: . bv - the basis vectors
2278: Output Parameter:
2279: . ld - the leading dimension
2281: Level: advanced
2283: Notes:
2284: The returned value may be different in different processes.
2286: The leading dimension must be used when accessing the internal array via
2287: BVGetArray() or BVGetArrayRead().
2289: .seealso: BVSetLeadingDimension(), BVGetArray(), BVGetArrayRead()
2290: @*/
2291: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2292: {
2293: PetscFunctionBegin;
2295: PetscAssertPointer(ld,2);
2296: *ld = bv->ld;
2297: PetscFunctionReturn(PETSC_SUCCESS);
2298: }