1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcbv.h" I*/
16: PetscBool BVRegisterAllCalled = PETSC_FALSE;
17: PetscFunctionList BVList = 0;
19: /*@C
20: BVSetType - Selects the type for the BV object.
22: Logically Collective on BV 24: Input Parameter:
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 ierr,(*r)(BV);
38: PetscBool match;
44: PetscObjectTypeCompare((PetscObject)bv,type,&match);
45: if (match) return(0);
46: PetscStrcmp(type,BVTENSOR,&match);
47: if (match) SETERRQ(PetscObjectComm((PetscObject)bv),1,"Use BVCreateTensor() to create a BV of type tensor");
49: PetscFunctionListFind(BVList,type,&r);
50: if (!r) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested BV type %s",type);
52: if (bv->ops->destroy) { (*bv->ops->destroy)(bv); }
53: PetscMemzero(bv->ops,sizeof(struct _BVOps));
55: PetscObjectChangeTypeName((PetscObject)bv,type);
56: if (bv->n < 0 && bv->N < 0) {
57: bv->ops->create = r;
58: } else {
59: PetscLogEventBegin(BV_Create,bv,0,0,0);
60: (*r)(bv);
61: PetscLogEventEnd(BV_Create,bv,0,0,0);
62: }
63: return(0);
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: . name - name of the type of basis vectors
77: Level: intermediate
79: .seealso: BVSetType()
80: @*/
81: PetscErrorCode BVGetType(BV bv,BVType *type) 82: {
86: *type = ((PetscObject)bv)->type_name;
87: return(0);
88: }
90: /*@
91: BVSetSizes - Sets the local and global sizes, and the number of columns.
93: Collective on BV 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: {
113: PetscInt ma;
119: if (N >= 0 && n > N) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local size %D cannot be larger than global size %D",n,N);
120: if (m <= 0) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %D must be positive",m);
121: if ((bv->n >= 0 || bv->N >= 0) && (bv->n != n || bv->N != N)) SETERRQ4(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change/reset vector sizes to %D local %D global after previously setting them to %D local %D global",n,N,bv->n,bv->N);
122: if (bv->m > 0 && bv->m != m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change the number of columns to %D after previously setting it to %D; use BVResize()",m,bv->m);
123: bv->n = n;
124: bv->N = N;
125: bv->m = m;
126: bv->k = m;
127: if (!bv->t) { /* create template vector and get actual dimensions */
128: VecCreate(PetscObjectComm((PetscObject)bv),&bv->t);
129: VecSetSizes(bv->t,bv->n,bv->N);
130: VecSetFromOptions(bv->t);
131: VecGetSize(bv->t,&bv->N);
132: VecGetLocalSize(bv->t,&bv->n);
133: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
134: MatGetLocalSize(bv->matrix,&ma,NULL);
135: if (bv->n!=ma) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %D does not match that of matrix given at BVSetMatrix %D",bv->n,ma);
136: }
137: }
138: if (bv->ops->create) {
139: PetscLogEventBegin(BV_Create,bv,0,0,0);
140: (*bv->ops->create)(bv);
141: PetscLogEventEnd(BV_Create,bv,0,0,0);
142: bv->ops->create = 0;
143: bv->defersfo = PETSC_FALSE;
144: }
145: return(0);
146: }
148: /*@
149: BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
150: Local and global sizes are specified indirectly by passing a template vector.
152: Collective on BV154: Input Parameters:
155: + bv - the basis vectors
156: . t - the template vector
157: - m - the number of columns
159: Level: beginner
161: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
162: @*/
163: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)164: {
166: PetscInt ma;
173: if (m <= 0) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %D must be positive",m);
174: if (bv->t) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Template vector was already set by a previous call to BVSetSizes/FromVec");
175: VecGetSize(t,&bv->N);
176: VecGetLocalSize(t,&bv->n);
177: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
178: MatGetLocalSize(bv->matrix,&ma,NULL);
179: if (bv->n!=ma) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %D does not match that of matrix given at BVSetMatrix %D",bv->n,ma);
180: }
181: bv->m = m;
182: bv->k = m;
183: bv->t = t;
184: PetscObjectReference((PetscObject)t);
185: if (bv->ops->create) {
186: (*bv->ops->create)(bv);
187: bv->ops->create = 0;
188: bv->defersfo = PETSC_FALSE;
189: }
190: return(0);
191: }
193: /*@
194: BVGetSizes - Returns the local and global sizes, and the number of columns.
196: Not Collective
198: Input Parameter:
199: . bv - the basis vectors
201: Output Parameters:
202: + n - the local size
203: . N - the global size
204: - m - the number of columns
206: Note:
207: Normal usage requires that bv has already been given its sizes, otherwise
208: the call fails. However, this function can also be used to determine if
209: a BV object has been initialized completely (sizes and type). For this,
210: call with n=NULL and N=NULL, then a return value of m=0 indicates that
211: the BV object is not ready for use yet.
213: Level: beginner
215: .seealso: BVSetSizes(), BVSetSizesFromVec()
216: @*/
217: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)218: {
220: if (!bv) {
221: if (m && !n && !N) *m = 0;
222: return(0);
223: }
225: if (n || N) BVCheckSizes(bv,1);
226: if (n) *n = bv->n;
227: if (N) *N = bv->N;
228: if (m) *m = bv->m;
229: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
230: return(0);
231: }
233: /*@
234: BVSetNumConstraints - Set the number of constraints.
236: Logically Collective on BV238: Input Parameters:
239: + V - basis vectors
240: - nc - number of constraints
242: Notes:
243: This function sets the number of constraints to nc and marks all remaining
244: columns as regular. Normal user would call BVInsertConstraints() instead.
246: If nc is smaller than the previously set value, then some of the constraints
247: are discarded. In particular, using nc=0 removes all constraints preserving
248: the content of regular columns.
250: Level: developer
252: .seealso: BVInsertConstraints()
253: @*/
254: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)255: {
257: PetscInt total,diff,i;
258: Vec x,y;
263: if (nc<0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %D) cannot be negative",nc);
265: BVCheckSizes(V,1);
266: if (V->ci[0]!=-V->nc-1 || V->ci[1]!=-V->nc-1) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
268: diff = nc-V->nc;
269: if (!diff) return(0);
270: total = V->nc+V->m;
271: if (total-nc<=0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
272: if (diff<0) { /* lessen constraints, shift contents of BV */
273: for (i=0;i<V->m;i++) {
274: BVGetColumn(V,i,&x);
275: BVGetColumn(V,i+diff,&y);
276: VecCopy(x,y);
277: BVRestoreColumn(V,i,&x);
278: BVRestoreColumn(V,i+diff,&y);
279: }
280: }
281: V->nc = nc;
282: V->ci[0] = -V->nc-1;
283: V->ci[1] = -V->nc-1;
284: V->l = 0;
285: V->m = total-nc;
286: V->k = V->m;
287: PetscObjectStateIncrease((PetscObject)V);
288: return(0);
289: }
291: /*@
292: BVGetNumConstraints - Returns the number of constraints.
294: Not Collective
296: Input Parameter:
297: . bv - the basis vectors
299: Output Parameters:
300: . nc - the number of constraints
302: Level: advanced
304: .seealso: BVGetSizes(), BVInsertConstraints()
305: @*/
306: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)307: {
311: *nc = bv->nc;
312: return(0);
313: }
315: /*@
316: BVResize - Change the number of columns.
318: Collective on BV320: Input Parameters:
321: + bv - the basis vectors
322: . m - the new number of columns
323: - copy - a flag indicating whether current values should be kept
325: Note:
326: Internal storage is reallocated. If the copy flag is set to true, then
327: the contents are copied to the leading part of the new space.
329: Level: advanced
331: .seealso: BVSetSizes(), BVSetSizesFromVec()
332: @*/
333: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)334: {
335: PetscErrorCode ierr;
336: PetscScalar *array;
337: const PetscScalar *omega;
338: Vec v;
345: if (m <= 0) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %D must be positive",m);
346: if (bv->nc && !bv->issplit) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
347: if (bv->m == m) return(0);
349: PetscLogEventBegin(BV_Create,bv,0,0,0);
350: (*bv->ops->resize)(bv,m,copy);
351: VecDestroy(&bv->buffer);
352: BVDestroy(&bv->cached);
353: PetscFree2(bv->h,bv->c);
354: if (bv->omega) {
355: if (bv->cuda) {
356: #if defined(PETSC_HAVE_CUDA)
357: VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v);
358: #else
359: SETERRQ(PetscObjectComm((PetscObject)bv),1,"Something wrong happened");
360: #endif
361: } else {
362: VecCreateSeq(PETSC_COMM_SELF,m,&v);
363: }
364: PetscLogObjectParent((PetscObject)bv,(PetscObject)v);
365: if (copy) {
366: VecGetArray(v,&array);
367: VecGetArrayRead(bv->omega,&omega);
368: PetscMemcpy(array,omega,PetscMin(m,bv->m)*sizeof(PetscScalar));
369: VecRestoreArrayRead(bv->omega,&omega);
370: VecRestoreArray(v,&array);
371: } else {
372: VecSet(v,1.0);
373: }
374: 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: PetscLogEventEnd(BV_Create,bv,0,0,0);
381: PetscObjectStateIncrease((PetscObject)bv);
382: return(0);
383: }
385: /*@
386: BVSetActiveColumns - Specify the columns that will be involved in operations.
388: Logically Collective on BV390: 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: {
418: BVCheckSizes(bv,1);
419: if (k==PETSC_DECIDE || k==PETSC_DEFAULT) {
420: bv->k = bv->m;
421: } else {
422: if (k<0 || k>bv->m) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k. Must be between 0 and m");
423: bv->k = k;
424: }
425: if (l==PETSC_DECIDE || l==PETSC_DEFAULT) {
426: bv->l = 0;
427: } else {
428: if (l<0 || l>bv->k) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l. Must be between 0 and k");
429: bv->l = l;
430: }
431: return(0);
432: }
434: /*@
435: BVGetActiveColumns - Returns the current active dimensions.
437: Not Collective
439: Input Parameter:
440: . bv - the basis vectors context
442: Output Parameter:
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: {
454: if (l) *l = bv->l;
455: if (k) *k = bv->k;
456: return(0);
457: }
459: /*@
460: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
462: Collective on BV464: 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()
487: @*/
488: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)489: {
491: PetscInt m,n;
496: if (B) {
498: MatGetLocalSize(B,&m,&n);
499: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
500: if (bv->m && bv->n!=n) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension BV %D, Mat %D",bv->n,n);
501: }
502: if (B) { PetscObjectReference((PetscObject)B); }
503: MatDestroy(&bv->matrix);
504: bv->matrix = B;
505: bv->indef = indef;
506: PetscObjectStateIncrease((PetscObject)bv);
507: return(0);
508: }
510: /*@
511: BVGetMatrix - Retrieves the matrix representation of the inner product.
513: Not collective, though a parallel Mat may be returned
515: Input Parameter:
516: . bv - the basis vectors context
518: Output Parameter:
519: + B - the matrix of the inner product (may be NULL)
520: - indef - the flag indicating if the matrix is indefinite
522: Level: advanced
524: .seealso: BVSetMatrix()
525: @*/
526: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)527: {
530: if (B) *B = bv->matrix;
531: if (indef) *indef = bv->indef;
532: return(0);
533: }
535: /*@
536: BVApplyMatrix - Multiplies a vector by the matrix representation of the
537: inner product.
539: Neighbor-wise Collective on BV and Vec
541: Input Parameter:
542: + bv - the basis vectors context
543: - x - the vector
545: Output Parameter:
546: . y - the result
548: Note:
549: If no matrix was specified this function copies the vector.
551: Level: advanced
553: .seealso: BVSetMatrix(), BVApplyMatrixBV()
554: @*/
555: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)556: {
563: if (bv->matrix) {
564: BV_IPMatMult(bv,x);
565: VecCopy(bv->Bx,y);
566: } else {
567: VecCopy(x,y);
568: }
569: return(0);
570: }
572: /*@
573: BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
574: of the inner product.
576: Neighbor-wise Collective on BV578: Input Parameter:
579: . X - the basis vectors context
581: Output Parameter:
582: . Y - the basis vectors to store the result (optional)
584: Note:
585: This function computes Y = B*X, where B is the matrix given with
586: BVSetMatrix(). This operation is computed as in BVMatMult().
587: If no matrix was specified, then it just copies Y = X.
589: If no Y is given, the result is stored internally in the cached BV.
591: Level: developer
593: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
594: @*/
595: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)596: {
601: if (Y) {
603: if (X->matrix) {
604: BVMatMult(X,X->matrix,Y);
605: } else {
606: BVCopy(X,Y);
607: }
608: } else {
609: BV_IPMatMultBV(X);
610: }
611: return(0);
612: }
614: /*@
615: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
617: Logically Collective on BV619: Input Parameter:
620: + bv - the basis vectors context
621: - omega - a vector representing the diagonal of the signature matrix
623: Note:
624: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
626: Level: developer
628: .seealso: BVSetMatrix(), BVGetSignature()
629: @*/
630: PetscErrorCode BVSetSignature(BV bv,Vec omega)631: {
632: PetscErrorCode ierr;
633: PetscInt i,n;
634: const PetscScalar *pomega;
635: PetscScalar *intern;
639: BVCheckSizes(bv,1);
643: VecGetSize(omega,&n);
644: if (n!=bv->k) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %D elements, should be %D",n,bv->k);
645: BV_AllocateSignature(bv);
646: if (bv->indef) {
647: VecGetArrayRead(omega,&pomega);
648: VecGetArray(bv->omega,&intern);
649: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
650: VecRestoreArray(bv->omega,&intern);
651: VecRestoreArrayRead(omega,&pomega);
652: } else {
653: PetscInfo(bv,"Ignoring signature because BV is not indefinite\n");
654: }
655: PetscObjectStateIncrease((PetscObject)bv);
656: return(0);
657: }
659: /*@
660: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
662: Not collective
664: Input Parameter:
665: . bv - the basis vectors context
667: Output Parameter:
668: . omega - a vector representing the diagonal of the signature matrix
670: Note:
671: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
673: Level: developer
675: .seealso: BVSetMatrix(), BVSetSignature()
676: @*/
677: PetscErrorCode BVGetSignature(BV bv,Vec omega)678: {
679: PetscErrorCode ierr;
680: PetscInt i,n;
681: PetscScalar *pomega;
682: const PetscScalar *intern;
686: BVCheckSizes(bv,1);
690: VecGetSize(omega,&n);
691: if (n!=bv->k) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %D elements, should be %D",n,bv->k);
692: if (bv->indef && bv->omega) {
693: VecGetArray(omega,&pomega);
694: VecGetArrayRead(bv->omega,&intern);
695: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
696: VecRestoreArrayRead(bv->omega,&intern);
697: VecRestoreArray(omega,&pomega);
698: } else {
699: VecSet(omega,1.0);
700: }
701: return(0);
702: }
704: /*@
705: BVSetBufferVec - Attach a vector object to be used as buffer space for
706: several operations.
708: Collective on BV710: Input Parameters:
711: + bv - the basis vectors context)
712: - buffer - the vector
714: Notes:
715: Use BVGetBufferVec() to retrieve the vector (for example, to free it
716: at the end of the computations).
718: The vector must be sequential of length (nc+m)*m, where m is the number
719: of columns of bv and nc is the number of constraints.
721: Level: developer
723: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
724: @*/
725: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)726: {
728: PetscInt ld,n;
729: PetscMPIInt size;
734: BVCheckSizes(bv,1);
735: VecGetSize(buffer,&n);
736: ld = bv->m+bv->nc;
737: if (n != ld*bv->m) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %d",ld*bv->m);
738: MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size);
739: if (size>1) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
740: 741: PetscObjectReference((PetscObject)buffer);
742: VecDestroy(&bv->buffer);
743: bv->buffer = buffer;
744: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
745: return(0);
746: }
748: /*@
749: BVGetBufferVec - Obtain the buffer vector associated with the BV object.
751: Not Collective
753: Input Parameters:
754: . bv - the basis vectors context
756: Output Parameter:
757: . buffer - vector
759: Notes:
760: The vector is created if not available previously. It is a sequential vector
761: of length (nc+m)*m, where m is the number of columns of bv and nc is the number
762: of constraints.
764: Developer Notes:
765: The buffer vector is viewed as a column-major matrix with leading dimension
766: ld=nc+m, and m columns at most. In the most common usage, it has the structure
767: .vb
768: | | C |
769: |s|---|
770: | | H |
771: .ve
772: where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
773: related to orthogonalization against constraints (first nc rows), and s is the
774: first column that contains scratch values computed during Gram-Schmidt
775: orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
776: store the coefficients.
778: Level: developer
780: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
781: @*/
782: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)783: {
785: PetscInt ld;
790: BVCheckSizes(bv,1);
791: if (!bv->buffer) {
792: ld = bv->m+bv->nc;
793: VecCreate(PETSC_COMM_SELF,&bv->buffer);
794: VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m);
795: VecSetType(bv->buffer,((PetscObject)bv->t)->type_name);
796: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
797: }
798: *buffer = bv->buffer;
799: return(0);
800: }
802: /*@
803: BVSetRandomContext - Sets the PetscRandom object associated with the BV,
804: to be used in operations that need random numbers.
806: Collective on BV808: Input Parameters:
809: + bv - the basis vectors context
810: - rand - the random number generator context
812: Level: advanced
814: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomColumn()
815: @*/
816: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)817: {
824: PetscObjectReference((PetscObject)rand);
825: PetscRandomDestroy(&bv->rand);
826: bv->rand = rand;
827: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
828: return(0);
829: }
831: /*@
832: BVGetRandomContext - Gets the PetscRandom object associated with the BV.
834: Not Collective
836: Input Parameter:
837: . bv - the basis vectors context
839: Output Parameter:
840: . rand - the random number generator context
842: Level: advanced
844: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn()
845: @*/
846: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)847: {
853: if (!bv->rand) {
854: PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand);
855: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
856: if (bv->rrandom) {
857: PetscRandomSetSeed(bv->rand,0x12345678);
858: PetscRandomSeed(bv->rand);
859: }
860: }
861: *rand = bv->rand;
862: return(0);
863: }
865: /*@
866: BVSetFromOptions - Sets BV options from the options database.
868: Collective on BV870: Input Parameter:
871: . bv - the basis vectors context
873: Level: beginner
874: @*/
875: PetscErrorCode BVSetFromOptions(BV bv)876: {
877: PetscErrorCode ierr;
878: char type[256];
879: PetscBool flg1,flg2,flg3,flg4;
880: PetscReal r;
881: BVOrthogType otype;
882: BVOrthogRefineType orefine;
883: BVOrthogBlockType oblock;
887: BVRegisterAll();
888: PetscObjectOptionsBegin((PetscObject)bv);
889: PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVSVEC),type,256,&flg1);
890: if (flg1) {
891: BVSetType(bv,type);
892: } else if (!((PetscObject)bv)->type_name) {
893: BVSetType(bv,BVSVEC);
894: }
896: otype = bv->orthog_type;
897: PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1);
898: orefine = bv->orthog_ref;
899: PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2);
900: r = bv->orthog_eta;
901: PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3);
902: oblock = bv->orthog_block;
903: PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4);
904: if (flg1 || flg2 || flg3 || flg4) { BVSetOrthogonalization(bv,otype,orefine,r,oblock); }
906: PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL);
908: /* undocumented option to generate random vectors that are independent of the number of processes */
909: PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL);
911: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
912: else if (bv->ops->setfromoptions) {
913: (*bv->ops->setfromoptions)(PetscOptionsObject,bv);
914: }
915: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)bv);
916: PetscOptionsEnd();
918: if (!bv->rand) { BVGetRandomContext(bv,&bv->rand); }
919: PetscRandomSetFromOptions(bv->rand);
920: return(0);
921: }
923: /*@
924: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
925: vectors (classical or modified Gram-Schmidt with or without refinement), and
926: for the block-orthogonalization (simultaneous orthogonalization of a set of
927: vectors).
929: Logically Collective on BV931: Input Parameters:
932: + bv - the basis vectors context
933: . type - the method of vector orthogonalization
934: . refine - type of refinement
935: . eta - parameter for selective refinement
936: - block - the method of block orthogonalization
938: Options Database Keys:
939: + -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
940: (default) or mgs for Modified Gram-Schmidt orthogonalization
941: . -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
942: . -bv_orthog_eta <eta> - For setting the value of eta
943: - -bv_orthog_block <block> - Where <block> is the block-orthogonalization method
945: Notes:
946: The default settings work well for most problems.
948: The parameter eta should be a real value between 0 and 1 (or PETSC_DEFAULT).
949: The value of eta is used only when the refinement type is "ifneeded".
951: When using several processors, MGS is likely to result in bad scalability.
953: If the method set for block orthogonalization is GS, then the computation
954: is done column by column with the vector orthogonalization.
956: Level: advanced
958: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType959: @*/
960: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)961: {
968: switch (type) {
969: case BV_ORTHOG_CGS:
970: case BV_ORTHOG_MGS:
971: bv->orthog_type = type;
972: break;
973: default:974: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
975: }
976: switch (refine) {
977: case BV_ORTHOG_REFINE_NEVER:
978: case BV_ORTHOG_REFINE_IFNEEDED:
979: case BV_ORTHOG_REFINE_ALWAYS:
980: bv->orthog_ref = refine;
981: break;
982: default:983: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
984: }
985: if (eta == PETSC_DEFAULT) {
986: bv->orthog_eta = 0.7071;
987: } else {
988: if (eta <= 0.0 || eta > 1.0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
989: bv->orthog_eta = eta;
990: }
991: switch (block) {
992: case BV_ORTHOG_BLOCK_GS:
993: case BV_ORTHOG_BLOCK_CHOL:
994: case BV_ORTHOG_BLOCK_TSQR:
995: case BV_ORTHOG_BLOCK_TSQRCHOL:
996: case BV_ORTHOG_BLOCK_SVQB:
997: bv->orthog_block = block;
998: break;
999: default:1000: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
1001: }
1002: return(0);
1003: }
1005: /*@
1006: BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.
1008: Not Collective
1010: Input Parameter:
1011: . bv - basis vectors context
1013: Output Parameter:
1014: + type - the method of vector orthogonalization
1015: . refine - type of refinement
1016: . eta - parameter for selective refinement
1017: - block - the method of block orthogonalization
1019: Level: advanced
1021: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType1022: @*/
1023: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)1024: {
1027: if (type) *type = bv->orthog_type;
1028: if (refine) *refine = bv->orthog_ref;
1029: if (eta) *eta = bv->orthog_eta;
1030: if (block) *block = bv->orthog_block;
1031: return(0);
1032: }
1034: /*@
1035: BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.
1037: Logically Collective on BV1039: Input Parameters:
1040: + bv - the basis vectors context
1041: - method - the method for the BVMatMult() operation
1043: Options Database Keys:
1044: . -bv_matmult <meth> - choose one of the methods: vecs, mat
1046: Notes:
1047: Allowed values are
1048: + BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1049: . BV_MATMULT_MAT - carry out a MatMatMult() product with a dense matrix
1050: - BV_MATMULT_MAT_SAVE - this case is deprecated
1052: The default is BV_MATMULT_MAT except in the case of BVVECS.
1054: Level: advanced
1056: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType1057: @*/
1058: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)1059: {
1065: switch (method) {
1066: case BV_MATMULT_VECS:
1067: case BV_MATMULT_MAT:
1068: bv->vmm = method;
1069: break;
1070: case BV_MATMULT_MAT_SAVE:
1071: PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n");
1072: bv->vmm = BV_MATMULT_MAT;
1073: break;
1074: default:1075: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1076: }
1077: return(0);
1078: }
1080: /*@
1081: BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.
1083: Not Collective
1085: Input Parameter:
1086: . bv - basis vectors context
1088: Output Parameter:
1089: . method - the method for the BVMatMult() operation
1091: Level: advanced
1093: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType1094: @*/
1095: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)1096: {
1100: *method = bv->vmm;
1101: return(0);
1102: }
1104: /*@
1105: BVGetColumn - Returns a Vec object that contains the entries of the
1106: requested column of the basis vectors object.
1108: Logically Collective on BV1110: Input Parameters:
1111: + bv - the basis vectors context
1112: - j - the index of the requested column
1114: Output Parameter:
1115: . v - vector containing the jth column
1117: Notes:
1118: The returned Vec must be seen as a reference (not a copy) of the BV1119: column, that is, modifying the Vec will change the BV entries as well.
1121: The returned Vec must not be destroyed. BVRestoreColumn() must be
1122: called when it is no longer needed. At most, two columns can be fetched,
1123: that is, this function can only be called twice before the corresponding
1124: BVRestoreColumn() is invoked.
1126: A negative index j selects the i-th constraint, where i=-j. Constraints
1127: should not be modified.
1129: Level: beginner
1131: .seealso: BVRestoreColumn(), BVInsertConstraints()
1132: @*/
1133: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)1134: {
1136: PetscInt l;
1141: BVCheckSizes(bv,1);
1143: if (j<0 && -j>bv->nc) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %D but only %D are available",-j,bv->nc);
1144: if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %D but only %D are available",j,bv->m);
1145: if (j==bv->ci[0] || j==bv->ci[1]) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Column %D already fetched in a previous call to BVGetColumn",j);
1146: l = BVAvailableVec;
1147: if (l==-1) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1148: (*bv->ops->getcolumn)(bv,j,v);
1149: bv->ci[l] = j;
1150: PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]);
1151: PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]);
1152: *v = bv->cv[l];
1153: return(0);
1154: }
1156: /*@
1157: BVRestoreColumn - Restore a column obtained with BVGetColumn().
1159: Logically Collective on BV1161: Input Parameters:
1162: + bv - the basis vectors context
1163: . j - the index of the column
1164: - v - vector obtained with BVGetColumn()
1166: Note:
1167: The arguments must match the corresponding call to BVGetColumn().
1169: Level: beginner
1171: .seealso: BVGetColumn()
1172: @*/
1173: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)1174: {
1175: PetscErrorCode ierr;
1176: PetscObjectId id;
1177: PetscObjectState st;
1178: PetscInt l;
1183: BVCheckSizes(bv,1);
1187: if (j<0 && -j>bv->nc) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %D but only %D are available",-j,bv->nc);
1188: if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %D but only %D are available",j,bv->m);
1189: if (j!=bv->ci[0] && j!=bv->ci[1]) SETERRQ1(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Column %D has not been fetched with a call to BVGetColumn",j);
1190: l = (j==bv->ci[0])? 0: 1;
1191: PetscObjectGetId((PetscObject)*v,&id);
1192: if (id!=bv->id[l]) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1193: PetscObjectStateGet((PetscObject)*v,&st);
1194: if (st!=bv->st[l]) {
1195: PetscObjectStateIncrease((PetscObject)bv);
1196: }
1197: if (bv->ops->restorecolumn) {
1198: (*bv->ops->restorecolumn)(bv,j,v);
1199: } else bv->cv[l] = NULL;
1200: bv->ci[l] = -bv->nc-1;
1201: bv->st[l] = -1;
1202: bv->id[l] = 0;
1203: *v = NULL;
1204: return(0);
1205: }
1207: /*@C
1208: BVGetArray - Returns a pointer to a contiguous array that contains this
1209: processor's portion of the BV data.
1211: Logically Collective on BV1213: Input Parameters:
1214: . bv - the basis vectors context
1216: Output Parameter:
1217: . a - location to put pointer to the array
1219: Notes:
1220: BVRestoreArray() must be called when access to the array is no longer needed.
1221: This operation may imply a data copy, for BV types that do not store
1222: data contiguously in memory.
1224: The pointer will normally point to the first entry of the first column,
1225: but if the BV has constraints then these go before the regular columns.
1227: Level: advanced
1229: .seealso: BVRestoreArray(), BVInsertConstraints()
1230: @*/
1231: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)1232: {
1238: BVCheckSizes(bv,1);
1239: (*bv->ops->getarray)(bv,a);
1240: return(0);
1241: }
1243: /*@C
1244: BVRestoreArray - Restore the BV object after BVGetArray() has been called.
1246: Logically Collective on BV1248: Input Parameters:
1249: + bv - the basis vectors context
1250: - a - location of pointer to array obtained from BVGetArray()
1252: Note:
1253: This operation may imply a data copy, for BV types that do not store
1254: data contiguously in memory.
1256: Level: advanced
1258: .seealso: BVGetColumn()
1259: @*/
1260: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)1261: {
1267: BVCheckSizes(bv,1);
1268: if (bv->ops->restorearray) {
1269: (*bv->ops->restorearray)(bv,a);
1270: }
1271: if (a) *a = NULL;
1272: PetscObjectStateIncrease((PetscObject)bv);
1273: return(0);
1274: }
1276: /*@C
1277: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1278: contains this processor's portion of the BV data.
1280: Not Collective
1282: Input Parameters:
1283: . bv - the basis vectors context
1285: Output Parameter:
1286: . a - location to put pointer to the array
1288: Notes:
1289: BVRestoreArrayRead() must be called when access to the array is no
1290: longer needed. This operation may imply a data copy, for BV types that
1291: do not store data contiguously in memory.
1293: The pointer will normally point to the first entry of the first column,
1294: but if the BV has constraints then these go before the regular columns.
1296: Level: advanced
1298: .seealso: BVRestoreArray(), BVInsertConstraints()
1299: @*/
1300: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)1301: {
1307: BVCheckSizes(bv,1);
1308: (*bv->ops->getarrayread)(bv,a);
1309: return(0);
1310: }
1312: /*@C
1313: BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1314: been called.
1316: Not Collective
1318: Input Parameters:
1319: + bv - the basis vectors context
1320: - a - location of pointer to array obtained from BVGetArrayRead()
1322: Level: advanced
1324: .seealso: BVGetColumn()
1325: @*/
1326: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)1327: {
1333: BVCheckSizes(bv,1);
1334: if (bv->ops->restorearrayread) {
1335: (*bv->ops->restorearrayread)(bv,a);
1336: }
1337: if (a) *a = NULL;
1338: return(0);
1339: }
1341: /*@
1342: BVCreateVec - Creates a new Vec object with the same type and dimensions
1343: as the columns of the basis vectors object.
1345: Collective on BV1347: Input Parameter:
1348: . bv - the basis vectors context
1350: Output Parameter:
1351: . v - the new vector
1353: Note:
1354: The user is responsible of destroying the returned vector.
1356: Level: beginner
1358: .seealso: BVCreateMat()
1359: @*/
1360: PetscErrorCode BVCreateVec(BV bv,Vec *v)1361: {
1366: BVCheckSizes(bv,1);
1368: VecDuplicate(bv->t,v);
1369: return(0);
1370: }
1372: /*@
1373: BVCreateMat - Creates a new Mat object of dense type and copies the contents
1374: of the BV object.
1376: Collective on BV1378: Input Parameter:
1379: . bv - the basis vectors context
1381: Output Parameter:
1382: . A - the new matrix
1384: Notes:
1385: The user is responsible of destroying the returned matrix.
1387: The matrix contains all columns of the BV, not just the active columns.
1389: Level: intermediate
1391: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1392: @*/
1393: PetscErrorCode BVCreateMat(BV bv,Mat *A)1394: {
1395: PetscErrorCode ierr;
1396: PetscScalar *aa;
1397: const PetscScalar *vv;
1401: BVCheckSizes(bv,1);
1404: MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,bv->N,bv->m,NULL,A);
1405: MatDenseGetArray(*A,&aa);
1406: BVGetArrayRead(bv,&vv);
1407: PetscMemcpy(aa,vv,bv->m*bv->n*sizeof(PetscScalar));
1408: BVRestoreArrayRead(bv,&vv);
1409: MatDenseRestoreArray(*A,&aa);
1410: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1411: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1412: return(0);
1413: }
1415: /*@
1416: BVGetMat - Returns a Mat object of dense type that shares the memory of
1417: the BV object.
1419: Collective on BV1421: Input Parameter:
1422: . bv - the basis vectors context
1424: Output Parameter:
1425: . A - the matrix
1427: Notes:
1428: The returned matrix contains only the active columns. If the content of
1429: the Mat is modified, these changes are also done in the BV object. The
1430: user must call BVRestoreMat() when no longer needed.
1432: This operation implies a call to BVGetArray(), which may result in data
1433: copies.
1435: Level: advanced
1437: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1438: @*/
1439: PetscErrorCode BVGetMat(BV bv,Mat *A)1440: {
1442: PetscScalar *vv,*aa;
1443: PetscBool create=PETSC_FALSE;
1444: PetscInt m,cols;
1448: BVCheckSizes(bv,1);
1450: m = bv->k-bv->l;
1451: if (!bv->Aget) create=PETSC_TRUE;
1452: else {
1453: MatDenseGetArray(bv->Aget,&aa);
1454: if (aa) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1455: MatGetSize(bv->Aget,NULL,&cols);
1456: if (cols!=m) {
1457: MatDestroy(&bv->Aget);
1458: create=PETSC_TRUE;
1459: }
1460: }
1461: BVGetArray(bv,&vv);
1462: if (create) {
1463: MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,vv,&bv->Aget); /* pass a pointer to avoid allocation of storage */
1464: MatDensePlaceArray(bv->Aget,NULL); /* replace with a null pointer, the value after BVRestoreMat */
1465: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Aget);
1466: MatAssemblyBegin(bv->Aget,MAT_FINAL_ASSEMBLY);
1467: MatAssemblyEnd(bv->Aget,MAT_FINAL_ASSEMBLY);
1468: }
1469: MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->n); /* set the actual pointer */
1470: *A = bv->Aget;
1471: return(0);
1472: }
1474: /*@
1475: BVRestoreMat - Restores the Mat obtained with BVGetMat().
1477: Collective on BV1479: Input Parameters:
1480: + bv - the basis vectors context
1481: - A - the fetched matrix
1483: Note:
1484: A call to this function must match a previous call of BVGetMat().
1485: The effect is that the contents of the Mat are copied back to the
1486: BV internal data structures.
1488: Level: advanced
1490: .seealso: BVGetMat(), BVRestoreArray()
1491: @*/
1492: PetscErrorCode BVRestoreMat(BV bv,Mat *A)1493: {
1495: PetscScalar *vv,*aa;
1499: BVCheckSizes(bv,1);
1501: if (!bv->Aget) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1502: if (bv->Aget!=*A) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1503: MatDenseGetArray(bv->Aget,&aa);
1504: vv = aa-(bv->nc+bv->l)*bv->n;
1505: MatDenseResetArray(bv->Aget);
1506: BVRestoreArray(bv,&vv);
1507: *A = NULL;
1508: return(0);
1509: }
1511: /*
1512: Copy all user-provided attributes of V to another BV object W
1513: */
1514: PETSC_STATIC_INLINE PetscErrorCode BVDuplicate_Private(BV V,BV W)1515: {
1519: BVSetType(W,((PetscObject)V)->type_name);
1520: W->orthog_type = V->orthog_type;
1521: W->orthog_ref = V->orthog_ref;
1522: W->orthog_eta = V->orthog_eta;
1523: W->orthog_block = V->orthog_block;
1524: if (V->matrix) { PetscObjectReference((PetscObject)V->matrix); }
1525: W->matrix = V->matrix;
1526: W->indef = V->indef;
1527: W->vmm = V->vmm;
1528: W->rrandom = V->rrandom;
1529: if (V->ops->duplicate) { (*V->ops->duplicate)(V,W); }
1530: PetscObjectStateIncrease((PetscObject)W);
1531: return(0);
1532: }
1534: /*@
1535: BVDuplicate - Creates a new basis vector object of the same type and
1536: dimensions as an existing one.
1538: Collective on BV1540: Input Parameter:
1541: . V - basis vectors context
1543: Output Parameter:
1544: . W - location to put the new BV1546: Notes:
1547: The new BV has the same type and dimensions as V, and it shares the same
1548: template vector. Also, the inner product matrix and orthogonalization
1549: options are copied.
1551: BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1552: for the new basis vectors. Use BVCopy() to copy the contents.
1554: Level: intermediate
1556: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1557: @*/
1558: PetscErrorCode BVDuplicate(BV V,BV *W)1559: {
1565: BVCheckSizes(V,1);
1567: BVCreate(PetscObjectComm((PetscObject)V),W);
1568: BVSetSizesFromVec(*W,V->t,V->m);
1569: BVDuplicate_Private(V,*W);
1570: return(0);
1571: }
1573: /*@
1574: BVDuplicateResize - Creates a new basis vector object of the same type and
1575: dimensions as an existing one, but with possibly different number of columns.
1577: Collective on BV1579: Input Parameter:
1580: + V - basis vectors context
1581: - m - the new number of columns
1583: Output Parameter:
1584: . W - location to put the new BV1586: Note:
1587: This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1588: contents of V are not copied to W.
1590: Level: intermediate
1592: .seealso: BVDuplicate(), BVResize()
1593: @*/
1594: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)1595: {
1601: BVCheckSizes(V,1);
1604: BVCreate(PetscObjectComm((PetscObject)V),W);
1605: BVSetSizesFromVec(*W,V->t,m);
1606: BVDuplicate_Private(V,*W);
1607: return(0);
1608: }
1610: /*@
1611: BVGetCachedBV - Returns a BV object stored internally that holds the
1612: result of B*X after a call to BVApplyMatrixBV().
1614: Not collective
1616: Input Parameter:
1617: . bv - the basis vectors context
1619: Output Parameter:
1620: . cached - the cached BV1622: Note:
1623: The cached BV is created if not available previously.
1625: Level: developer
1627: .seealso: BVApplyMatrixBV()
1628: @*/
1629: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)1630: {
1636: BVCheckSizes(bv,1);
1637: if (!bv->cached) {
1638: BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached);
1639: BVSetSizesFromVec(bv->cached,bv->t,bv->m);
1640: BVDuplicate_Private(bv,bv->cached);
1641: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->cached);
1642: }
1643: *cached = bv->cached;
1644: return(0);
1645: }
1647: /*@
1648: BVCopy - Copies a basis vector object into another one, W <- V.
1650: Logically Collective on BV1652: Input Parameter:
1653: . V - basis vectors context
1655: Output Parameter:
1656: . W - the copy
1658: Note:
1659: Both V and W must be distributed in the same manner; local copies are
1660: done. Only active columns (excluding the leading ones) are copied.
1661: In the destination W, columns are overwritten starting from the leading ones.
1662: Constraints are not copied.
1664: Level: beginner
1666: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1667: @*/
1668: PetscErrorCode BVCopy(BV V,BV W)1669: {
1670: PetscErrorCode ierr;
1671: PetscScalar *womega;
1672: const PetscScalar *vomega;
1677: BVCheckSizes(V,1);
1680: BVCheckSizes(W,2);
1682: if (V->n!=W->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, W %D",V->n,W->n);
1683: if (V->k-V->l>W->m-W->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"W has %D non-leading columns, not enough to store %D columns",W->m-W->l,V->k-V->l);
1684: if (!V->n) return(0);
1686: PetscLogEventBegin(BV_Copy,V,W,0,0);
1687: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1688: /* copy signature */
1689: BV_AllocateSignature(W);
1690: VecGetArrayRead(V->omega,&vomega);
1691: VecGetArray(W->omega,&womega);
1692: PetscMemcpy(womega+W->nc+W->l,vomega+V->nc+V->l,(V->k-V->l)*sizeof(PetscScalar));
1693: VecRestoreArray(W->omega,&womega);
1694: VecRestoreArrayRead(V->omega,&vomega);
1695: }
1696: (*V->ops->copy)(V,W);
1697: PetscLogEventEnd(BV_Copy,V,W,0,0);
1698: PetscObjectStateIncrease((PetscObject)W);
1699: return(0);
1700: }
1702: /*@
1703: BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1705: Logically Collective on BV1707: Input Parameter:
1708: + V - basis vectors context
1709: - j - the column number to be copied
1711: Output Parameter:
1712: . w - the copied column
1714: Note:
1715: Both V and w must be distributed in the same manner; local copies are done.
1717: Level: beginner
1719: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1720: @*/
1721: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)1722: {
1724: PetscInt n,N;
1725: Vec z;
1730: BVCheckSizes(V,1);
1735: VecGetSize(w,&N);
1736: VecGetLocalSize(w,&n);
1737: if (N!=V->N || n!=V->n) SETERRQ4(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %D, local %D) do not match BV sizes (global %D, local %D)",N,n,V->N,V->n);
1739: PetscLogEventBegin(BV_Copy,V,w,0,0);
1740: BVGetColumn(V,j,&z);
1741: VecCopy(z,w);
1742: BVRestoreColumn(V,j,&z);
1743: PetscLogEventEnd(BV_Copy,V,w,0,0);
1744: return(0);
1745: }
1747: /*@
1748: BVCopyColumn - Copies the values from one of the columns to another one.
1750: Logically Collective on BV1752: Input Parameter:
1753: + V - basis vectors context
1754: . j - the number of the source column
1755: - i - the number of the destination column
1757: Level: beginner
1759: .seealso: BVCopy(), BVCopyVec()
1760: @*/
1761: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)1762: {
1764: Vec z,w;
1765: PetscScalar *omega;
1770: BVCheckSizes(V,1);
1773: if (j==i) return(0);
1775: PetscLogEventBegin(BV_Copy,V,0,0,0);
1776: if (V->omega) {
1777: VecGetArray(V->omega,&omega);
1778: omega[i] = omega[j];
1779: VecRestoreArray(V->omega,&omega);
1780: }
1781: if (V->ops->copycolumn) {
1782: (*V->ops->copycolumn)(V,j,i);
1783: } else {
1784: BVGetColumn(V,j,&z);
1785: BVGetColumn(V,i,&w);
1786: VecCopy(z,w);
1787: BVRestoreColumn(V,j,&z);
1788: BVRestoreColumn(V,i,&w);
1789: }
1790: PetscLogEventEnd(BV_Copy,V,0,0,0);
1791: PetscObjectStateIncrease((PetscObject)V);
1792: return(0);
1793: }
1795: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)1796: {
1798: PetscInt ncols;
1801: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1802: if (!*split) {
1803: BVCreate(PetscObjectComm((PetscObject)bv),split);
1804: PetscLogObjectParent((PetscObject)bv,(PetscObject)*split);
1805: (*split)->issplit = left? 1: 2;
1806: (*split)->splitparent = bv;
1807: BVSetSizesFromVec(*split,bv->t,ncols);
1808: BVDuplicate_Private(bv,*split);
1809: } else {
1810: BVResize(*split,ncols,PETSC_FALSE);
1811: }
1812: (*split)->l = 0;
1813: (*split)->k = left? bv->l: bv->k-bv->l;
1814: (*split)->nc = left? bv->nc: 0;
1815: (*split)->m = ncols-(*split)->nc;
1816: if ((*split)->nc) {
1817: (*split)->ci[0] = -(*split)->nc-1;
1818: (*split)->ci[1] = -(*split)->nc-1;
1819: }
1820: if (left) {
1821: PetscObjectStateGet((PetscObject)*split,&bv->lstate);
1822: } else {
1823: PetscObjectStateGet((PetscObject)*split,&bv->rstate);
1824: }
1825: return(0);
1826: }
1828: /*@
1829: BVGetSplit - Splits the BV object into two BV objects that share the
1830: internal data, one of them containing the leading columns and the other
1831: one containing the remaining columns.
1833: Logically Collective on BV1835: Input Parameters:
1836: . bv - the basis vectors context
1838: Output Parameters:
1839: + L - left BV containing leading columns (can be NULL)
1840: - R - right BV containing remaining columns (can be NULL)
1842: Notes:
1843: The columns are split in two sets. The leading columns (including the
1844: constraints) are assigned to the left BV and the remaining columns
1845: are assigned to the right BV. The number of leading columns, as
1846: specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1847: guarantee that both L and R have at least one column).
1849: The returned BV's must be seen as references (not copies) of the input
1850: BV, that is, modifying them will change the entries of bv as well.
1851: The returned BV's must not be destroyed. BVRestoreSplit() must be called
1852: when they are no longer needed.
1854: Pass NULL for any of the output BV's that is not needed.
1856: Level: advanced
1858: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1859: @*/
1860: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)1861: {
1867: BVCheckSizes(bv,1);
1868: if (!bv->l) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1869: if (bv->lsplit) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1870: bv->lsplit = bv->nc+bv->l;
1871: BVGetSplit_Private(bv,PETSC_TRUE,&bv->L);
1872: BVGetSplit_Private(bv,PETSC_FALSE,&bv->R);
1873: if (L) *L = bv->L;
1874: if (R) *R = bv->R;
1875: return(0);
1876: }
1878: /*@
1879: BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1881: Logically Collective on BV1883: Input Parameters:
1884: + bv - the basis vectors context
1885: . L - left BV obtained with BVGetSplit()
1886: - R - right BV obtained with BVGetSplit()
1888: Note:
1889: The arguments must match the corresponding call to BVGetSplit().
1891: Level: advanced
1893: .seealso: BVGetSplit()
1894: @*/
1895: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)1896: {
1902: BVCheckSizes(bv,1);
1903: if (!bv->lsplit) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1904: if (L && *L!=bv->L) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1905: if (R && *R!=bv->R) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
1906: if (L && ((*L)->ci[0]>(*L)->nc-1 || (*L)->ci[1]>(*L)->nc-1)) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 2 has unrestored columns, use BVRestoreColumn()");
1907: if (R && ((*R)->ci[0]>(*R)->nc-1 || (*R)->ci[1]>(*R)->nc-1)) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 3 has unrestored columns, use BVRestoreColumn()");
1909: if (bv->ops->restoresplit) {
1910: (*bv->ops->restoresplit)(bv,L,R);
1911: }
1912: bv->lsplit = 0;
1913: if (L) *L = NULL;
1914: if (R) *R = NULL;
1915: return(0);
1916: }