Actual source code: bvbasic.c
slepc-main 2024-05-17
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 PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
367: if (copy) {
368: PetscCall(VecGetArray(v,&array));
369: PetscCall(VecGetArrayRead(bv->omega,&omega));
370: PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
371: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
372: PetscCall(VecRestoreArray(v,&array));
373: } else PetscCall(VecSet(v,1.0));
374: PetscCall(VecDestroy(&bv->omega));
375: bv->omega = v;
376: }
377: bv->m = m;
378: bv->k = PetscMin(bv->k,m);
379: bv->l = PetscMin(bv->l,m);
380: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
381: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: BVSetActiveColumns - Specify the columns that will be involved in operations.
388: Logically Collective
390: Input Parameters:
391: + bv - the basis vectors context
392: . l - number of leading columns
393: - k - number of active columns
395: Notes:
396: In operations such as BVMult() or BVDot(), only the first k columns are
397: considered. This is useful when the BV is filled from left to right, so
398: the last m-k columns do not have relevant information.
400: Also in operations such as BVMult() or BVDot(), the first l columns are
401: normally not included in the computation. See the manpage of each
402: operation.
404: In orthogonalization operations, the first l columns are treated
405: differently, they participate in the orthogonalization but the computed
406: coefficients are not stored.
408: Level: intermediate
410: .seealso: BVGetActiveColumns(), BVSetSizes()
411: @*/
412: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
413: {
414: PetscFunctionBegin;
418: BVCheckSizes(bv,1);
419: if (PetscUnlikely(k==PETSC_DECIDE || k==PETSC_DEFAULT)) {
420: bv->k = bv->m;
421: } else {
422: 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);
423: bv->k = k;
424: }
425: if (PetscUnlikely(l==PETSC_DECIDE || l==PETSC_DEFAULT)) {
426: bv->l = 0;
427: } else {
428: 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);
429: bv->l = l;
430: }
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: BVGetActiveColumns - Returns the current active dimensions.
437: Not Collective
439: Input Parameter:
440: . bv - the basis vectors context
442: Output Parameters:
443: + l - number of leading columns
444: - k - number of active columns
446: Level: intermediate
448: .seealso: BVSetActiveColumns()
449: @*/
450: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
451: {
452: PetscFunctionBegin;
454: if (l) *l = bv->l;
455: if (k) *k = bv->k;
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /*@
460: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
462: Collective
464: Input Parameters:
465: + bv - the basis vectors context
466: . B - a symmetric matrix (may be NULL)
467: - indef - a flag indicating if the matrix is indefinite
469: Notes:
470: This is used to specify a non-standard inner product, whose matrix
471: representation is given by B. Then, all inner products required during
472: orthogonalization are computed as (x,y)_B=y^H*B*x rather than the
473: standard form (x,y)=y^H*x.
475: Matrix B must be real symmetric (or complex Hermitian). A genuine inner
476: product requires that B is also positive (semi-)definite. However, we
477: also allow for an indefinite B (setting indef=PETSC_TRUE), in which
478: case the orthogonalization uses an indefinite inner product.
480: This affects operations BVDot(), BVNorm(), BVOrthogonalize(), and variants.
482: Setting B=NULL has the same effect as if the identity matrix was passed.
484: Level: advanced
486: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
487: @*/
488: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
489: {
490: PetscInt m,n;
492: PetscFunctionBegin;
495: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
496: if (B) {
498: PetscCall(MatGetLocalSize(B,&m,&n));
499: PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
500: 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);
501: }
502: if (B) PetscCall(PetscObjectReference((PetscObject)B));
503: PetscCall(MatDestroy(&bv->matrix));
504: bv->matrix = B;
505: bv->indef = indef;
506: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
507: if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
508: if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
509: }
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: /*@
514: BVGetMatrix - Retrieves the matrix representation of the inner product.
516: Not Collective
518: Input Parameter:
519: . bv - the basis vectors context
521: Output Parameters:
522: + B - the matrix of the inner product (may be NULL)
523: - indef - the flag indicating if the matrix is indefinite
525: Level: advanced
527: .seealso: BVSetMatrix()
528: @*/
529: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
530: {
531: PetscFunctionBegin;
533: if (B) *B = bv->matrix;
534: if (indef) *indef = bv->indef;
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: /*@
539: BVApplyMatrix - Multiplies a vector by the matrix representation of the
540: inner product.
542: Neighbor-wise Collective
544: Input Parameters:
545: + bv - the basis vectors context
546: - x - the vector
548: Output Parameter:
549: . y - the result
551: Note:
552: If no matrix was specified this function copies the vector.
554: Level: advanced
556: .seealso: BVSetMatrix(), BVApplyMatrixBV()
557: @*/
558: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
559: {
560: PetscFunctionBegin;
564: if (bv->matrix) {
565: PetscCall(BV_IPMatMult(bv,x));
566: PetscCall(VecCopy(bv->Bx,y));
567: } else PetscCall(VecCopy(x,y));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: /*@
572: BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
573: of the inner product.
575: Neighbor-wise Collective
577: Input Parameter:
578: . X - the basis vectors context
580: Output Parameter:
581: . Y - the basis vectors to store the result (optional)
583: Note:
584: This function computes Y = B*X, where B is the matrix given with
585: BVSetMatrix(). This operation is computed as in BVMatMult().
586: If no matrix was specified, then it just copies Y = X.
588: If no Y is given, the result is stored internally in the cached BV.
590: Level: developer
592: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
593: @*/
594: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
595: {
596: PetscFunctionBegin;
598: if (Y) {
600: if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
601: else PetscCall(BVCopy(X,Y));
602: } else PetscCall(BV_IPMatMultBV(X));
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /*@
607: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
609: Logically Collective
611: Input Parameters:
612: + bv - the basis vectors context
613: - omega - a vector representing the diagonal of the signature matrix
615: Note:
616: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
618: Level: developer
620: .seealso: BVSetMatrix(), BVGetSignature()
621: @*/
622: PetscErrorCode BVSetSignature(BV bv,Vec omega)
623: {
624: PetscInt i,n;
625: const PetscScalar *pomega;
626: PetscScalar *intern;
628: PetscFunctionBegin;
630: BVCheckSizes(bv,1);
634: PetscCall(VecGetSize(omega,&n));
635: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
636: PetscCall(BV_AllocateSignature(bv));
637: if (bv->indef) {
638: PetscCall(VecGetArrayRead(omega,&pomega));
639: PetscCall(VecGetArray(bv->omega,&intern));
640: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
641: PetscCall(VecRestoreArray(bv->omega,&intern));
642: PetscCall(VecRestoreArrayRead(omega,&pomega));
643: } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
644: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: /*@
649: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
651: Not Collective
653: Input Parameter:
654: . bv - the basis vectors context
656: Output Parameter:
657: . omega - a vector representing the diagonal of the signature matrix
659: Note:
660: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
662: Level: developer
664: .seealso: BVSetMatrix(), BVSetSignature()
665: @*/
666: PetscErrorCode BVGetSignature(BV bv,Vec omega)
667: {
668: PetscInt i,n;
669: PetscScalar *pomega;
670: const PetscScalar *intern;
672: PetscFunctionBegin;
674: BVCheckSizes(bv,1);
678: PetscCall(VecGetSize(omega,&n));
679: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
680: if (bv->indef && bv->omega) {
681: PetscCall(VecGetArray(omega,&pomega));
682: PetscCall(VecGetArrayRead(bv->omega,&intern));
683: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
684: PetscCall(VecRestoreArrayRead(bv->omega,&intern));
685: PetscCall(VecRestoreArray(omega,&pomega));
686: } else PetscCall(VecSet(omega,1.0));
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: /*@
691: BVSetBufferVec - Attach a vector object to be used as buffer space for
692: several operations.
694: Collective
696: Input Parameters:
697: + bv - the basis vectors context)
698: - buffer - the vector
700: Notes:
701: Use BVGetBufferVec() to retrieve the vector (for example, to free it
702: at the end of the computations).
704: The vector must be sequential of length (nc+m)*m, where m is the number
705: of columns of bv and nc is the number of constraints.
707: Level: developer
709: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
710: @*/
711: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
712: {
713: PetscInt ld,n;
714: PetscMPIInt size;
716: PetscFunctionBegin;
719: BVCheckSizes(bv,1);
720: PetscCall(VecGetSize(buffer,&n));
721: ld = bv->m+bv->nc;
722: PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
723: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
724: PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
726: PetscCall(PetscObjectReference((PetscObject)buffer));
727: PetscCall(VecDestroy(&bv->buffer));
728: bv->buffer = buffer;
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: /*@
733: BVGetBufferVec - Obtain the buffer vector associated with the BV object.
735: Collective
737: Input Parameters:
738: . bv - the basis vectors context
740: Output Parameter:
741: . buffer - vector
743: Notes:
744: The vector is created if not available previously. It is a sequential vector
745: of length (nc+m)*m, where m is the number of columns of bv and nc is the number
746: of constraints.
748: Developer Notes:
749: The buffer vector is viewed as a column-major matrix with leading dimension
750: ld=nc+m, and m columns at most. In the most common usage, it has the structure
751: .vb
752: | | C |
753: |s|---|
754: | | H |
755: .ve
756: where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
757: related to orthogonalization against constraints (first nc rows), and s is the
758: first column that contains scratch values computed during Gram-Schmidt
759: orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
760: store the coefficients.
762: Level: developer
764: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
765: @*/
766: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
767: {
768: PetscInt ld;
770: PetscFunctionBegin;
772: PetscAssertPointer(buffer,2);
773: BVCheckSizes(bv,1);
774: if (!bv->buffer) {
775: ld = bv->m+bv->nc;
776: PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
777: PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
778: PetscCall(VecSetType(bv->buffer,bv->vtype));
779: }
780: *buffer = bv->buffer;
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: /*@
785: BVSetRandomContext - Sets the PetscRandom object associated with the BV,
786: to be used in operations that need random numbers.
788: Collective
790: Input Parameters:
791: + bv - the basis vectors context
792: - rand - the random number generator context
794: Level: advanced
796: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
797: @*/
798: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
799: {
800: PetscFunctionBegin;
803: PetscCheckSameComm(bv,1,rand,2);
804: PetscCall(PetscObjectReference((PetscObject)rand));
805: PetscCall(PetscRandomDestroy(&bv->rand));
806: bv->rand = rand;
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /*@
811: BVGetRandomContext - Gets the PetscRandom object associated with the BV.
813: Collective
815: Input Parameter:
816: . bv - the basis vectors context
818: Output Parameter:
819: . rand - the random number generator context
821: Level: advanced
823: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
824: @*/
825: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
826: {
827: PetscFunctionBegin;
829: PetscAssertPointer(rand,2);
830: if (!bv->rand) {
831: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
832: if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
833: if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
834: if (bv->rrandom) {
835: PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
836: PetscCall(PetscRandomSeed(bv->rand));
837: }
838: }
839: *rand = bv->rand;
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: /*@
844: BVSetFromOptions - Sets BV options from the options database.
846: Collective
848: Input Parameter:
849: . bv - the basis vectors context
851: Level: beginner
853: .seealso: BVSetOptionsPrefix()
854: @*/
855: PetscErrorCode BVSetFromOptions(BV bv)
856: {
857: char type[256];
858: PetscBool flg1,flg2,flg3,flg4;
859: PetscReal r;
860: BVOrthogType otype;
861: BVOrthogRefineType orefine;
862: BVOrthogBlockType oblock;
864: PetscFunctionBegin;
866: PetscCall(BVRegisterAll());
867: PetscObjectOptionsBegin((PetscObject)bv);
868: PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
869: if (flg1) PetscCall(BVSetType(bv,type));
870: else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));
872: otype = bv->orthog_type;
873: PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
874: orefine = bv->orthog_ref;
875: PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
876: r = bv->orthog_eta;
877: PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
878: oblock = bv->orthog_block;
879: PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
880: if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
882: PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
884: PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
885: if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
887: /* undocumented option to generate random vectors that are independent of the number of processes */
888: PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
890: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
891: else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
892: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
893: PetscOptionsEnd();
894: bv->sfocalled = PETSC_TRUE;
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: /*@
899: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
900: vectors (classical or modified Gram-Schmidt with or without refinement), and
901: for the block-orthogonalization (simultaneous orthogonalization of a set of
902: vectors).
904: Logically Collective
906: Input Parameters:
907: + bv - the basis vectors context
908: . type - the method of vector orthogonalization
909: . refine - type of refinement
910: . eta - parameter for selective refinement
911: - block - the method of block orthogonalization
913: Options Database Keys:
914: + -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
915: (default) or mgs for Modified Gram-Schmidt orthogonalization
916: . -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
917: . -bv_orthog_eta <eta> - For setting the value of eta
918: - -bv_orthog_block <block> - Where <block> is the block-orthogonalization method
920: Notes:
921: The default settings work well for most problems.
923: The parameter eta should be a real value between 0 and 1 (or PETSC_DEFAULT).
924: The value of eta is used only when the refinement type is "ifneeded".
926: When using several processors, MGS is likely to result in bad scalability.
928: If the method set for block orthogonalization is GS, then the computation
929: is done column by column with the vector orthogonalization.
931: Level: advanced
933: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
934: @*/
935: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
936: {
937: PetscFunctionBegin;
943: switch (type) {
944: case BV_ORTHOG_CGS:
945: case BV_ORTHOG_MGS:
946: bv->orthog_type = type;
947: break;
948: default:
949: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
950: }
951: switch (refine) {
952: case BV_ORTHOG_REFINE_NEVER:
953: case BV_ORTHOG_REFINE_IFNEEDED:
954: case BV_ORTHOG_REFINE_ALWAYS:
955: bv->orthog_ref = refine;
956: break;
957: default:
958: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
959: }
960: if (eta == (PetscReal)PETSC_DEFAULT) {
961: bv->orthog_eta = 0.7071;
962: } else {
963: PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
964: bv->orthog_eta = eta;
965: }
966: switch (block) {
967: case BV_ORTHOG_BLOCK_GS:
968: case BV_ORTHOG_BLOCK_CHOL:
969: case BV_ORTHOG_BLOCK_TSQR:
970: case BV_ORTHOG_BLOCK_TSQRCHOL:
971: case BV_ORTHOG_BLOCK_SVQB:
972: bv->orthog_block = block;
973: break;
974: default:
975: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
976: }
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@
981: BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.
983: Not Collective
985: Input Parameter:
986: . bv - basis vectors context
988: Output Parameters:
989: + type - the method of vector orthogonalization
990: . refine - type of refinement
991: . eta - parameter for selective refinement
992: - block - the method of block orthogonalization
994: Level: advanced
996: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
997: @*/
998: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
999: {
1000: PetscFunctionBegin;
1002: if (type) *type = bv->orthog_type;
1003: if (refine) *refine = bv->orthog_ref;
1004: if (eta) *eta = bv->orthog_eta;
1005: if (block) *block = bv->orthog_block;
1006: PetscFunctionReturn(PETSC_SUCCESS);
1007: }
1009: /*@
1010: BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.
1012: Logically Collective
1014: Input Parameters:
1015: + bv - the basis vectors context
1016: - method - the method for the BVMatMult() operation
1018: Options Database Keys:
1019: . -bv_matmult <meth> - choose one of the methods: vecs, mat
1021: Notes:
1022: Allowed values are
1023: + BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1024: . BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1025: - BV_MATMULT_MAT_SAVE - this case is deprecated
1027: The default is BV_MATMULT_MAT except in the case of BVVECS.
1029: Level: advanced
1031: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1032: @*/
1033: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1034: {
1035: PetscFunctionBegin;
1038: switch (method) {
1039: case BV_MATMULT_VECS:
1040: case BV_MATMULT_MAT:
1041: bv->vmm = method;
1042: break;
1043: case BV_MATMULT_MAT_SAVE:
1044: PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1045: bv->vmm = BV_MATMULT_MAT;
1046: break;
1047: default:
1048: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1049: }
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*@
1054: BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.
1056: Not Collective
1058: Input Parameter:
1059: . bv - basis vectors context
1061: Output Parameter:
1062: . method - the method for the BVMatMult() operation
1064: Level: advanced
1066: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1067: @*/
1068: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1069: {
1070: PetscFunctionBegin;
1072: PetscAssertPointer(method,2);
1073: *method = bv->vmm;
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@
1078: BVGetColumn - Returns a Vec object that contains the entries of the
1079: requested column of the basis vectors object.
1081: Logically Collective
1083: Input Parameters:
1084: + bv - the basis vectors context
1085: - j - the index of the requested column
1087: Output Parameter:
1088: . v - vector containing the jth column
1090: Notes:
1091: The returned Vec must be seen as a reference (not a copy) of the BV
1092: column, that is, modifying the Vec will change the BV entries as well.
1094: The returned Vec must not be destroyed. BVRestoreColumn() must be
1095: called when it is no longer needed. At most, two columns can be fetched,
1096: that is, this function can only be called twice before the corresponding
1097: BVRestoreColumn() is invoked.
1099: A negative index j selects the i-th constraint, where i=-j. Constraints
1100: should not be modified.
1102: Level: beginner
1104: .seealso: BVRestoreColumn(), BVInsertConstraints()
1105: @*/
1106: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1107: {
1108: PetscInt l;
1110: PetscFunctionBegin;
1113: BVCheckSizes(bv,1);
1114: BVCheckOp(bv,1,getcolumn);
1116: 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);
1117: 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);
1118: 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);
1119: l = BVAvailableVec;
1120: PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1121: PetscUseTypeMethod(bv,getcolumn,j,v);
1122: bv->ci[l] = j;
1123: PetscCall(PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]));
1124: PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1125: *v = bv->cv[l];
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1129: /*@
1130: BVRestoreColumn - Restore a column obtained with BVGetColumn().
1132: Logically Collective
1134: Input Parameters:
1135: + bv - the basis vectors context
1136: . j - the index of the column
1137: - v - vector obtained with BVGetColumn()
1139: Note:
1140: The arguments must match the corresponding call to BVGetColumn().
1142: Level: beginner
1144: .seealso: BVGetColumn()
1145: @*/
1146: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1147: {
1148: PetscObjectId id;
1149: PetscObjectState st;
1150: PetscInt l;
1152: PetscFunctionBegin;
1155: BVCheckSizes(bv,1);
1157: PetscAssertPointer(v,3);
1159: 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);
1160: 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);
1161: 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);
1162: l = (j==bv->ci[0])? 0: 1;
1163: PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1164: PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1165: PetscCall(PetscObjectStateGet((PetscObject)*v,&st));
1166: if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1167: PetscUseTypeMethod(bv,restorecolumn,j,v);
1168: bv->ci[l] = -bv->nc-1;
1169: bv->st[l] = -1;
1170: bv->id[l] = 0;
1171: *v = NULL;
1172: PetscFunctionReturn(PETSC_SUCCESS);
1173: }
1175: /*@C
1176: BVGetArray - Returns a pointer to a contiguous array that contains this
1177: processor's portion of the BV data.
1179: Logically Collective
1181: Input Parameters:
1182: . bv - the basis vectors context
1184: Output Parameter:
1185: . a - location to put pointer to the array
1187: Notes:
1188: BVRestoreArray() must be called when access to the array is no longer needed.
1189: This operation may imply a data copy, for BV types that do not store
1190: data contiguously in memory.
1192: The pointer will normally point to the first entry of the first column,
1193: but if the BV has constraints then these go before the regular columns.
1195: Note that for manipulating the pointer to the BV array, one must take into
1196: account the leading dimension, which might be different from the local
1197: number of rows, see BVGetLeadingDimension().
1199: Use BVGetArrayRead() for read-only access.
1201: Level: advanced
1203: .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArrayRead()
1204: @*/
1205: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1206: {
1207: PetscFunctionBegin;
1210: BVCheckSizes(bv,1);
1211: BVCheckOp(bv,1,getarray);
1212: PetscUseTypeMethod(bv,getarray,a);
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }
1216: /*@C
1217: BVRestoreArray - Restore the BV object after BVGetArray() has been called.
1219: Logically Collective
1221: Input Parameters:
1222: + bv - the basis vectors context
1223: - a - location of pointer to array obtained from BVGetArray()
1225: Note:
1226: This operation may imply a data copy, for BV types that do not store
1227: data contiguously in memory.
1229: Level: advanced
1231: .seealso: BVGetColumn()
1232: @*/
1233: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1234: {
1235: PetscFunctionBegin;
1238: BVCheckSizes(bv,1);
1239: PetscTryTypeMethod(bv,restorearray,a);
1240: if (a) *a = NULL;
1241: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1242: PetscFunctionReturn(PETSC_SUCCESS);
1243: }
1245: /*@C
1246: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1247: contains this processor's portion of the BV data.
1249: Not Collective
1251: Input Parameters:
1252: . bv - the basis vectors context
1254: Output Parameter:
1255: . a - location to put pointer to the array
1257: Notes:
1258: BVRestoreArrayRead() must be called when access to the array is no
1259: longer needed. This operation may imply a data copy, for BV types that
1260: do not store data contiguously in memory.
1262: The pointer will normally point to the first entry of the first column,
1263: but if the BV has constraints then these go before the regular columns.
1265: Level: advanced
1267: .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArray()
1268: @*/
1269: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1270: {
1271: PetscFunctionBegin;
1274: BVCheckSizes(bv,1);
1275: BVCheckOp(bv,1,getarrayread);
1276: PetscUseTypeMethod(bv,getarrayread,a);
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: /*@C
1281: BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1282: been called.
1284: Not Collective
1286: Input Parameters:
1287: + bv - the basis vectors context
1288: - a - location of pointer to array obtained from BVGetArrayRead()
1290: Level: advanced
1292: .seealso: BVGetColumn()
1293: @*/
1294: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1295: {
1296: PetscFunctionBegin;
1299: BVCheckSizes(bv,1);
1300: PetscTryTypeMethod(bv,restorearrayread,a);
1301: if (a) *a = NULL;
1302: PetscFunctionReturn(PETSC_SUCCESS);
1303: }
1305: /*@
1306: BVCreateVec - Creates a new Vec object with the same type and dimensions
1307: as the columns of the basis vectors object.
1309: Collective
1311: Input Parameter:
1312: . bv - the basis vectors context
1314: Output Parameter:
1315: . v - the new vector
1317: Note:
1318: The user is responsible of destroying the returned vector.
1320: Level: beginner
1322: .seealso: BVCreateMat(), BVCreateVecEmpty()
1323: @*/
1324: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1325: {
1326: PetscFunctionBegin;
1328: BVCheckSizes(bv,1);
1329: PetscAssertPointer(v,2);
1330: PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1331: PetscCall(VecSetLayout(*v,bv->map));
1332: PetscCall(VecSetType(*v,bv->vtype));
1333: PetscCall(VecSetUp(*v));
1334: PetscFunctionReturn(PETSC_SUCCESS);
1335: }
1337: /*@
1338: BVCreateVecEmpty - Creates a new Vec object with the same type and dimensions
1339: as the columns of the basis vectors object, but without internal array.
1341: Collective
1343: Input Parameter:
1344: . bv - the basis vectors context
1346: Output Parameter:
1347: . v - the new vector
1349: Note:
1350: This works as BVCreateVec(), but the new vector does not have the array allocated,
1351: so the intended usage is with VecPlaceArray().
1353: Level: developer
1355: .seealso: BVCreateVec()
1356: @*/
1357: PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1358: {
1359: PetscBool standard,cuda,mpi;
1360: PetscInt N,nloc,bs;
1362: PetscFunctionBegin;
1364: BVCheckSizes(bv,1);
1365: PetscAssertPointer(v,2);
1367: PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1368: PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1369: if (standard || cuda) {
1370: PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,""));
1371: PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1372: PetscCall(PetscLayoutGetSize(bv->map,&N));
1373: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1374: if (cuda) {
1375: #if defined(PETSC_HAVE_CUDA)
1376: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1377: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1378: #endif
1379: } else {
1380: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1381: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1382: }
1383: } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1384: PetscFunctionReturn(PETSC_SUCCESS);
1385: }
1387: /*@
1388: BVSetVecType - Set the vector type to be used when creating vectors via BVCreateVec().
1390: Collective
1392: Input Parameters:
1393: + bv - the basis vectors context
1394: - vtype - the vector type
1396: Level: advanced
1398: Note:
1399: This is not needed if the BV object is set up with BVSetSizesFromVec(), but may be
1400: required in the case of BVSetSizes() if one wants to work with non-standard vectors.
1402: .seealso: BVGetVecType(), BVSetSizesFromVec(), BVSetSizes()
1403: @*/
1404: PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1405: {
1406: PetscBool std;
1407: PetscMPIInt size;
1409: PetscFunctionBegin;
1411: PetscCall(PetscFree(bv->vtype));
1412: PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1413: if (std) {
1414: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1415: PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1416: } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1417: PetscFunctionReturn(PETSC_SUCCESS);
1418: }
1420: /*@
1421: BVGetVecType - Get the vector type to be used when creating vectors via BVCreateVec().
1423: Not Collective
1425: Input Parameter:
1426: . bv - the basis vectors context
1428: Output Parameter:
1429: . vtype - the vector type
1431: Level: advanced
1433: .seealso: BVSetVecType()
1434: @*/
1435: PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1436: {
1437: PetscFunctionBegin;
1439: PetscAssertPointer(vtype,2);
1440: *vtype = bv->vtype;
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: /*@
1445: BVCreateMat - Creates a new Mat object of dense type and copies the contents
1446: of the BV object.
1448: Collective
1450: Input Parameter:
1451: . bv - the basis vectors context
1453: Output Parameter:
1454: . A - the new matrix
1456: Notes:
1457: The user is responsible of destroying the returned matrix.
1459: The matrix contains all columns of the BV, not just the active columns.
1461: Level: intermediate
1463: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1464: @*/
1465: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1466: {
1467: PetscInt ksave,lsave;
1468: Mat B;
1470: PetscFunctionBegin;
1472: BVCheckSizes(bv,1);
1473: PetscAssertPointer(A,2);
1475: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1476: lsave = bv->l;
1477: ksave = bv->k;
1478: bv->l = 0;
1479: bv->k = bv->m;
1480: PetscCall(BVGetMat(bv,&B));
1481: PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1482: PetscCall(BVRestoreMat(bv,&B));
1483: bv->l = lsave;
1484: bv->k = ksave;
1485: PetscFunctionReturn(PETSC_SUCCESS);
1486: }
1488: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1489: {
1490: PetscScalar *vv,*aa;
1491: PetscBool create=PETSC_FALSE;
1492: PetscInt m,cols;
1494: PetscFunctionBegin;
1495: m = bv->k-bv->l;
1496: if (!bv->Aget) create=PETSC_TRUE;
1497: else {
1498: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1499: PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1500: PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1501: if (cols!=m) {
1502: PetscCall(MatDestroy(&bv->Aget));
1503: create=PETSC_TRUE;
1504: }
1505: }
1506: PetscCall(BVGetArray(bv,&vv));
1507: if (create) {
1508: 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 */
1509: PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1510: }
1511: PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld)); /* set the actual pointer */
1512: *A = bv->Aget;
1513: PetscFunctionReturn(PETSC_SUCCESS);
1514: }
1516: /*@
1517: BVGetMat - Returns a Mat object of dense type that shares the memory of
1518: the BV object.
1520: Collective
1522: Input Parameter:
1523: . bv - the basis vectors context
1525: Output Parameter:
1526: . A - the matrix
1528: Notes:
1529: The returned matrix contains only the active columns. If the content of
1530: the Mat is modified, these changes are also done in the BV object. The
1531: user must call BVRestoreMat() when no longer needed.
1533: This operation implies a call to BVGetArray(), which may result in data
1534: copies.
1536: Level: advanced
1538: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1539: @*/
1540: PetscErrorCode BVGetMat(BV bv,Mat *A)
1541: {
1542: char name[64];
1544: PetscFunctionBegin;
1546: BVCheckSizes(bv,1);
1547: PetscAssertPointer(A,2);
1548: PetscUseTypeMethod(bv,getmat,A);
1549: if (((PetscObject)bv)->name) { /* set A's name based on BV name */
1550: PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1551: PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1552: PetscCall(PetscObjectSetName((PetscObject)*A,name));
1553: }
1554: PetscFunctionReturn(PETSC_SUCCESS);
1555: }
1557: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1558: {
1559: PetscScalar *vv,*aa;
1561: PetscFunctionBegin;
1562: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1563: vv = aa-(bv->nc+bv->l)*bv->ld;
1564: PetscCall(MatDenseResetArray(bv->Aget));
1565: PetscCall(BVRestoreArray(bv,&vv));
1566: *A = NULL;
1567: PetscFunctionReturn(PETSC_SUCCESS);
1568: }
1570: /*@
1571: BVRestoreMat - Restores the Mat obtained with BVGetMat().
1573: Logically Collective
1575: Input Parameters:
1576: + bv - the basis vectors context
1577: - A - the fetched matrix
1579: Note:
1580: A call to this function must match a previous call of BVGetMat().
1581: The effect is that the contents of the Mat are copied back to the
1582: BV internal data structures.
1584: Level: advanced
1586: .seealso: BVGetMat(), BVRestoreArray()
1587: @*/
1588: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1589: {
1590: PetscFunctionBegin;
1592: BVCheckSizes(bv,1);
1593: PetscAssertPointer(A,2);
1594: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1595: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1596: PetscUseTypeMethod(bv,restoremat,A);
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /*
1601: Copy all user-provided attributes of V to another BV object W
1602: */
1603: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1604: {
1605: PetscFunctionBegin;
1606: PetscCall(PetscLayoutReference(V->map,&W->map));
1607: PetscCall(BVSetVecType(W,V->vtype));
1608: W->ld = V->ld;
1609: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1610: W->orthog_type = V->orthog_type;
1611: W->orthog_ref = V->orthog_ref;
1612: W->orthog_eta = V->orthog_eta;
1613: W->orthog_block = V->orthog_block;
1614: if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1615: W->matrix = V->matrix;
1616: W->indef = V->indef;
1617: W->vmm = V->vmm;
1618: W->rrandom = V->rrandom;
1619: W->deftol = V->deftol;
1620: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1621: W->rand = V->rand;
1622: W->sfocalled = V->sfocalled;
1623: PetscTryTypeMethod(V,duplicate,W);
1624: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1625: PetscFunctionReturn(PETSC_SUCCESS);
1626: }
1628: /*@
1629: BVDuplicate - Creates a new basis vector object of the same type and
1630: dimensions as an existing one.
1632: Collective
1634: Input Parameter:
1635: . V - basis vectors context
1637: Output Parameter:
1638: . W - location to put the new BV
1640: Notes:
1641: The new BV has the same type and dimensions as V. Also, the inner
1642: product matrix and orthogonalization options are copied.
1644: BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1645: for the new basis vectors. Use BVCopy() to copy the contents.
1647: Level: intermediate
1649: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1650: @*/
1651: PetscErrorCode BVDuplicate(BV V,BV *W)
1652: {
1653: PetscFunctionBegin;
1656: BVCheckSizes(V,1);
1657: PetscAssertPointer(W,2);
1658: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1659: (*W)->N = V->N;
1660: (*W)->n = V->n;
1661: (*W)->m = V->m;
1662: (*W)->k = V->m;
1663: PetscCall(BVDuplicate_Private(V,*W));
1664: PetscFunctionReturn(PETSC_SUCCESS);
1665: }
1667: /*@
1668: BVDuplicateResize - Creates a new basis vector object of the same type and
1669: dimensions as an existing one, but with possibly different number of columns.
1671: Collective
1673: Input Parameters:
1674: + V - basis vectors context
1675: - m - the new number of columns
1677: Output Parameter:
1678: . W - location to put the new BV
1680: Note:
1681: This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1682: contents of V are not copied to W.
1684: Level: intermediate
1686: .seealso: BVDuplicate(), BVResize()
1687: @*/
1688: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1689: {
1690: PetscFunctionBegin;
1693: BVCheckSizes(V,1);
1695: PetscAssertPointer(W,3);
1696: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1697: (*W)->N = V->N;
1698: (*W)->n = V->n;
1699: (*W)->m = m;
1700: (*W)->k = m;
1701: PetscCall(BVDuplicate_Private(V,*W));
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: /*@
1706: BVGetCachedBV - Returns a BV object stored internally that holds the
1707: result of B*X after a call to BVApplyMatrixBV().
1709: Collective
1711: Input Parameter:
1712: . bv - the basis vectors context
1714: Output Parameter:
1715: . cached - the cached BV
1717: Note:
1718: The cached BV is created if not available previously.
1720: Level: developer
1722: .seealso: BVApplyMatrixBV()
1723: @*/
1724: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1725: {
1726: PetscFunctionBegin;
1728: PetscAssertPointer(cached,2);
1729: BVCheckSizes(bv,1);
1730: if (!bv->cached) {
1731: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1732: bv->cached->N = bv->N;
1733: bv->cached->n = bv->n;
1734: bv->cached->m = bv->m;
1735: bv->cached->k = bv->m;
1736: PetscCall(BVDuplicate_Private(bv,bv->cached));
1737: }
1738: *cached = bv->cached;
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: /*@
1743: BVCopy - Copies a basis vector object into another one, W <- V.
1745: Logically Collective
1747: Input Parameter:
1748: . V - basis vectors context
1750: Output Parameter:
1751: . W - the copy
1753: Note:
1754: Both V and W must be distributed in the same manner; local copies are
1755: done. Only active columns (excluding the leading ones) are copied.
1756: In the destination W, columns are overwritten starting from the leading ones.
1757: Constraints are not copied.
1759: Level: beginner
1761: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1762: @*/
1763: PetscErrorCode BVCopy(BV V,BV W)
1764: {
1765: PetscScalar *womega;
1766: const PetscScalar *vomega;
1768: PetscFunctionBegin;
1771: BVCheckSizes(V,1);
1772: BVCheckOp(V,1,copy);
1775: BVCheckSizes(W,2);
1776: PetscCheckSameTypeAndComm(V,1,W,2);
1777: 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);
1778: 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);
1779: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1781: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1782: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1783: /* copy signature */
1784: PetscCall(BV_AllocateSignature(W));
1785: PetscCall(VecGetArrayRead(V->omega,&vomega));
1786: PetscCall(VecGetArray(W->omega,&womega));
1787: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1788: PetscCall(VecRestoreArray(W->omega,&womega));
1789: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1790: }
1791: PetscUseTypeMethod(V,copy,W);
1792: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1793: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: /*@
1798: BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1800: Logically Collective
1802: Input Parameters:
1803: + V - basis vectors context
1804: - j - the column number to be copied
1806: Output Parameter:
1807: . w - the copied column
1809: Note:
1810: Both V and w must be distributed in the same manner; local copies are done.
1812: Level: beginner
1814: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1815: @*/
1816: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1817: {
1818: PetscInt n,N;
1819: Vec z;
1821: PetscFunctionBegin;
1824: BVCheckSizes(V,1);
1827: PetscCheckSameComm(V,1,w,3);
1829: PetscCall(VecGetSize(w,&N));
1830: PetscCall(VecGetLocalSize(w,&n));
1831: 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);
1833: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1834: PetscCall(BVGetColumn(V,j,&z));
1835: PetscCall(VecCopy(z,w));
1836: PetscCall(BVRestoreColumn(V,j,&z));
1837: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1838: PetscFunctionReturn(PETSC_SUCCESS);
1839: }
1841: /*@
1842: BVCopyColumn - Copies the values from one of the columns to another one.
1844: Logically Collective
1846: Input Parameters:
1847: + V - basis vectors context
1848: . j - the number of the source column
1849: - i - the number of the destination column
1851: Level: beginner
1853: .seealso: BVCopy(), BVCopyVec()
1854: @*/
1855: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1856: {
1857: PetscScalar *omega;
1859: PetscFunctionBegin;
1862: BVCheckSizes(V,1);
1865: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1867: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1868: if (V->omega) {
1869: PetscCall(VecGetArray(V->omega,&omega));
1870: omega[i] = omega[j];
1871: PetscCall(VecRestoreArray(V->omega,&omega));
1872: }
1873: PetscUseTypeMethod(V,copycolumn,j,i);
1874: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1875: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1876: PetscFunctionReturn(PETSC_SUCCESS);
1877: }
1879: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1880: {
1881: PetscInt ncols;
1883: PetscFunctionBegin;
1884: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1885: if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1886: if (!*split) {
1887: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1888: (*split)->issplit = left? 1: 2;
1889: (*split)->splitparent = bv;
1890: (*split)->N = bv->N;
1891: (*split)->n = bv->n;
1892: (*split)->m = bv->m;
1893: (*split)->k = bv->m;
1894: PetscCall(BVDuplicate_Private(bv,*split));
1895: }
1896: (*split)->l = 0;
1897: (*split)->k = left? bv->l: bv->k-bv->l;
1898: (*split)->nc = left? bv->nc: 0;
1899: (*split)->m = ncols-(*split)->nc;
1900: if ((*split)->nc) {
1901: (*split)->ci[0] = -(*split)->nc-1;
1902: (*split)->ci[1] = -(*split)->nc-1;
1903: }
1904: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1905: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1906: PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1909: /*@
1910: BVGetSplit - Splits the BV object into two BV objects that share the
1911: internal data, one of them containing the leading columns and the other
1912: one containing the remaining columns.
1914: Collective
1916: Input Parameter:
1917: . bv - the basis vectors context
1919: Output Parameters:
1920: + L - left BV containing leading columns (can be NULL)
1921: - R - right BV containing remaining columns (can be NULL)
1923: Notes:
1924: The columns are split in two sets. The leading columns (including the
1925: constraints) are assigned to the left BV and the remaining columns
1926: are assigned to the right BV. The number of leading columns, as
1927: specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1928: guarantee that both L and R have at least one column).
1930: The returned BV's must be seen as references (not copies) of the input
1931: BV, that is, modifying them will change the entries of bv as well.
1932: The returned BV's must not be destroyed. BVRestoreSplit() must be called
1933: when they are no longer needed.
1935: Pass NULL for any of the output BV's that is not needed.
1937: Level: advanced
1939: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints(), BVGetSplitRows()
1940: @*/
1941: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1942: {
1943: PetscFunctionBegin;
1946: BVCheckSizes(bv,1);
1947: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1948: PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1949: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1950: bv->lsplit = bv->nc+bv->l;
1951: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1952: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1953: if (L) *L = bv->L;
1954: if (R) *R = bv->R;
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@
1959: BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1961: Logically Collective
1963: Input Parameters:
1964: + bv - the basis vectors context
1965: . L - left BV obtained with BVGetSplit()
1966: - R - right BV obtained with BVGetSplit()
1968: Note:
1969: The arguments must match the corresponding call to BVGetSplit().
1971: Level: advanced
1973: .seealso: BVGetSplit()
1974: @*/
1975: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1976: {
1977: PetscFunctionBegin;
1980: BVCheckSizes(bv,1);
1981: PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1982: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1983: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
1984: 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()");
1985: 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()");
1987: PetscTryTypeMethod(bv,restoresplit,L,R);
1988: bv->lsplit = 0;
1989: if (L) *L = NULL;
1990: if (R) *R = NULL;
1991: PetscFunctionReturn(PETSC_SUCCESS);
1992: }
1994: /*
1995: Copy all user-provided attributes of V to another BV object W with different layout
1996: */
1997: static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
1998: {
1999: PetscFunctionBegin;
2000: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2001: PetscCall(BVSetVecType(W,V->vtype));
2002: W->ld = V->ld;
2003: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2004: W->orthog_type = V->orthog_type;
2005: W->orthog_ref = V->orthog_ref;
2006: W->orthog_eta = V->orthog_eta;
2007: W->orthog_block = V->orthog_block;
2008: W->vmm = V->vmm;
2009: W->rrandom = V->rrandom;
2010: W->deftol = V->deftol;
2011: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
2012: W->rand = V->rand;
2013: W->sfocalled = V->sfocalled;
2014: PetscTryTypeMethod(V,duplicate,W);
2015: PetscCall(PetscObjectStateIncrease((PetscObject)W));
2016: PetscFunctionReturn(PETSC_SUCCESS);
2017: }
2019: static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2020: {
2021: PetscInt rstart,rend,lstart;
2022: PetscInt N,n;
2023: PetscBool contig;
2025: PetscFunctionBegin;
2026: PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2027: PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2028: PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2029: if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2030: else bv->lsplit = -lstart;
2031: PetscCall(ISGetSize(is,&N));
2032: PetscCall(ISGetLocalSize(is,&n));
2033: if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2034: if (!*split) {
2035: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2036: (*split)->issplit = top? -1: -2;
2037: (*split)->splitparent = bv;
2038: (*split)->N = N;
2039: (*split)->n = n;
2040: (*split)->m = bv->m;
2041: PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2042: }
2043: (*split)->k = bv->k;
2044: (*split)->l = bv->l;
2045: (*split)->nc = bv->nc;
2046: if ((*split)->nc) {
2047: (*split)->ci[0] = -(*split)->nc-1;
2048: (*split)->ci[1] = -(*split)->nc-1;
2049: }
2050: if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2051: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2052: PetscFunctionReturn(PETSC_SUCCESS);
2053: }
2055: /*@
2056: BVGetSplitRows - Splits the BV object into two BV objects that share the
2057: internal data, using a disjoint horizontal splitting.
2059: Collective
2061: Input Parameters:
2062: + bv - the basis vectors context
2063: . isup - the index set that defines the upper part of the horizontal splitting
2064: - islo - the index set that defines the lower part of the horizontal splitting
2066: Output Parameters:
2067: + U - the resulting BV containing the upper rows
2068: - L - the resulting BV containing the lower rows
2070: Notes:
2071: The index sets must be such that every MPI process can extract the selected
2072: rows from its local part of the input BV, and this part must be contiguous.
2073: With one process, isup will list contiguous indices starting from 0, and islo
2074: will contain the remaining indices, hence we refer to upper and lower part.
2075: However, with several processes the indices will be interleaved because
2076: isup will refer to the upper part of the local array.
2078: The intended use of this function is with matrices of MATNEST type, where
2079: MatNestGetISs() will return the appropriate index sets.
2081: The returned BV's must be seen as references (not copies) of the input
2082: BV, that is, modifying them will change the entries of bv as well.
2083: The returned BV's must not be destroyed. BVRestoreSplitRows() must be called
2084: when they are no longer needed.
2086: Pass NULL for any of the output BV's that is not needed.
2088: Level: advanced
2090: .seealso: BVRestoreSplitRows(), BVGetSplit()
2091: @*/
2092: PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2093: {
2094: PetscFunctionBegin;
2097: BVCheckSizes(bv,1);
2100: PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2101: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2102: if (U) {
2103: PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2104: *U = bv->R;
2105: }
2106: if (L) {
2107: PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2108: *L = bv->L;
2109: }
2110: PetscFunctionReturn(PETSC_SUCCESS);
2111: }
2113: /*@
2114: BVRestoreSplitRows - Restore the BV objects obtained with BVGetSplitRows().
2116: Logically Collective
2118: Input Parameters:
2119: + bv - the basis vectors context
2120: . isup - the index set that defines the upper part of the horizontal splitting
2121: . islo - the index set that defines the lower part of the horizontal splitting
2122: . U - the BV containing the upper rows
2123: - L - the BV containing the lower rows
2125: Note:
2126: The arguments must match the corresponding call to BVGetSplitRows().
2128: Level: advanced
2130: .seealso: BVGetSplitRows()
2131: @*/
2132: PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2133: {
2134: PetscFunctionBegin;
2137: BVCheckSizes(bv,1);
2138: PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2139: 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()");
2140: 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()");
2141: PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2142: bv->lsplit = 0;
2143: PetscFunctionReturn(PETSC_SUCCESS);
2144: }
2146: /*@
2147: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2148: definite inner product.
2150: Logically Collective
2152: Input Parameters:
2153: + bv - basis vectors
2154: - deftol - the tolerance
2156: Options Database Key:
2157: . -bv_definite_tol <deftol> - the tolerance
2159: Notes:
2160: When using a non-standard inner product, see BVSetMatrix(), the solver needs
2161: to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
2162: been declared indefinite, the value z'*B*z must be positive, but due to
2163: rounding error a tiny value may become negative. A tolerance is used to
2164: detect this situation. Likewise, in complex arithmetic z'*B*z should be
2165: real, and we use the same tolerance to check whether a nonzero imaginary part
2166: can be considered negligible.
2168: This function sets this tolerance, which defaults to 10*eps, where eps is
2169: the machine epsilon. The default value should be good for most applications.
2171: Level: advanced
2173: .seealso: BVSetMatrix()
2174: @*/
2175: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2176: {
2177: PetscFunctionBegin;
2180: if (deftol == (PetscReal)PETSC_DEFAULT) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2181: else {
2182: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2183: bv->deftol = deftol;
2184: }
2185: PetscFunctionReturn(PETSC_SUCCESS);
2186: }
2188: /*@
2189: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2190: inner product.
2192: Not Collective
2194: Input Parameter:
2195: . bv - the basis vectors
2197: Output Parameter:
2198: . deftol - the tolerance
2200: Level: advanced
2202: .seealso: BVSetDefiniteTolerance()
2203: @*/
2204: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2205: {
2206: PetscFunctionBegin;
2208: PetscAssertPointer(deftol,2);
2209: *deftol = bv->deftol;
2210: PetscFunctionReturn(PETSC_SUCCESS);
2211: }
2213: /*@
2214: BVSetLeadingDimension - Set the leading dimension to be used for storing the BV data.
2216: Not Collective
2218: Input Parameters:
2219: + bv - basis vectors
2220: - ld - the leading dimension
2222: Notes:
2223: This parameter is relevant for BVMAT, though it might be employed in other types
2224: as well.
2226: When the internal data of the BV is stored as a dense matrix, the leading dimension
2227: has the same meaning as in MatDenseSetLDA(), i.e., the distance in number of
2228: elements from one entry of the matrix to the one in the next column at the same
2229: row. The leading dimension refers to the local array, and hence can be different
2230: in different processes.
2232: The user does not need to change this parameter. The default value is equal to the
2233: number of local rows, but this value may be increased a little to guarantee alignment
2234: (especially in the case of GPU storage).
2236: Level: advanced
2238: .seealso: BVGetLeadingDimension()
2239: @*/
2240: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2241: {
2242: PetscFunctionBegin;
2245: 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");
2246: if (ld == PETSC_DEFAULT) bv->ld = 0;
2247: else {
2248: PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2249: bv->ld = ld;
2250: }
2251: PetscFunctionReturn(PETSC_SUCCESS);
2252: }
2254: /*@
2255: BVGetLeadingDimension - Returns the leading dimension of the BV.
2257: Not Collective
2259: Input Parameter:
2260: . bv - the basis vectors
2262: Output Parameter:
2263: . ld - the leading dimension
2265: Level: advanced
2267: Notes:
2268: The returned value may be different in different processes.
2270: The leading dimension must be used when accessing the internal array via
2271: BVGetArray() or BVGetArrayRead().
2273: .seealso: BVSetLeadingDimension(), BVGetArray(), BVGetArrayRead()
2274: @*/
2275: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2276: {
2277: PetscFunctionBegin;
2279: PetscAssertPointer(ld,2);
2280: *ld = bv->ld;
2281: PetscFunctionReturn(PETSC_SUCCESS);
2282: }