Actual source code: bvbasic.c
slepc-3.19.2 2023-09-05
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;
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;
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;
114: PetscFunctionBegin;
118: 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);
119: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
120: 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);
121: 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);
122: bv->n = n;
123: bv->N = N;
124: bv->m = m;
125: bv->k = m;
126: if (!bv->t) { /* create template vector and get actual dimensions */
127: PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),&bv->t));
128: PetscCall(VecSetSizes(bv->t,bv->n,bv->N));
129: PetscCall(VecSetFromOptions(bv->t));
130: PetscCall(VecGetSize(bv->t,&bv->N));
131: PetscCall(VecGetLocalSize(bv->t,&bv->n));
132: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
133: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
134: 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);
135: }
136: }
137: if (bv->ops->create) {
138: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
139: PetscUseTypeMethod(bv,create);
140: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
141: bv->ops->create = NULL;
142: bv->defersfo = PETSC_FALSE;
143: }
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*@
148: BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
149: Local and global sizes are specified indirectly by passing a template vector.
151: Collective
153: Input Parameters:
154: + bv - the basis vectors
155: . t - the template vector
156: - m - the number of columns
158: Level: beginner
160: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
161: @*/
162: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
163: {
164: PetscInt ma;
166: PetscFunctionBegin;
169: PetscCheckSameComm(bv,1,t,2);
171: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
172: PetscCheck(!bv->t,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Template vector was already set by a previous call to BVSetSizes/FromVec");
173: PetscCall(VecGetSize(t,&bv->N));
174: PetscCall(VecGetLocalSize(t,&bv->n));
175: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
176: PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
177: 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);
178: }
179: bv->m = m;
180: bv->k = m;
181: bv->t = t;
182: PetscCall(PetscObjectReference((PetscObject)t));
183: if (bv->ops->create) {
184: PetscUseTypeMethod(bv,create);
185: bv->ops->create = NULL;
186: bv->defersfo = PETSC_FALSE;
187: }
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*@
192: BVGetSizes - Returns the local and global sizes, and the number of columns.
194: Not Collective
196: Input Parameter:
197: . bv - the basis vectors
199: Output Parameters:
200: + n - the local size
201: . N - the global size
202: - m - the number of columns
204: Note:
205: Normal usage requires that bv has already been given its sizes, otherwise
206: the call fails. However, this function can also be used to determine if
207: a BV object has been initialized completely (sizes and type). For this,
208: call with n=NULL and N=NULL, then a return value of m=0 indicates that
209: the BV object is not ready for use yet.
211: Level: beginner
213: .seealso: BVSetSizes(), BVSetSizesFromVec()
214: @*/
215: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
216: {
217: PetscFunctionBegin;
218: if (!bv) {
219: if (m && !n && !N) *m = 0;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: if (n || N) BVCheckSizes(bv,1);
224: if (n) *n = bv->n;
225: if (N) *N = bv->N;
226: if (m) *m = bv->m;
227: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: BVSetNumConstraints - Set the number of constraints.
234: Logically Collective
236: Input Parameters:
237: + V - basis vectors
238: - nc - number of constraints
240: Notes:
241: This function sets the number of constraints to nc and marks all remaining
242: columns as regular. Normal user would call BVInsertConstraints() instead.
244: If nc is smaller than the previously set value, then some of the constraints
245: are discarded. In particular, using nc=0 removes all constraints preserving
246: the content of regular columns.
248: Level: developer
250: .seealso: BVInsertConstraints()
251: @*/
252: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
253: {
254: PetscInt total,diff,i;
255: Vec x,y;
257: PetscFunctionBegin;
260: PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
262: BVCheckSizes(V,1);
263: PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
265: diff = nc-V->nc;
266: if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
267: total = V->nc+V->m;
268: PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
269: if (diff<0) { /* lessen constraints, shift contents of BV */
270: for (i=0;i<V->m;i++) {
271: PetscCall(BVGetColumn(V,i,&x));
272: PetscCall(BVGetColumn(V,i+diff,&y));
273: PetscCall(VecCopy(x,y));
274: PetscCall(BVRestoreColumn(V,i,&x));
275: PetscCall(BVRestoreColumn(V,i+diff,&y));
276: }
277: }
278: V->nc = nc;
279: V->ci[0] = -V->nc-1;
280: V->ci[1] = -V->nc-1;
281: V->m = total-nc;
282: V->l = PetscMin(V->l,V->m);
283: V->k = PetscMin(V->k,V->m);
284: PetscCall(PetscObjectStateIncrease((PetscObject)V));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: BVGetNumConstraints - Returns the number of constraints.
291: Not Collective
293: Input Parameter:
294: . bv - the basis vectors
296: Output Parameters:
297: . nc - the number of constraints
299: Level: advanced
301: .seealso: BVGetSizes(), BVInsertConstraints()
302: @*/
303: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
304: {
305: PetscFunctionBegin;
308: *nc = bv->nc;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@
313: BVResize - Change the number of columns.
315: Collective
317: Input Parameters:
318: + bv - the basis vectors
319: . m - the new number of columns
320: - copy - a flag indicating whether current values should be kept
322: Note:
323: Internal storage is reallocated. If the copy flag is set to true, then
324: the contents are copied to the leading part of the new space.
326: Level: advanced
328: .seealso: BVSetSizes(), BVSetSizesFromVec()
329: @*/
330: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
331: {
332: PetscScalar *array;
333: const PetscScalar *omega;
334: Vec v;
336: PetscFunctionBegin;
341: PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
342: PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
343: if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
344: BVCheckOp(bv,1,resize);
346: PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
347: PetscUseTypeMethod(bv,resize,m,copy);
348: PetscCall(VecDestroy(&bv->buffer));
349: PetscCall(BVDestroy(&bv->cached));
350: PetscCall(PetscFree2(bv->h,bv->c));
351: if (bv->omega) {
352: if (bv->cuda) {
353: #if defined(PETSC_HAVE_CUDA)
354: PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
355: #else
356: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
357: #endif
358: } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
359: if (copy) {
360: PetscCall(VecGetArray(v,&array));
361: PetscCall(VecGetArrayRead(bv->omega,&omega));
362: PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
363: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
364: PetscCall(VecRestoreArray(v,&array));
365: } else PetscCall(VecSet(v,1.0));
366: PetscCall(VecDestroy(&bv->omega));
367: bv->omega = v;
368: }
369: bv->m = m;
370: bv->k = PetscMin(bv->k,m);
371: bv->l = PetscMin(bv->l,m);
372: PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
373: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: /*@
378: BVSetActiveColumns - Specify the columns that will be involved in operations.
380: Logically Collective
382: Input Parameters:
383: + bv - the basis vectors context
384: . l - number of leading columns
385: - k - number of active columns
387: Notes:
388: In operations such as BVMult() or BVDot(), only the first k columns are
389: considered. This is useful when the BV is filled from left to right, so
390: the last m-k columns do not have relevant information.
392: Also in operations such as BVMult() or BVDot(), the first l columns are
393: normally not included in the computation. See the manpage of each
394: operation.
396: In orthogonalization operations, the first l columns are treated
397: differently, they participate in the orthogonalization but the computed
398: coefficients are not stored.
400: Level: intermediate
402: .seealso: BVGetActiveColumns(), BVSetSizes()
403: @*/
404: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
405: {
406: PetscFunctionBegin;
410: BVCheckSizes(bv,1);
411: if (PetscUnlikely(k==PETSC_DECIDE || k==PETSC_DEFAULT)) {
412: bv->k = bv->m;
413: } else {
414: 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);
415: bv->k = k;
416: }
417: if (PetscUnlikely(l==PETSC_DECIDE || l==PETSC_DEFAULT)) {
418: bv->l = 0;
419: } else {
420: 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);
421: bv->l = l;
422: }
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*@
427: BVGetActiveColumns - Returns the current active dimensions.
429: Not Collective
431: Input Parameter:
432: . bv - the basis vectors context
434: Output Parameters:
435: + l - number of leading columns
436: - k - number of active columns
438: Level: intermediate
440: .seealso: BVSetActiveColumns()
441: @*/
442: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
443: {
444: PetscFunctionBegin;
446: if (l) *l = bv->l;
447: if (k) *k = bv->k;
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*@
452: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
454: Collective
456: Input Parameters:
457: + bv - the basis vectors context
458: . B - a symmetric matrix (may be NULL)
459: - indef - a flag indicating if the matrix is indefinite
461: Notes:
462: This is used to specify a non-standard inner product, whose matrix
463: representation is given by B. Then, all inner products required during
464: orthogonalization are computed as (x,y)_B=y^H*B*x rather than the
465: standard form (x,y)=y^H*x.
467: Matrix B must be real symmetric (or complex Hermitian). A genuine inner
468: product requires that B is also positive (semi-)definite. However, we
469: also allow for an indefinite B (setting indef=PETSC_TRUE), in which
470: case the orthogonalization uses an indefinite inner product.
472: This affects operations BVDot(), BVNorm(), BVOrthogonalize(), and variants.
474: Setting B=NULL has the same effect as if the identity matrix was passed.
476: Level: advanced
478: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
479: @*/
480: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
481: {
482: PetscInt m,n;
484: PetscFunctionBegin;
487: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
488: if (B) {
490: PetscCall(MatGetLocalSize(B,&m,&n));
491: PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
492: 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);
493: }
494: if (B) PetscCall(PetscObjectReference((PetscObject)B));
495: PetscCall(MatDestroy(&bv->matrix));
496: bv->matrix = B;
497: bv->indef = indef;
498: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
499: if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
500: if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
501: }
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*@
506: BVGetMatrix - Retrieves the matrix representation of the inner product.
508: Not Collective
510: Input Parameter:
511: . bv - the basis vectors context
513: Output Parameters:
514: + B - the matrix of the inner product (may be NULL)
515: - indef - the flag indicating if the matrix is indefinite
517: Level: advanced
519: .seealso: BVSetMatrix()
520: @*/
521: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
522: {
523: PetscFunctionBegin;
525: if (B) *B = bv->matrix;
526: if (indef) *indef = bv->indef;
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*@
531: BVApplyMatrix - Multiplies a vector by the matrix representation of the
532: inner product.
534: Neighbor-wise Collective
536: Input Parameters:
537: + bv - the basis vectors context
538: - x - the vector
540: Output Parameter:
541: . y - the result
543: Note:
544: If no matrix was specified this function copies the vector.
546: Level: advanced
548: .seealso: BVSetMatrix(), BVApplyMatrixBV()
549: @*/
550: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
551: {
552: PetscFunctionBegin;
556: if (bv->matrix) {
557: PetscCall(BV_IPMatMult(bv,x));
558: PetscCall(VecCopy(bv->Bx,y));
559: } else PetscCall(VecCopy(x,y));
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
565: of the inner product.
567: Neighbor-wise Collective
569: Input Parameter:
570: . X - the basis vectors context
572: Output Parameter:
573: . Y - the basis vectors to store the result (optional)
575: Note:
576: This function computes Y = B*X, where B is the matrix given with
577: BVSetMatrix(). This operation is computed as in BVMatMult().
578: If no matrix was specified, then it just copies Y = X.
580: If no Y is given, the result is stored internally in the cached BV.
582: Level: developer
584: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
585: @*/
586: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
587: {
588: PetscFunctionBegin;
590: if (Y) {
592: if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
593: else PetscCall(BVCopy(X,Y));
594: } else PetscCall(BV_IPMatMultBV(X));
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: /*@
599: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
601: Logically Collective
603: Input Parameters:
604: + bv - the basis vectors context
605: - omega - a vector representing the diagonal of the signature matrix
607: Note:
608: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
610: Level: developer
612: .seealso: BVSetMatrix(), BVGetSignature()
613: @*/
614: PetscErrorCode BVSetSignature(BV bv,Vec omega)
615: {
616: PetscInt i,n;
617: const PetscScalar *pomega;
618: PetscScalar *intern;
620: PetscFunctionBegin;
622: BVCheckSizes(bv,1);
626: PetscCall(VecGetSize(omega,&n));
627: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
628: PetscCall(BV_AllocateSignature(bv));
629: if (bv->indef) {
630: PetscCall(VecGetArrayRead(omega,&pomega));
631: PetscCall(VecGetArray(bv->omega,&intern));
632: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
633: PetscCall(VecRestoreArray(bv->omega,&intern));
634: PetscCall(VecRestoreArrayRead(omega,&pomega));
635: } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
636: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: /*@
641: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
643: Not Collective
645: Input Parameter:
646: . bv - the basis vectors context
648: Output Parameter:
649: . omega - a vector representing the diagonal of the signature matrix
651: Note:
652: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
654: Level: developer
656: .seealso: BVSetMatrix(), BVSetSignature()
657: @*/
658: PetscErrorCode BVGetSignature(BV bv,Vec omega)
659: {
660: PetscInt i,n;
661: PetscScalar *pomega;
662: const PetscScalar *intern;
664: PetscFunctionBegin;
666: BVCheckSizes(bv,1);
670: PetscCall(VecGetSize(omega,&n));
671: PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
672: if (bv->indef && bv->omega) {
673: PetscCall(VecGetArray(omega,&pomega));
674: PetscCall(VecGetArrayRead(bv->omega,&intern));
675: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
676: PetscCall(VecRestoreArrayRead(bv->omega,&intern));
677: PetscCall(VecRestoreArray(omega,&pomega));
678: } else PetscCall(VecSet(omega,1.0));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: /*@
683: BVSetBufferVec - Attach a vector object to be used as buffer space for
684: several operations.
686: Collective
688: Input Parameters:
689: + bv - the basis vectors context)
690: - buffer - the vector
692: Notes:
693: Use BVGetBufferVec() to retrieve the vector (for example, to free it
694: at the end of the computations).
696: The vector must be sequential of length (nc+m)*m, where m is the number
697: of columns of bv and nc is the number of constraints.
699: Level: developer
701: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
702: @*/
703: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
704: {
705: PetscInt ld,n;
706: PetscMPIInt size;
708: PetscFunctionBegin;
711: BVCheckSizes(bv,1);
712: PetscCall(VecGetSize(buffer,&n));
713: ld = bv->m+bv->nc;
714: PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
715: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
716: PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
718: PetscCall(PetscObjectReference((PetscObject)buffer));
719: PetscCall(VecDestroy(&bv->buffer));
720: bv->buffer = buffer;
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: BVGetBufferVec - Obtain the buffer vector associated with the BV object.
727: Collective
729: Input Parameters:
730: . bv - the basis vectors context
732: Output Parameter:
733: . buffer - vector
735: Notes:
736: The vector is created if not available previously. It is a sequential vector
737: of length (nc+m)*m, where m is the number of columns of bv and nc is the number
738: of constraints.
740: Developer Notes:
741: The buffer vector is viewed as a column-major matrix with leading dimension
742: ld=nc+m, and m columns at most. In the most common usage, it has the structure
743: .vb
744: | | C |
745: |s|---|
746: | | H |
747: .ve
748: where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
749: related to orthogonalization against constraints (first nc rows), and s is the
750: first column that contains scratch values computed during Gram-Schmidt
751: orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
752: store the coefficients.
754: Level: developer
756: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
757: @*/
758: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
759: {
760: PetscInt ld;
762: PetscFunctionBegin;
765: BVCheckSizes(bv,1);
766: if (!bv->buffer) {
767: ld = bv->m+bv->nc;
768: PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
769: PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
770: PetscCall(VecSetType(bv->buffer,((PetscObject)bv->t)->type_name));
771: }
772: *buffer = bv->buffer;
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: /*@
777: BVSetRandomContext - Sets the PetscRandom object associated with the BV,
778: to be used in operations that need random numbers.
780: Collective
782: Input Parameters:
783: + bv - the basis vectors context
784: - rand - the random number generator context
786: Level: advanced
788: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
789: @*/
790: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
791: {
792: PetscFunctionBegin;
795: PetscCheckSameComm(bv,1,rand,2);
796: PetscCall(PetscObjectReference((PetscObject)rand));
797: PetscCall(PetscRandomDestroy(&bv->rand));
798: bv->rand = rand;
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: /*@
803: BVGetRandomContext - Gets the PetscRandom object associated with the BV.
805: Collective
807: Input Parameter:
808: . bv - the basis vectors context
810: Output Parameter:
811: . rand - the random number generator context
813: Level: advanced
815: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
816: @*/
817: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
818: {
819: PetscFunctionBegin;
822: if (!bv->rand) {
823: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
824: if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
825: if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
826: if (bv->rrandom) {
827: PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
828: PetscCall(PetscRandomSeed(bv->rand));
829: }
830: }
831: *rand = bv->rand;
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /*@
836: BVSetFromOptions - Sets BV options from the options database.
838: Collective
840: Input Parameter:
841: . bv - the basis vectors context
843: Level: beginner
845: .seealso: BVSetOptionsPrefix()
846: @*/
847: PetscErrorCode BVSetFromOptions(BV bv)
848: {
849: char type[256];
850: PetscBool flg1,flg2,flg3,flg4;
851: PetscReal r;
852: BVOrthogType otype;
853: BVOrthogRefineType orefine;
854: BVOrthogBlockType oblock;
856: PetscFunctionBegin;
858: PetscCall(BVRegisterAll());
859: PetscObjectOptionsBegin((PetscObject)bv);
860: PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVSVEC),type,sizeof(type),&flg1));
861: if (flg1) PetscCall(BVSetType(bv,type));
862: else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVSVEC));
864: otype = bv->orthog_type;
865: PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
866: orefine = bv->orthog_ref;
867: PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
868: r = bv->orthog_eta;
869: PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
870: oblock = bv->orthog_block;
871: PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
872: if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
874: PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
876: PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
877: if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
879: /* undocumented option to generate random vectors that are independent of the number of processes */
880: PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
882: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
883: else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
884: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
885: PetscOptionsEnd();
886: bv->sfocalled = PETSC_TRUE;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: /*@
891: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
892: vectors (classical or modified Gram-Schmidt with or without refinement), and
893: for the block-orthogonalization (simultaneous orthogonalization of a set of
894: vectors).
896: Logically Collective
898: Input Parameters:
899: + bv - the basis vectors context
900: . type - the method of vector orthogonalization
901: . refine - type of refinement
902: . eta - parameter for selective refinement
903: - block - the method of block orthogonalization
905: Options Database Keys:
906: + -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
907: (default) or mgs for Modified Gram-Schmidt orthogonalization
908: . -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
909: . -bv_orthog_eta <eta> - For setting the value of eta
910: - -bv_orthog_block <block> - Where <block> is the block-orthogonalization method
912: Notes:
913: The default settings work well for most problems.
915: The parameter eta should be a real value between 0 and 1 (or PETSC_DEFAULT).
916: The value of eta is used only when the refinement type is "ifneeded".
918: When using several processors, MGS is likely to result in bad scalability.
920: If the method set for block orthogonalization is GS, then the computation
921: is done column by column with the vector orthogonalization.
923: Level: advanced
925: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
926: @*/
927: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
928: {
929: PetscFunctionBegin;
935: switch (type) {
936: case BV_ORTHOG_CGS:
937: case BV_ORTHOG_MGS:
938: bv->orthog_type = type;
939: break;
940: default:
941: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
942: }
943: switch (refine) {
944: case BV_ORTHOG_REFINE_NEVER:
945: case BV_ORTHOG_REFINE_IFNEEDED:
946: case BV_ORTHOG_REFINE_ALWAYS:
947: bv->orthog_ref = refine;
948: break;
949: default:
950: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
951: }
952: if (eta == (PetscReal)PETSC_DEFAULT) {
953: bv->orthog_eta = 0.7071;
954: } else {
955: PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
956: bv->orthog_eta = eta;
957: }
958: switch (block) {
959: case BV_ORTHOG_BLOCK_GS:
960: case BV_ORTHOG_BLOCK_CHOL:
961: case BV_ORTHOG_BLOCK_TSQR:
962: case BV_ORTHOG_BLOCK_TSQRCHOL:
963: case BV_ORTHOG_BLOCK_SVQB:
964: bv->orthog_block = block;
965: break;
966: default:
967: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
968: }
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: /*@
973: BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.
975: Not Collective
977: Input Parameter:
978: . bv - basis vectors context
980: Output Parameters:
981: + type - the method of vector orthogonalization
982: . refine - type of refinement
983: . eta - parameter for selective refinement
984: - block - the method of block orthogonalization
986: Level: advanced
988: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
989: @*/
990: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
991: {
992: PetscFunctionBegin;
994: if (type) *type = bv->orthog_type;
995: if (refine) *refine = bv->orthog_ref;
996: if (eta) *eta = bv->orthog_eta;
997: if (block) *block = bv->orthog_block;
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: /*@
1002: BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.
1004: Logically Collective
1006: Input Parameters:
1007: + bv - the basis vectors context
1008: - method - the method for the BVMatMult() operation
1010: Options Database Keys:
1011: . -bv_matmult <meth> - choose one of the methods: vecs, mat
1013: Notes:
1014: Allowed values are
1015: + BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1016: . BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1017: - BV_MATMULT_MAT_SAVE - this case is deprecated
1019: The default is BV_MATMULT_MAT except in the case of BVVECS.
1021: Level: advanced
1023: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1024: @*/
1025: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1026: {
1027: PetscFunctionBegin;
1030: switch (method) {
1031: case BV_MATMULT_VECS:
1032: case BV_MATMULT_MAT:
1033: bv->vmm = method;
1034: break;
1035: case BV_MATMULT_MAT_SAVE:
1036: PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1037: bv->vmm = BV_MATMULT_MAT;
1038: break;
1039: default:
1040: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1041: }
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /*@
1046: BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.
1048: Not Collective
1050: Input Parameter:
1051: . bv - basis vectors context
1053: Output Parameter:
1054: . method - the method for the BVMatMult() operation
1056: Level: advanced
1058: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1059: @*/
1060: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1061: {
1062: PetscFunctionBegin;
1065: *method = bv->vmm;
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@
1070: BVGetColumn - Returns a Vec object that contains the entries of the
1071: requested column of the basis vectors object.
1073: Logically Collective
1075: Input Parameters:
1076: + bv - the basis vectors context
1077: - j - the index of the requested column
1079: Output Parameter:
1080: . v - vector containing the jth column
1082: Notes:
1083: The returned Vec must be seen as a reference (not a copy) of the BV
1084: column, that is, modifying the Vec will change the BV entries as well.
1086: The returned Vec must not be destroyed. BVRestoreColumn() must be
1087: called when it is no longer needed. At most, two columns can be fetched,
1088: that is, this function can only be called twice before the corresponding
1089: BVRestoreColumn() is invoked.
1091: A negative index j selects the i-th constraint, where i=-j. Constraints
1092: should not be modified.
1094: Level: beginner
1096: .seealso: BVRestoreColumn(), BVInsertConstraints()
1097: @*/
1098: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1099: {
1100: PetscInt l;
1102: PetscFunctionBegin;
1105: BVCheckSizes(bv,1);
1106: BVCheckOp(bv,1,getcolumn);
1108: 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);
1109: 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);
1110: 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);
1111: l = BVAvailableVec;
1112: PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1113: PetscUseTypeMethod(bv,getcolumn,j,v);
1114: bv->ci[l] = j;
1115: PetscCall(PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]));
1116: PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1117: *v = bv->cv[l];
1118: PetscFunctionReturn(PETSC_SUCCESS);
1119: }
1121: /*@
1122: BVRestoreColumn - Restore a column obtained with BVGetColumn().
1124: Logically Collective
1126: Input Parameters:
1127: + bv - the basis vectors context
1128: . j - the index of the column
1129: - v - vector obtained with BVGetColumn()
1131: Note:
1132: The arguments must match the corresponding call to BVGetColumn().
1134: Level: beginner
1136: .seealso: BVGetColumn()
1137: @*/
1138: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1139: {
1140: PetscObjectId id;
1141: PetscObjectState st;
1142: PetscInt l;
1144: PetscFunctionBegin;
1147: BVCheckSizes(bv,1);
1151: 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);
1152: 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);
1153: 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);
1154: l = (j==bv->ci[0])? 0: 1;
1155: PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1156: PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1157: PetscCall(PetscObjectStateGet((PetscObject)*v,&st));
1158: if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1159: PetscUseTypeMethod(bv,restorecolumn,j,v);
1160: bv->ci[l] = -bv->nc-1;
1161: bv->st[l] = -1;
1162: bv->id[l] = 0;
1163: *v = NULL;
1164: PetscFunctionReturn(PETSC_SUCCESS);
1165: }
1167: /*@C
1168: BVGetArray - Returns a pointer to a contiguous array that contains this
1169: processor's portion of the BV data.
1171: Logically Collective
1173: Input Parameters:
1174: . bv - the basis vectors context
1176: Output Parameter:
1177: . a - location to put pointer to the array
1179: Notes:
1180: BVRestoreArray() must be called when access to the array is no longer needed.
1181: This operation may imply a data copy, for BV types that do not store
1182: data contiguously in memory.
1184: The pointer will normally point to the first entry of the first column,
1185: but if the BV has constraints then these go before the regular columns.
1187: Level: advanced
1189: .seealso: BVRestoreArray(), BVInsertConstraints()
1190: @*/
1191: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1192: {
1193: PetscFunctionBegin;
1196: BVCheckSizes(bv,1);
1197: BVCheckOp(bv,1,getarray);
1198: PetscUseTypeMethod(bv,getarray,a);
1199: PetscFunctionReturn(PETSC_SUCCESS);
1200: }
1202: /*@C
1203: BVRestoreArray - Restore the BV object after BVGetArray() has been called.
1205: Logically Collective
1207: Input Parameters:
1208: + bv - the basis vectors context
1209: - a - location of pointer to array obtained from BVGetArray()
1211: Note:
1212: This operation may imply a data copy, for BV types that do not store
1213: data contiguously in memory.
1215: Level: advanced
1217: .seealso: BVGetColumn()
1218: @*/
1219: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1220: {
1221: PetscFunctionBegin;
1224: BVCheckSizes(bv,1);
1225: PetscTryTypeMethod(bv,restorearray,a);
1226: if (a) *a = NULL;
1227: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: /*@C
1232: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1233: contains this processor's portion of the BV data.
1235: Not Collective
1237: Input Parameters:
1238: . bv - the basis vectors context
1240: Output Parameter:
1241: . a - location to put pointer to the array
1243: Notes:
1244: BVRestoreArrayRead() must be called when access to the array is no
1245: longer needed. This operation may imply a data copy, for BV types that
1246: do not store data contiguously in memory.
1248: The pointer will normally point to the first entry of the first column,
1249: but if the BV has constraints then these go before the regular columns.
1251: Level: advanced
1253: .seealso: BVRestoreArray(), BVInsertConstraints()
1254: @*/
1255: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1256: {
1257: PetscFunctionBegin;
1260: BVCheckSizes(bv,1);
1261: BVCheckOp(bv,1,getarrayread);
1262: PetscUseTypeMethod(bv,getarrayread,a);
1263: PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1266: /*@C
1267: BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1268: been called.
1270: Not Collective
1272: Input Parameters:
1273: + bv - the basis vectors context
1274: - a - location of pointer to array obtained from BVGetArrayRead()
1276: Level: advanced
1278: .seealso: BVGetColumn()
1279: @*/
1280: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1281: {
1282: PetscFunctionBegin;
1285: BVCheckSizes(bv,1);
1286: PetscTryTypeMethod(bv,restorearrayread,a);
1287: if (a) *a = NULL;
1288: PetscFunctionReturn(PETSC_SUCCESS);
1289: }
1291: /*@
1292: BVCreateVec - Creates a new Vec object with the same type and dimensions
1293: as the columns of the basis vectors object.
1295: Collective
1297: Input Parameter:
1298: . bv - the basis vectors context
1300: Output Parameter:
1301: . v - the new vector
1303: Note:
1304: The user is responsible of destroying the returned vector.
1306: Level: beginner
1308: .seealso: BVCreateMat()
1309: @*/
1310: PetscErrorCode BVCreateVec(BV bv,Vec *v)
1311: {
1312: PetscFunctionBegin;
1314: BVCheckSizes(bv,1);
1316: PetscCall(VecDuplicate(bv->t,v));
1317: PetscFunctionReturn(PETSC_SUCCESS);
1318: }
1320: /*@
1321: BVCreateMat - Creates a new Mat object of dense type and copies the contents
1322: of the BV object.
1324: Collective
1326: Input Parameter:
1327: . bv - the basis vectors context
1329: Output Parameter:
1330: . A - the new matrix
1332: Notes:
1333: The user is responsible of destroying the returned matrix.
1335: The matrix contains all columns of the BV, not just the active columns.
1337: Level: intermediate
1339: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1340: @*/
1341: PetscErrorCode BVCreateMat(BV bv,Mat *A)
1342: {
1343: PetscScalar *aa;
1344: const PetscScalar *vv;
1346: PetscFunctionBegin;
1348: BVCheckSizes(bv,1);
1351: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,bv->N,bv->m,NULL,A));
1352: PetscCall(MatDenseGetArrayWrite(*A,&aa));
1353: PetscCall(BVGetArrayRead(bv,&vv));
1354: PetscCall(PetscArraycpy(aa,vv,bv->m*bv->n));
1355: PetscCall(BVRestoreArrayRead(bv,&vv));
1356: PetscCall(MatDenseRestoreArrayWrite(*A,&aa));
1357: PetscFunctionReturn(PETSC_SUCCESS);
1358: }
1360: PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1361: {
1362: PetscScalar *vv,*aa;
1363: PetscBool create=PETSC_FALSE;
1364: PetscInt m,cols;
1366: PetscFunctionBegin;
1367: m = bv->k-bv->l;
1368: if (!bv->Aget) create=PETSC_TRUE;
1369: else {
1370: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1371: PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1372: PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1373: if (cols!=m) {
1374: PetscCall(MatDestroy(&bv->Aget));
1375: create=PETSC_TRUE;
1376: }
1377: }
1378: PetscCall(BVGetArray(bv,&vv));
1379: if (create) {
1380: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,vv,&bv->Aget)); /* pass a pointer to avoid allocation of storage */
1381: PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1382: }
1383: PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->n)); /* set the actual pointer */
1384: *A = bv->Aget;
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: /*@
1389: BVGetMat - Returns a Mat object of dense type that shares the memory of
1390: the BV object.
1392: Collective
1394: Input Parameter:
1395: . bv - the basis vectors context
1397: Output Parameter:
1398: . A - the matrix
1400: Notes:
1401: The returned matrix contains only the active columns. If the content of
1402: the Mat is modified, these changes are also done in the BV object. The
1403: user must call BVRestoreMat() when no longer needed.
1405: This operation implies a call to BVGetArray(), which may result in data
1406: copies.
1408: Level: advanced
1410: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1411: @*/
1412: PetscErrorCode BVGetMat(BV bv,Mat *A)
1413: {
1414: PetscFunctionBegin;
1416: BVCheckSizes(bv,1);
1418: PetscUseTypeMethod(bv,getmat,A);
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1423: {
1424: PetscScalar *vv,*aa;
1426: PetscFunctionBegin;
1427: PetscCall(MatDenseGetArray(bv->Aget,&aa));
1428: vv = aa-(bv->nc+bv->l)*bv->n;
1429: PetscCall(MatDenseResetArray(bv->Aget));
1430: PetscCall(BVRestoreArray(bv,&vv));
1431: *A = NULL;
1432: PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1435: /*@
1436: BVRestoreMat - Restores the Mat obtained with BVGetMat().
1438: Logically Collective
1440: Input Parameters:
1441: + bv - the basis vectors context
1442: - A - the fetched matrix
1444: Note:
1445: A call to this function must match a previous call of BVGetMat().
1446: The effect is that the contents of the Mat are copied back to the
1447: BV internal data structures.
1449: Level: advanced
1451: .seealso: BVGetMat(), BVRestoreArray()
1452: @*/
1453: PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1454: {
1455: PetscFunctionBegin;
1457: BVCheckSizes(bv,1);
1459: PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1460: PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1461: PetscUseTypeMethod(bv,restoremat,A);
1462: PetscFunctionReturn(PETSC_SUCCESS);
1463: }
1465: /*
1466: Copy all user-provided attributes of V to another BV object W
1467: */
1468: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1469: {
1470: PetscFunctionBegin;
1471: PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1472: W->orthog_type = V->orthog_type;
1473: W->orthog_ref = V->orthog_ref;
1474: W->orthog_eta = V->orthog_eta;
1475: W->orthog_block = V->orthog_block;
1476: if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1477: W->matrix = V->matrix;
1478: W->indef = V->indef;
1479: W->vmm = V->vmm;
1480: W->rrandom = V->rrandom;
1481: W->deftol = V->deftol;
1482: if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1483: W->rand = V->rand;
1484: W->sfocalled = V->sfocalled;
1485: PetscTryTypeMethod(V,duplicate,W);
1486: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1487: PetscFunctionReturn(PETSC_SUCCESS);
1488: }
1490: /*@
1491: BVDuplicate - Creates a new basis vector object of the same type and
1492: dimensions as an existing one.
1494: Collective
1496: Input Parameter:
1497: . V - basis vectors context
1499: Output Parameter:
1500: . W - location to put the new BV
1502: Notes:
1503: The new BV has the same type and dimensions as V, and it shares the same
1504: template vector. Also, the inner product matrix and orthogonalization
1505: options are copied.
1507: BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1508: for the new basis vectors. Use BVCopy() to copy the contents.
1510: Level: intermediate
1512: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1513: @*/
1514: PetscErrorCode BVDuplicate(BV V,BV *W)
1515: {
1516: PetscFunctionBegin;
1519: BVCheckSizes(V,1);
1521: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1522: PetscCall(BVSetSizesFromVec(*W,V->t,V->m));
1523: PetscCall(BVDuplicate_Private(V,*W));
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: /*@
1528: BVDuplicateResize - Creates a new basis vector object of the same type and
1529: dimensions as an existing one, but with possibly different number of columns.
1531: Collective
1533: Input Parameters:
1534: + V - basis vectors context
1535: - m - the new number of columns
1537: Output Parameter:
1538: . W - location to put the new BV
1540: Note:
1541: This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1542: contents of V are not copied to W.
1544: Level: intermediate
1546: .seealso: BVDuplicate(), BVResize()
1547: @*/
1548: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1549: {
1550: PetscFunctionBegin;
1553: BVCheckSizes(V,1);
1556: PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1557: PetscCall(BVSetSizesFromVec(*W,V->t,m));
1558: PetscCall(BVDuplicate_Private(V,*W));
1559: PetscFunctionReturn(PETSC_SUCCESS);
1560: }
1562: /*@
1563: BVGetCachedBV - Returns a BV object stored internally that holds the
1564: result of B*X after a call to BVApplyMatrixBV().
1566: Collective
1568: Input Parameter:
1569: . bv - the basis vectors context
1571: Output Parameter:
1572: . cached - the cached BV
1574: Note:
1575: The cached BV is created if not available previously.
1577: Level: developer
1579: .seealso: BVApplyMatrixBV()
1580: @*/
1581: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1582: {
1583: PetscFunctionBegin;
1586: BVCheckSizes(bv,1);
1587: if (!bv->cached) {
1588: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1589: PetscCall(BVSetSizesFromVec(bv->cached,bv->t,bv->m));
1590: PetscCall(BVDuplicate_Private(bv,bv->cached));
1591: }
1592: *cached = bv->cached;
1593: PetscFunctionReturn(PETSC_SUCCESS);
1594: }
1596: /*@
1597: BVCopy - Copies a basis vector object into another one, W <- V.
1599: Logically Collective
1601: Input Parameter:
1602: . V - basis vectors context
1604: Output Parameter:
1605: . W - the copy
1607: Note:
1608: Both V and W must be distributed in the same manner; local copies are
1609: done. Only active columns (excluding the leading ones) are copied.
1610: In the destination W, columns are overwritten starting from the leading ones.
1611: Constraints are not copied.
1613: Level: beginner
1615: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1616: @*/
1617: PetscErrorCode BVCopy(BV V,BV W)
1618: {
1619: PetscScalar *womega;
1620: const PetscScalar *vomega;
1622: PetscFunctionBegin;
1625: BVCheckSizes(V,1);
1626: BVCheckOp(V,1,copy);
1629: BVCheckSizes(W,2);
1630: PetscCheckSameTypeAndComm(V,1,W,2);
1631: 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);
1632: 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);
1633: if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1635: PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1636: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1637: /* copy signature */
1638: PetscCall(BV_AllocateSignature(W));
1639: PetscCall(VecGetArrayRead(V->omega,&vomega));
1640: PetscCall(VecGetArray(W->omega,&womega));
1641: PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1642: PetscCall(VecRestoreArray(W->omega,&womega));
1643: PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1644: }
1645: PetscUseTypeMethod(V,copy,W);
1646: PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1647: PetscCall(PetscObjectStateIncrease((PetscObject)W));
1648: PetscFunctionReturn(PETSC_SUCCESS);
1649: }
1651: /*@
1652: BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1654: Logically Collective
1656: Input Parameters:
1657: + V - basis vectors context
1658: - j - the column number to be copied
1660: Output Parameter:
1661: . w - the copied column
1663: Note:
1664: Both V and w must be distributed in the same manner; local copies are done.
1666: Level: beginner
1668: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1669: @*/
1670: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1671: {
1672: PetscInt n,N;
1673: Vec z;
1675: PetscFunctionBegin;
1678: BVCheckSizes(V,1);
1681: PetscCheckSameComm(V,1,w,3);
1683: PetscCall(VecGetSize(w,&N));
1684: PetscCall(VecGetLocalSize(w,&n));
1685: 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);
1687: PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1688: PetscCall(BVGetColumn(V,j,&z));
1689: PetscCall(VecCopy(z,w));
1690: PetscCall(BVRestoreColumn(V,j,&z));
1691: PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1692: PetscFunctionReturn(PETSC_SUCCESS);
1693: }
1695: /*@
1696: BVCopyColumn - Copies the values from one of the columns to another one.
1698: Logically Collective
1700: Input Parameters:
1701: + V - basis vectors context
1702: . j - the number of the source column
1703: - i - the number of the destination column
1705: Level: beginner
1707: .seealso: BVCopy(), BVCopyVec()
1708: @*/
1709: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1710: {
1711: PetscScalar *omega;
1713: PetscFunctionBegin;
1716: BVCheckSizes(V,1);
1719: if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1721: PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1722: if (V->omega) {
1723: PetscCall(VecGetArray(V->omega,&omega));
1724: omega[i] = omega[j];
1725: PetscCall(VecRestoreArray(V->omega,&omega));
1726: }
1727: PetscUseTypeMethod(V,copycolumn,j,i);
1728: PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1729: PetscCall(PetscObjectStateIncrease((PetscObject)V));
1730: PetscFunctionReturn(PETSC_SUCCESS);
1731: }
1733: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1734: {
1735: PetscInt ncols;
1737: PetscFunctionBegin;
1738: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1739: if (*split && ncols!=(*split)->m) PetscCall(BVDestroy(split));
1740: if (!*split) {
1741: PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1742: (*split)->issplit = left? 1: 2;
1743: (*split)->splitparent = bv;
1744: PetscCall(BVSetSizesFromVec(*split,bv->t,ncols));
1745: PetscCall(BVDuplicate_Private(bv,*split));
1746: }
1747: (*split)->l = 0;
1748: (*split)->k = left? bv->l: bv->k-bv->l;
1749: (*split)->nc = left? bv->nc: 0;
1750: (*split)->m = ncols-(*split)->nc;
1751: if ((*split)->nc) {
1752: (*split)->ci[0] = -(*split)->nc-1;
1753: (*split)->ci[1] = -(*split)->nc-1;
1754: }
1755: if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1756: else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1757: PetscFunctionReturn(PETSC_SUCCESS);
1758: }
1760: /*@
1761: BVGetSplit - Splits the BV object into two BV objects that share the
1762: internal data, one of them containing the leading columns and the other
1763: one containing the remaining columns.
1765: Collective
1767: Input Parameter:
1768: . bv - the basis vectors context
1770: Output Parameters:
1771: + L - left BV containing leading columns (can be NULL)
1772: - R - right BV containing remaining columns (can be NULL)
1774: Notes:
1775: The columns are split in two sets. The leading columns (including the
1776: constraints) are assigned to the left BV and the remaining columns
1777: are assigned to the right BV. The number of leading columns, as
1778: specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1779: guarantee that both L and R have at least one column).
1781: The returned BV's must be seen as references (not copies) of the input
1782: BV, that is, modifying them will change the entries of bv as well.
1783: The returned BV's must not be destroyed. BVRestoreSplit() must be called
1784: when they are no longer needed.
1786: Pass NULL for any of the output BV's that is not needed.
1788: Level: advanced
1790: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1791: @*/
1792: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1793: {
1794: PetscFunctionBegin;
1797: BVCheckSizes(bv,1);
1798: PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1799: PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1800: bv->lsplit = bv->nc+bv->l;
1801: PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1802: PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1803: if (L) *L = bv->L;
1804: if (R) *R = bv->R;
1805: PetscFunctionReturn(PETSC_SUCCESS);
1806: }
1808: /*@
1809: BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1811: Logically Collective
1813: Input Parameters:
1814: + bv - the basis vectors context
1815: . L - left BV obtained with BVGetSplit()
1816: - R - right BV obtained with BVGetSplit()
1818: Note:
1819: The arguments must match the corresponding call to BVGetSplit().
1821: Level: advanced
1823: .seealso: BVGetSplit()
1824: @*/
1825: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1826: {
1827: PetscFunctionBegin;
1830: BVCheckSizes(bv,1);
1831: PetscCheck(bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1832: PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1833: PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
1834: 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()");
1835: 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()");
1837: PetscTryTypeMethod(bv,restoresplit,L,R);
1838: bv->lsplit = 0;
1839: if (L) *L = NULL;
1840: if (R) *R = NULL;
1841: PetscFunctionReturn(PETSC_SUCCESS);
1842: }
1844: /*@
1845: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
1846: definite inner product.
1848: Logically Collective
1850: Input Parameters:
1851: + bv - basis vectors
1852: - deftol - the tolerance
1854: Options Database Key:
1855: . -bv_definite_tol <deftol> - the tolerance
1857: Notes:
1858: When using a non-standard inner product, see BVSetMatrix(), the solver needs
1859: to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
1860: been declared indefinite, the value z'*B*z must be positive, but due to
1861: rounding error a tiny value may become negative. A tolerance is used to
1862: detect this situation. Likewise, in complex arithmetic z'*B*z should be
1863: real, and we use the same tolerance to check whether a nonzero imaginary part
1864: can be considered negligible.
1866: This function sets this tolerance, which defaults to 10*eps, where eps is
1867: the machine epsilon. The default value should be good for most applications.
1869: Level: advanced
1871: .seealso: BVSetMatrix()
1872: @*/
1873: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
1874: {
1875: PetscFunctionBegin;
1878: if (deftol == (PetscReal)PETSC_DEFAULT) bv->deftol = 10*PETSC_MACHINE_EPSILON;
1879: else {
1880: PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
1881: bv->deftol = deftol;
1882: }
1883: PetscFunctionReturn(PETSC_SUCCESS);
1884: }
1886: /*@
1887: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
1888: inner product.
1890: Not Collective
1892: Input Parameter:
1893: . bv - the basis vectors
1895: Output Parameters:
1896: . deftol - the tolerance
1898: Level: advanced
1900: .seealso: BVSetDefiniteTolerance()
1901: @*/
1902: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
1903: {
1904: PetscFunctionBegin;
1907: *deftol = bv->deftol;
1908: PetscFunctionReturn(PETSC_SUCCESS);
1909: }