Actual source code: bvbasic.c
slepc-3.21.0 2024-03-30
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: /*@C
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: /*@C
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: /*@C
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: /*@C
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: PetscFunctionBegin;
1544: BVCheckSizes(bv,1);
1545: PetscAssertPointer(A,2);
1546: PetscUseTypeMethod(bv,getmat,A);
1547: PetscFunctionReturn(PETSC_SUCCESS);
1548: }
1550: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1551: {
1552: PetscScalar *vv,*aa;
1554: PetscFunctionBegin;
1555: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1556: vv = aa-(bv->nc+bv->l)*bv->ld;
1557: PetscCall(MatDenseResetArray(bv->Aget));
1558: PetscCall(BVRestoreArray(bv,&vv));
1559: *A = NULL;
1560: PetscFunctionReturn(PETSC_SUCCESS);
1561: }
1563: /*@
1564: BVRestoreMat - Restores the Mat obtained with BVGetMat().
1566: Logically Collective
1568: Input Parameters:
1569: + bv - the basis vectors context
1570: - A - the fetched matrix
1572: Note:
1573: A call to this function must match a previous call of BVGetMat().
1574: The effect is that the contents of the Mat are copied back to the
1575: BV internal data structures.
1577: Level: advanced
1579: .seealso: BVGetMat(), BVRestoreArray()
1580: @*/
1581: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1582: {
1583: PetscFunctionBegin;
1585: BVCheckSizes(bv,1);
1586: PetscAssertPointer(A,2);
1587: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1588: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1589: PetscUseTypeMethod(bv,restoremat,A);
1590: PetscFunctionReturn(PETSC_SUCCESS);
1591: }
1593: /*
1594: Copy all user-provided attributes of V to another BV object W
1595: */
1596: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1597: {
1598: PetscFunctionBegin;
1599: PetscCall(PetscLayoutReference(V->map,&W->map));
1600: PetscCall(BVSetVecType(W,V->vtype));
1601: W->ld = V->ld;
1602: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1603: W->orthog_type = V->orthog_type;
1604: W->orthog_ref = V->orthog_ref;
1605: W->orthog_eta = V->orthog_eta;
1606: W->orthog_block = V->orthog_block;
1607: if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1608: W->matrix = V->matrix;
1609: W->indef = V->indef;
1610: W->vmm = V->vmm;
1611: W->rrandom = V->rrandom;
1612: W->deftol = V->deftol;
1613: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1614: W->rand = V->rand;
1615: W->sfocalled = V->sfocalled;
1616: PetscTryTypeMethod(V,duplicate,W);
1617: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1618: PetscFunctionReturn(PETSC_SUCCESS);
1619: }
1621: /*@
1622: BVDuplicate - Creates a new basis vector object of the same type and
1623: dimensions as an existing one.
1625: Collective
1627: Input Parameter:
1628: . V - basis vectors context
1630: Output Parameter:
1631: . W - location to put the new BV
1633: Notes:
1634: The new BV has the same type and dimensions as V. Also, the inner
1635: product matrix and orthogonalization options are copied.
1637: BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1638: for the new basis vectors. Use BVCopy() to copy the contents.
1640: Level: intermediate
1642: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1643: @*/
1644: PetscErrorCode BVDuplicate(BV V,BV *W)
1645: {
1646: PetscFunctionBegin;
1649: BVCheckSizes(V,1);
1650: PetscAssertPointer(W,2);
1651: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1652: (*W)->N = V->N;
1653: (*W)->n = V->n;
1654: (*W)->m = V->m;
1655: (*W)->k = V->m;
1656: PetscCall(BVDuplicate_Private(V,*W));
1657: PetscFunctionReturn(PETSC_SUCCESS);
1658: }
1660: /*@
1661: BVDuplicateResize - Creates a new basis vector object of the same type and
1662: dimensions as an existing one, but with possibly different number of columns.
1664: Collective
1666: Input Parameters:
1667: + V - basis vectors context
1668: - m - the new number of columns
1670: Output Parameter:
1671: . W - location to put the new BV
1673: Note:
1674: This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1675: contents of V are not copied to W.
1677: Level: intermediate
1679: .seealso: BVDuplicate(), BVResize()
1680: @*/
1681: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1682: {
1683: PetscFunctionBegin;
1686: BVCheckSizes(V,1);
1688: PetscAssertPointer(W,3);
1689: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1690: (*W)->N = V->N;
1691: (*W)->n = V->n;
1692: (*W)->m = m;
1693: (*W)->k = m;
1694: PetscCall(BVDuplicate_Private(V,*W));
1695: PetscFunctionReturn(PETSC_SUCCESS);
1696: }
1698: /*@
1699: BVGetCachedBV - Returns a BV object stored internally that holds the
1700: result of B*X after a call to BVApplyMatrixBV().
1702: Collective
1704: Input Parameter:
1705: . bv - the basis vectors context
1707: Output Parameter:
1708: . cached - the cached BV
1710: Note:
1711: The cached BV is created if not available previously.
1713: Level: developer
1715: .seealso: BVApplyMatrixBV()
1716: @*/
1717: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1718: {
1719: PetscFunctionBegin;
1721: PetscAssertPointer(cached,2);
1722: BVCheckSizes(bv,1);
1723: if (!bv->cached) {
1724: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1725: bv->cached->N = bv->N;
1726: bv->cached->n = bv->n;
1727: bv->cached->m = bv->m;
1728: bv->cached->k = bv->m;
1729: PetscCall(BVDuplicate_Private(bv,bv->cached));
1730: }
1731: *cached = bv->cached;
1732: PetscFunctionReturn(PETSC_SUCCESS);
1733: }
1735: /*@
1736: BVCopy - Copies a basis vector object into another one, W <- V.
1738: Logically Collective
1740: Input Parameter:
1741: . V - basis vectors context
1743: Output Parameter:
1744: . W - the copy
1746: Note:
1747: Both V and W must be distributed in the same manner; local copies are
1748: done. Only active columns (excluding the leading ones) are copied.
1749: In the destination W, columns are overwritten starting from the leading ones.
1750: Constraints are not copied.
1752: Level: beginner
1754: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1755: @*/
1756: PetscErrorCode BVCopy(BV V,BV W)
1757: {
1758: PetscScalar *womega;
1759: const PetscScalar *vomega;
1761: PetscFunctionBegin;
1764: BVCheckSizes(V,1);
1765: BVCheckOp(V,1,copy);
1768: BVCheckSizes(W,2);
1769: PetscCheckSameTypeAndComm(V,1,W,2);
1770: 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);
1771: 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);
1772: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1774: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1775: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1776: /* copy signature */
1777: PetscCall(BV_AllocateSignature(W));
1778: PetscCall(VecGetArrayRead(V->omega,&vomega));
1779: PetscCall(VecGetArray(W->omega,&womega));
1780: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1781: PetscCall(VecRestoreArray(W->omega,&womega));
1782: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1783: }
1784: PetscUseTypeMethod(V,copy,W);
1785: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1786: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1787: PetscFunctionReturn(PETSC_SUCCESS);
1788: }
1790: /*@
1791: BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1793: Logically Collective
1795: Input Parameters:
1796: + V - basis vectors context
1797: - j - the column number to be copied
1799: Output Parameter:
1800: . w - the copied column
1802: Note:
1803: Both V and w must be distributed in the same manner; local copies are done.
1805: Level: beginner
1807: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1808: @*/
1809: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1810: {
1811: PetscInt n,N;
1812: Vec z;
1814: PetscFunctionBegin;
1817: BVCheckSizes(V,1);
1820: PetscCheckSameComm(V,1,w,3);
1822: PetscCall(VecGetSize(w,&N));
1823: PetscCall(VecGetLocalSize(w,&n));
1824: 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);
1826: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1827: PetscCall(BVGetColumn(V,j,&z));
1828: PetscCall(VecCopy(z,w));
1829: PetscCall(BVRestoreColumn(V,j,&z));
1830: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1831: PetscFunctionReturn(PETSC_SUCCESS);
1832: }
1834: /*@
1835: BVCopyColumn - Copies the values from one of the columns to another one.
1837: Logically Collective
1839: Input Parameters:
1840: + V - basis vectors context
1841: . j - the number of the source column
1842: - i - the number of the destination column
1844: Level: beginner
1846: .seealso: BVCopy(), BVCopyVec()
1847: @*/
1848: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1849: {
1850: PetscScalar *omega;
1852: PetscFunctionBegin;
1855: BVCheckSizes(V,1);
1858: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1860: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1861: if (V->omega) {
1862: PetscCall(VecGetArray(V->omega,&omega));
1863: omega[i] = omega[j];
1864: PetscCall(VecRestoreArray(V->omega,&omega));
1865: }
1866: PetscUseTypeMethod(V,copycolumn,j,i);
1867: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1868: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1869: PetscFunctionReturn(PETSC_SUCCESS);
1870: }
1872: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1873: {
1874: PetscInt ncols;
1876: PetscFunctionBegin;
1877: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1878: if (*split && ncols!=(*split)->m) PetscCall(BVDestroy(split));
1879: if (!*split) {
1880: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1881: (*split)->issplit = left? 1: 2;
1882: (*split)->splitparent = bv;
1883: (*split)->N = bv->N;
1884: (*split)->n = bv->n;
1885: (*split)->m = bv->m;
1886: (*split)->k = bv->m;
1887: PetscCall(BVDuplicate_Private(bv,*split));
1888: }
1889: (*split)->l = 0;
1890: (*split)->k = left? bv->l: bv->k-bv->l;
1891: (*split)->nc = left? bv->nc: 0;
1892: (*split)->m = ncols-(*split)->nc;
1893: if ((*split)->nc) {
1894: (*split)->ci[0] = -(*split)->nc-1;
1895: (*split)->ci[1] = -(*split)->nc-1;
1896: }
1897: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1898: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1899: PetscFunctionReturn(PETSC_SUCCESS);
1900: }
1902: /*@
1903: BVGetSplit - Splits the BV object into two BV objects that share the
1904: internal data, one of them containing the leading columns and the other
1905: one containing the remaining columns.
1907: Collective
1909: Input Parameter:
1910: . bv - the basis vectors context
1912: Output Parameters:
1913: + L - left BV containing leading columns (can be NULL)
1914: - R - right BV containing remaining columns (can be NULL)
1916: Notes:
1917: The columns are split in two sets. The leading columns (including the
1918: constraints) are assigned to the left BV and the remaining columns
1919: are assigned to the right BV. The number of leading columns, as
1920: specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1921: guarantee that both L and R have at least one column).
1923: The returned BV's must be seen as references (not copies) of the input
1924: BV, that is, modifying them will change the entries of bv as well.
1925: The returned BV's must not be destroyed. BVRestoreSplit() must be called
1926: when they are no longer needed.
1928: Pass NULL for any of the output BV's that is not needed.
1930: Level: advanced
1932: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1933: @*/
1934: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1935: {
1936: PetscFunctionBegin;
1939: BVCheckSizes(bv,1);
1940: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1941: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1942: bv->lsplit = bv->nc+bv->l;
1943: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1944: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1945: if (L) *L = bv->L;
1946: if (R) *R = bv->R;
1947: PetscFunctionReturn(PETSC_SUCCESS);
1948: }
1950: /*@
1951: BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1953: Logically Collective
1955: Input Parameters:
1956: + bv - the basis vectors context
1957: . L - left BV obtained with BVGetSplit()
1958: - R - right BV obtained with BVGetSplit()
1960: Note:
1961: The arguments must match the corresponding call to BVGetSplit().
1963: Level: advanced
1965: .seealso: BVGetSplit()
1966: @*/
1967: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1968: {
1969: PetscFunctionBegin;
1972: BVCheckSizes(bv,1);
1973: PetscCheck(bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1974: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1975: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
1976: 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()");
1977: 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()");
1979: PetscTryTypeMethod(bv,restoresplit,L,R);
1980: bv->lsplit = 0;
1981: if (L) *L = NULL;
1982: if (R) *R = NULL;
1983: PetscFunctionReturn(PETSC_SUCCESS);
1984: }
1986: /*@
1987: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
1988: definite inner product.
1990: Logically Collective
1992: Input Parameters:
1993: + bv - basis vectors
1994: - deftol - the tolerance
1996: Options Database Key:
1997: . -bv_definite_tol <deftol> - the tolerance
1999: Notes:
2000: When using a non-standard inner product, see BVSetMatrix(), the solver needs
2001: to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
2002: been declared indefinite, the value z'*B*z must be positive, but due to
2003: rounding error a tiny value may become negative. A tolerance is used to
2004: detect this situation. Likewise, in complex arithmetic z'*B*z should be
2005: real, and we use the same tolerance to check whether a nonzero imaginary part
2006: can be considered negligible.
2008: This function sets this tolerance, which defaults to 10*eps, where eps is
2009: the machine epsilon. The default value should be good for most applications.
2011: Level: advanced
2013: .seealso: BVSetMatrix()
2014: @*/
2015: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2016: {
2017: PetscFunctionBegin;
2020: if (deftol == (PetscReal)PETSC_DEFAULT) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2021: else {
2022: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2023: bv->deftol = deftol;
2024: }
2025: PetscFunctionReturn(PETSC_SUCCESS);
2026: }
2028: /*@
2029: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2030: inner product.
2032: Not Collective
2034: Input Parameter:
2035: . bv - the basis vectors
2037: Output Parameter:
2038: . deftol - the tolerance
2040: Level: advanced
2042: .seealso: BVSetDefiniteTolerance()
2043: @*/
2044: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2045: {
2046: PetscFunctionBegin;
2048: PetscAssertPointer(deftol,2);
2049: *deftol = bv->deftol;
2050: PetscFunctionReturn(PETSC_SUCCESS);
2051: }
2053: /*@
2054: BVSetLeadingDimension - Set the leading dimension to be used for storing the BV data.
2056: Not Collective
2058: Input Parameters:
2059: + bv - basis vectors
2060: - ld - the leading dimension
2062: Notes:
2063: This parameter is relevant for BVMAT, though it might be employed in other types
2064: as well.
2066: When the internal data of the BV is stored as a dense matrix, the leading dimension
2067: has the same meaning as in MatDenseSetLDA(), i.e., the distance in number of
2068: elements from one entry of the matrix to the one in the next column at the same
2069: row. The leading dimension refers to the local array, and hence can be different
2070: in different processes.
2072: The user does not need to change this parameter. The default value is equal to the
2073: number of local rows, but this value may be increased a little to guarantee alignment
2074: (especially in the case of GPU storage).
2076: Level: advanced
2078: .seealso: BVGetLeadingDimension()
2079: @*/
2080: PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2081: {
2082: PetscFunctionBegin;
2085: 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");
2086: if (ld == PETSC_DEFAULT) bv->ld = 0;
2087: else {
2088: PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2089: bv->ld = ld;
2090: }
2091: PetscFunctionReturn(PETSC_SUCCESS);
2092: }
2094: /*@
2095: BVGetLeadingDimension - Returns the leading dimension of the BV.
2097: Not Collective
2099: Input Parameter:
2100: . bv - the basis vectors
2102: Output Parameter:
2103: . ld - the leading dimension
2105: Level: advanced
2107: Notes:
2108: The returned value may be different in different processes.
2110: The leading dimension must be used when accessing the internal array via
2111: BVGetArray() or BVGetArrayRead().
2113: .seealso: BVSetLeadingDimension(), BVGetArray(), BVGetArrayRead()
2114: @*/
2115: PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2116: {
2117: PetscFunctionBegin;
2119: PetscAssertPointer(ld,2);
2120: *ld = bv->ld;
2121: PetscFunctionReturn(PETSC_SUCCESS);
2122: }