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 = 0;
19: /*@C
20: BVSetType - Selects the type for the BV object.
22: Logically Collective on bv
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;
43: PetscObjectTypeCompare((PetscObject)bv,type,&match);
44: if (match) PetscFunctionReturn(0);
45: PetscStrcmp(type,BVTENSOR,&match);
48: PetscFunctionListFind(BVList,type,&r);
51: if (bv->ops->destroy) (*bv->ops->destroy)(bv);
52: PetscMemzero(bv->ops,sizeof(struct _BVOps));
54: PetscObjectChangeTypeName((PetscObject)bv,type);
55: if (bv->n < 0 && bv->N < 0) {
56: bv->ops->create = r;
57: } else {
58: PetscLogEventBegin(BV_Create,bv,0,0,0);
59: (*r)(bv);
60: PetscLogEventEnd(BV_Create,bv,0,0,0);
61: }
62: PetscFunctionReturn(0);
63: }
65: /*@C
66: BVGetType - Gets the BV type name (as a string) from the BV context.
68: Not Collective
70: Input Parameter:
71: . bv - the basis vectors context
73: Output Parameter:
74: . type - name of the type of basis vectors
76: Level: intermediate
78: .seealso: BVSetType()
79: @*/
80: PetscErrorCode BVGetType(BV bv,BVType *type) 81: {
84: *type = ((PetscObject)bv)->type_name;
85: PetscFunctionReturn(0);
86: }
88: /*@
89: BVSetSizes - Sets the local and global sizes, and the number of columns.
91: Collective on bv
93: Input Parameters:
94: + bv - the basis vectors
95: . n - the local size (or PETSC_DECIDE to have it set)
96: . N - the global size (or PETSC_DECIDE)
97: - m - the number of columns
99: Notes:
100: n and N cannot be both PETSC_DECIDE.
101: If one processor calls this with N of PETSC_DECIDE then all processors must,
102: otherwise the program will hang.
104: Level: beginner
106: .seealso: BVSetSizesFromVec(), BVGetSizes(), BVResize()
107: @*/
108: PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)109: {
110: PetscInt ma;
119: bv->n = n;
120: bv->N = N;
121: bv->m = m;
122: bv->k = m;
123: if (!bv->t) { /* create template vector and get actual dimensions */
124: VecCreate(PetscObjectComm((PetscObject)bv),&bv->t);
125: VecSetSizes(bv->t,bv->n,bv->N);
126: VecSetFromOptions(bv->t);
127: VecGetSize(bv->t,&bv->N);
128: VecGetLocalSize(bv->t,&bv->n);
129: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
130: MatGetLocalSize(bv->matrix,&ma,NULL);
132: }
133: }
134: if (bv->ops->create) {
135: PetscLogEventBegin(BV_Create,bv,0,0,0);
136: (*bv->ops->create)(bv);
137: PetscLogEventEnd(BV_Create,bv,0,0,0);
138: bv->ops->create = 0;
139: bv->defersfo = PETSC_FALSE;
140: }
141: PetscFunctionReturn(0);
142: }
144: /*@
145: BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
146: Local and global sizes are specified indirectly by passing a template vector.
148: Collective on bv
150: Input Parameters:
151: + bv - the basis vectors
152: . t - the template vector
153: - m - the number of columns
155: Level: beginner
157: .seealso: BVSetSizes(), BVGetSizes(), BVResize()
158: @*/
159: PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)160: {
161: PetscInt ma;
169: VecGetSize(t,&bv->N);
170: VecGetLocalSize(t,&bv->n);
171: if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
172: MatGetLocalSize(bv->matrix,&ma,NULL);
174: }
175: bv->m = m;
176: bv->k = m;
177: bv->t = t;
178: PetscObjectReference((PetscObject)t);
179: if (bv->ops->create) {
180: (*bv->ops->create)(bv);
181: bv->ops->create = 0;
182: bv->defersfo = PETSC_FALSE;
183: }
184: PetscFunctionReturn(0);
185: }
187: /*@
188: BVGetSizes - Returns the local and global sizes, and the number of columns.
190: Not Collective
192: Input Parameter:
193: . bv - the basis vectors
195: Output Parameters:
196: + n - the local size
197: . N - the global size
198: - m - the number of columns
200: Note:
201: Normal usage requires that bv has already been given its sizes, otherwise
202: the call fails. However, this function can also be used to determine if
203: a BV object has been initialized completely (sizes and type). For this,
204: call with n=NULL and N=NULL, then a return value of m=0 indicates that
205: the BV object is not ready for use yet.
207: Level: beginner
209: .seealso: BVSetSizes(), BVSetSizesFromVec()
210: @*/
211: PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)212: {
213: if (!bv) {
214: if (m && !n && !N) *m = 0;
215: PetscFunctionReturn(0);
216: }
218: if (n || N) BVCheckSizes(bv,1);
219: if (n) *n = bv->n;
220: if (N) *N = bv->N;
221: if (m) *m = bv->m;
222: if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
223: PetscFunctionReturn(0);
224: }
226: /*@
227: BVSetNumConstraints - Set the number of constraints.
229: Logically Collective on V
231: Input Parameters:
232: + V - basis vectors
233: - nc - number of constraints
235: Notes:
236: This function sets the number of constraints to nc and marks all remaining
237: columns as regular. Normal user would call BVInsertConstraints() instead.
239: If nc is smaller than the previously set value, then some of the constraints
240: are discarded. In particular, using nc=0 removes all constraints preserving
241: the content of regular columns.
243: Level: developer
245: .seealso: BVInsertConstraints()
246: @*/
247: PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)248: {
249: PetscInt total,diff,i;
250: Vec x,y;
256: BVCheckSizes(V,1);
259: diff = nc-V->nc;
260: if (!diff) PetscFunctionReturn(0);
261: total = V->nc+V->m;
263: if (diff<0) { /* lessen constraints, shift contents of BV */
264: for (i=0;i<V->m;i++) {
265: BVGetColumn(V,i,&x);
266: BVGetColumn(V,i+diff,&y);
267: VecCopy(x,y);
268: BVRestoreColumn(V,i,&x);
269: BVRestoreColumn(V,i+diff,&y);
270: }
271: }
272: V->nc = nc;
273: V->ci[0] = -V->nc-1;
274: V->ci[1] = -V->nc-1;
275: V->m = total-nc;
276: V->l = PetscMin(V->l,V->m);
277: V->k = PetscMin(V->k,V->m);
278: PetscObjectStateIncrease((PetscObject)V);
279: PetscFunctionReturn(0);
280: }
282: /*@
283: BVGetNumConstraints - Returns the number of constraints.
285: Not Collective
287: Input Parameter:
288: . bv - the basis vectors
290: Output Parameters:
291: . nc - the number of constraints
293: Level: advanced
295: .seealso: BVGetSizes(), BVInsertConstraints()
296: @*/
297: PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)298: {
301: *nc = bv->nc;
302: PetscFunctionReturn(0);
303: }
305: /*@
306: BVResize - Change the number of columns.
308: Collective on bv
310: Input Parameters:
311: + bv - the basis vectors
312: . m - the new number of columns
313: - copy - a flag indicating whether current values should be kept
315: Note:
316: Internal storage is reallocated. If the copy flag is set to true, then
317: the contents are copied to the leading part of the new space.
319: Level: advanced
321: .seealso: BVSetSizes(), BVSetSizesFromVec()
322: @*/
323: PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)324: {
325: PetscScalar *array;
326: const PetscScalar *omega;
327: Vec v;
335: if (bv->m == m) PetscFunctionReturn(0);
336: BVCheckOp(bv,1,resize);
338: PetscLogEventBegin(BV_Create,bv,0,0,0);
339: (*bv->ops->resize)(bv,m,copy);
340: VecDestroy(&bv->buffer);
341: BVDestroy(&bv->cached);
342: PetscFree2(bv->h,bv->c);
343: if (bv->omega) {
344: if (bv->cuda) {
345: #if defined(PETSC_HAVE_CUDA)
346: VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v);
347: #else
348: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
349: #endif
350: } else VecCreateSeq(PETSC_COMM_SELF,m,&v);
351: PetscLogObjectParent((PetscObject)bv,(PetscObject)v);
352: if (copy) {
353: VecGetArray(v,&array);
354: VecGetArrayRead(bv->omega,&omega);
355: PetscArraycpy(array,omega,PetscMin(m,bv->m));
356: VecRestoreArrayRead(bv->omega,&omega);
357: VecRestoreArray(v,&array);
358: } else VecSet(v,1.0);
359: VecDestroy(&bv->omega);
360: bv->omega = v;
361: }
362: bv->m = m;
363: bv->k = PetscMin(bv->k,m);
364: bv->l = PetscMin(bv->l,m);
365: PetscLogEventEnd(BV_Create,bv,0,0,0);
366: PetscObjectStateIncrease((PetscObject)bv);
367: PetscFunctionReturn(0);
368: }
370: /*@
371: BVSetActiveColumns - Specify the columns that will be involved in operations.
373: Logically Collective on bv
375: Input Parameters:
376: + bv - the basis vectors context
377: . l - number of leading columns
378: - k - number of active columns
380: Notes:
381: In operations such as BVMult() or BVDot(), only the first k columns are
382: considered. This is useful when the BV is filled from left to right, so
383: the last m-k columns do not have relevant information.
385: Also in operations such as BVMult() or BVDot(), the first l columns are
386: normally not included in the computation. See the manpage of each
387: operation.
389: In orthogonalization operations, the first l columns are treated
390: differently, they participate in the orthogonalization but the computed
391: coefficients are not stored.
393: Level: intermediate
395: .seealso: BVGetActiveColumns(), BVSetSizes()
396: @*/
397: PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)398: {
402: BVCheckSizes(bv,1);
403: if (PetscUnlikely(k==PETSC_DECIDE || k==PETSC_DEFAULT)) {
404: bv->k = bv->m;
405: } else {
407: bv->k = k;
408: }
409: if (PetscUnlikely(l==PETSC_DECIDE || l==PETSC_DEFAULT)) {
410: bv->l = 0;
411: } else {
413: bv->l = l;
414: }
415: PetscFunctionReturn(0);
416: }
418: /*@
419: BVGetActiveColumns - Returns the current active dimensions.
421: Not Collective
423: Input Parameter:
424: . bv - the basis vectors context
426: Output Parameters:
427: + l - number of leading columns
428: - k - number of active columns
430: Level: intermediate
432: .seealso: BVSetActiveColumns()
433: @*/
434: PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)435: {
437: if (l) *l = bv->l;
438: if (k) *k = bv->k;
439: PetscFunctionReturn(0);
440: }
442: /*@
443: BVSetMatrix - Specifies the inner product to be used in orthogonalization.
445: Collective on bv
447: Input Parameters:
448: + bv - the basis vectors context
449: . B - a symmetric matrix (may be NULL)
450: - indef - a flag indicating if the matrix is indefinite
452: Notes:
453: This is used to specify a non-standard inner product, whose matrix
454: representation is given by B. Then, all inner products required during
455: orthogonalization are computed as (x,y)_B=y^H*B*x rather than the
456: standard form (x,y)=y^H*x.
458: Matrix B must be real symmetric (or complex Hermitian). A genuine inner
459: product requires that B is also positive (semi-)definite. However, we
460: also allow for an indefinite B (setting indef=PETSC_TRUE), in which
461: case the orthogonalization uses an indefinite inner product.
463: This affects operations BVDot(), BVNorm(), BVOrthogonalize(), and variants.
465: Setting B=NULL has the same effect as if the identity matrix was passed.
467: Level: advanced
469: .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
470: @*/
471: PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)472: {
473: PetscInt m,n;
477: if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
478: if (B) {
480: MatGetLocalSize(B,&m,&n);
483: }
484: if (B) PetscObjectReference((PetscObject)B);
485: MatDestroy(&bv->matrix);
486: bv->matrix = B;
487: bv->indef = indef;
488: PetscObjectStateIncrease((PetscObject)bv);
489: if (bv->Bx) PetscObjectStateIncrease((PetscObject)bv->Bx);
490: if (bv->cached) PetscObjectStateIncrease((PetscObject)bv->cached);
491: }
492: PetscFunctionReturn(0);
493: }
495: /*@
496: BVGetMatrix - Retrieves the matrix representation of the inner product.
498: Not collective, though a parallel Mat may be returned
500: Input Parameter:
501: . bv - the basis vectors context
503: Output Parameters:
504: + B - the matrix of the inner product (may be NULL)
505: - indef - the flag indicating if the matrix is indefinite
507: Level: advanced
509: .seealso: BVSetMatrix()
510: @*/
511: PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)512: {
514: if (B) *B = bv->matrix;
515: if (indef) *indef = bv->indef;
516: PetscFunctionReturn(0);
517: }
519: /*@
520: BVApplyMatrix - Multiplies a vector by the matrix representation of the
521: inner product.
523: Neighbor-wise Collective on bv
525: Input Parameters:
526: + bv - the basis vectors context
527: - x - the vector
529: Output Parameter:
530: . y - the result
532: Note:
533: If no matrix was specified this function copies the vector.
535: Level: advanced
537: .seealso: BVSetMatrix(), BVApplyMatrixBV()
538: @*/
539: PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)540: {
544: if (bv->matrix) {
545: BV_IPMatMult(bv,x);
546: VecCopy(bv->Bx,y);
547: } else VecCopy(x,y);
548: PetscFunctionReturn(0);
549: }
551: /*@
552: BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
553: of the inner product.
555: Neighbor-wise Collective on X
557: Input Parameter:
558: . X - the basis vectors context
560: Output Parameter:
561: . Y - the basis vectors to store the result (optional)
563: Note:
564: This function computes Y = B*X, where B is the matrix given with
565: BVSetMatrix(). This operation is computed as in BVMatMult().
566: If no matrix was specified, then it just copies Y = X.
568: If no Y is given, the result is stored internally in the cached BV.
570: Level: developer
572: .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
573: @*/
574: PetscErrorCode BVApplyMatrixBV(BV X,BV Y)575: {
577: if (Y) {
579: if (X->matrix) BVMatMult(X,X->matrix,Y);
580: else BVCopy(X,Y);
581: } else BV_IPMatMultBV(X);
582: PetscFunctionReturn(0);
583: }
585: /*@
586: BVSetSignature - Sets the signature matrix to be used in orthogonalization.
588: Logically Collective on bv
590: Input Parameters:
591: + bv - the basis vectors context
592: - omega - a vector representing the diagonal of the signature matrix
594: Note:
595: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
597: Level: developer
599: .seealso: BVSetMatrix(), BVGetSignature()
600: @*/
601: PetscErrorCode BVSetSignature(BV bv,Vec omega)602: {
603: PetscInt i,n;
604: const PetscScalar *pomega;
605: PetscScalar *intern;
608: BVCheckSizes(bv,1);
612: VecGetSize(omega,&n);
614: BV_AllocateSignature(bv);
615: if (bv->indef) {
616: VecGetArrayRead(omega,&pomega);
617: VecGetArray(bv->omega,&intern);
618: for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
619: VecRestoreArray(bv->omega,&intern);
620: VecRestoreArrayRead(omega,&pomega);
621: } else PetscInfo(bv,"Ignoring signature because BV is not indefinite\n");
622: PetscObjectStateIncrease((PetscObject)bv);
623: PetscFunctionReturn(0);
624: }
626: /*@
627: BVGetSignature - Retrieves the signature matrix from last orthogonalization.
629: Not collective
631: Input Parameter:
632: . bv - the basis vectors context
634: Output Parameter:
635: . omega - a vector representing the diagonal of the signature matrix
637: Note:
638: The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
640: Level: developer
642: .seealso: BVSetMatrix(), BVSetSignature()
643: @*/
644: PetscErrorCode BVGetSignature(BV bv,Vec omega)645: {
646: PetscInt i,n;
647: PetscScalar *pomega;
648: const PetscScalar *intern;
651: BVCheckSizes(bv,1);
655: VecGetSize(omega,&n);
657: if (bv->indef && bv->omega) {
658: VecGetArray(omega,&pomega);
659: VecGetArrayRead(bv->omega,&intern);
660: for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
661: VecRestoreArrayRead(bv->omega,&intern);
662: VecRestoreArray(omega,&pomega);
663: } else VecSet(omega,1.0);
664: PetscFunctionReturn(0);
665: }
667: /*@
668: BVSetBufferVec - Attach a vector object to be used as buffer space for
669: several operations.
671: Collective on bv
673: Input Parameters:
674: + bv - the basis vectors context)
675: - buffer - the vector
677: Notes:
678: Use BVGetBufferVec() to retrieve the vector (for example, to free it
679: at the end of the computations).
681: The vector must be sequential of length (nc+m)*m, where m is the number
682: of columns of bv and nc is the number of constraints.
684: Level: developer
686: .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
687: @*/
688: PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)689: {
690: PetscInt ld,n;
691: PetscMPIInt size;
695: BVCheckSizes(bv,1);
696: VecGetSize(buffer,&n);
697: ld = bv->m+bv->nc;
699: MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size);
702: PetscObjectReference((PetscObject)buffer);
703: VecDestroy(&bv->buffer);
704: bv->buffer = buffer;
705: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
706: PetscFunctionReturn(0);
707: }
709: /*@
710: BVGetBufferVec - Obtain the buffer vector associated with the BV object.
712: Not Collective, but Vec returned is parallel if BV is parallel
714: Input Parameters:
715: . bv - the basis vectors context
717: Output Parameter:
718: . buffer - vector
720: Notes:
721: The vector is created if not available previously. It is a sequential vector
722: of length (nc+m)*m, where m is the number of columns of bv and nc is the number
723: of constraints.
725: Developer Notes:
726: The buffer vector is viewed as a column-major matrix with leading dimension
727: ld=nc+m, and m columns at most. In the most common usage, it has the structure
728: .vb
729: | | C |
730: |s|---|
731: | | H |
732: .ve
733: where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
734: related to orthogonalization against constraints (first nc rows), and s is the
735: first column that contains scratch values computed during Gram-Schmidt
736: orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
737: store the coefficients.
739: Level: developer
741: .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
742: @*/
743: PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)744: {
745: PetscInt ld;
749: BVCheckSizes(bv,1);
750: if (!bv->buffer) {
751: ld = bv->m+bv->nc;
752: VecCreate(PETSC_COMM_SELF,&bv->buffer);
753: VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m);
754: VecSetType(bv->buffer,((PetscObject)bv->t)->type_name);
755: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->buffer);
756: }
757: *buffer = bv->buffer;
758: PetscFunctionReturn(0);
759: }
761: /*@
762: BVSetRandomContext - Sets the PetscRandom object associated with the BV,
763: to be used in operations that need random numbers.
765: Collective on bv
767: Input Parameters:
768: + bv - the basis vectors context
769: - rand - the random number generator context
771: Level: advanced
773: .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
774: @*/
775: PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)776: {
780: PetscObjectReference((PetscObject)rand);
781: PetscRandomDestroy(&bv->rand);
782: bv->rand = rand;
783: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
784: PetscFunctionReturn(0);
785: }
787: /*@
788: BVGetRandomContext - Gets the PetscRandom object associated with the BV.
790: Not Collective
792: Input Parameter:
793: . bv - the basis vectors context
795: Output Parameter:
796: . rand - the random number generator context
798: Level: advanced
800: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
801: @*/
802: PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)803: {
806: if (!bv->rand) {
807: PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand);
808: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->rand);
809: if (bv->cuda) PetscRandomSetType(bv->rand,PETSCCURAND);
810: if (bv->sfocalled) PetscRandomSetFromOptions(bv->rand);
811: if (bv->rrandom) {
812: PetscRandomSetSeed(bv->rand,0x12345678);
813: PetscRandomSeed(bv->rand);
814: }
815: }
816: *rand = bv->rand;
817: PetscFunctionReturn(0);
818: }
820: /*@
821: BVSetFromOptions - Sets BV options from the options database.
823: Collective on bv
825: Input Parameter:
826: . bv - the basis vectors context
828: Level: beginner
830: .seealso: BVSetOptionsPrefix()
831: @*/
832: PetscErrorCode BVSetFromOptions(BV bv)833: {
834: PetscErrorCode ierr;
835: char type[256];
836: PetscBool flg1,flg2,flg3,flg4;
837: PetscReal r;
838: BVOrthogType otype;
839: BVOrthogRefineType orefine;
840: BVOrthogBlockType oblock;
843: BVRegisterAll();
844: ierr = PetscObjectOptionsBegin((PetscObject)bv);
845: PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVSVEC),type,sizeof(type),&flg1);
846: if (flg1) BVSetType(bv,type);
847: else if (!((PetscObject)bv)->type_name) BVSetType(bv,BVSVEC);
849: otype = bv->orthog_type;
850: PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1);
851: orefine = bv->orthog_ref;
852: PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2);
853: r = bv->orthog_eta;
854: PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3);
855: oblock = bv->orthog_block;
856: PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4);
857: if (flg1 || flg2 || flg3 || flg4) BVSetOrthogonalization(bv,otype,orefine,r,oblock);
859: PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL);
861: PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1);
862: if (flg1) BVSetDefiniteTolerance(bv,r);
864: /* undocumented option to generate random vectors that are independent of the number of processes */
865: PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL);
867: if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
868: else if (bv->ops->setfromoptions) (*bv->ops->setfromoptions)(PetscOptionsObject,bv);
869: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)bv);
870: ierr = PetscOptionsEnd();
871: bv->sfocalled = PETSC_TRUE;
872: PetscFunctionReturn(0);
873: }
875: /*@
876: BVSetOrthogonalization - Specifies the method used for the orthogonalization of
877: vectors (classical or modified Gram-Schmidt with or without refinement), and
878: for the block-orthogonalization (simultaneous orthogonalization of a set of
879: vectors).
881: Logically Collective on bv
883: Input Parameters:
884: + bv - the basis vectors context
885: . type - the method of vector orthogonalization
886: . refine - type of refinement
887: . eta - parameter for selective refinement
888: - block - the method of block orthogonalization
890: Options Database Keys:
891: + -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
892: (default) or mgs for Modified Gram-Schmidt orthogonalization
893: . -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
894: . -bv_orthog_eta <eta> - For setting the value of eta
895: - -bv_orthog_block <block> - Where <block> is the block-orthogonalization method
897: Notes:
898: The default settings work well for most problems.
900: The parameter eta should be a real value between 0 and 1 (or PETSC_DEFAULT).
901: The value of eta is used only when the refinement type is "ifneeded".
903: When using several processors, MGS is likely to result in bad scalability.
905: If the method set for block orthogonalization is GS, then the computation
906: is done column by column with the vector orthogonalization.
908: Level: advanced
910: .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType911: @*/
912: PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)913: {
919: switch (type) {
920: case BV_ORTHOG_CGS:
921: case BV_ORTHOG_MGS:
922: bv->orthog_type = type;
923: break;
924: default:925: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
926: }
927: switch (refine) {
928: case BV_ORTHOG_REFINE_NEVER:
929: case BV_ORTHOG_REFINE_IFNEEDED:
930: case BV_ORTHOG_REFINE_ALWAYS:
931: bv->orthog_ref = refine;
932: break;
933: default:934: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
935: }
936: if (eta == PETSC_DEFAULT) {
937: bv->orthog_eta = 0.7071;
938: } else {
940: bv->orthog_eta = eta;
941: }
942: switch (block) {
943: case BV_ORTHOG_BLOCK_GS:
944: case BV_ORTHOG_BLOCK_CHOL:
945: case BV_ORTHOG_BLOCK_TSQR:
946: case BV_ORTHOG_BLOCK_TSQRCHOL:
947: case BV_ORTHOG_BLOCK_SVQB:
948: bv->orthog_block = block;
949: break;
950: default:951: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
952: }
953: PetscFunctionReturn(0);
954: }
956: /*@
957: BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.
959: Not Collective
961: Input Parameter:
962: . bv - basis vectors context
964: Output Parameters:
965: + type - the method of vector orthogonalization
966: . refine - type of refinement
967: . eta - parameter for selective refinement
968: - block - the method of block orthogonalization
970: Level: advanced
972: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType973: @*/
974: PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)975: {
977: if (type) *type = bv->orthog_type;
978: if (refine) *refine = bv->orthog_ref;
979: if (eta) *eta = bv->orthog_eta;
980: if (block) *block = bv->orthog_block;
981: PetscFunctionReturn(0);
982: }
984: /*@
985: BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.
987: Logically Collective on bv
989: Input Parameters:
990: + bv - the basis vectors context
991: - method - the method for the BVMatMult() operation
993: Options Database Keys:
994: . -bv_matmult <meth> - choose one of the methods: vecs, mat
996: Notes:
997: Allowed values are
998: + BV_MATMULT_VECS - perform a matrix-vector multiply per each column
999: . BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1000: - BV_MATMULT_MAT_SAVE - this case is deprecated
1002: The default is BV_MATMULT_MAT except in the case of BVVECS.
1004: Level: advanced
1006: .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType1007: @*/
1008: PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)1009: {
1012: switch (method) {
1013: case BV_MATMULT_VECS:
1014: case BV_MATMULT_MAT:
1015: bv->vmm = method;
1016: break;
1017: case BV_MATMULT_MAT_SAVE:
1018: PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n");
1019: bv->vmm = BV_MATMULT_MAT;
1020: break;
1021: default:1022: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1023: }
1024: PetscFunctionReturn(0);
1025: }
1027: /*@
1028: BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.
1030: Not Collective
1032: Input Parameter:
1033: . bv - basis vectors context
1035: Output Parameter:
1036: . method - the method for the BVMatMult() operation
1038: Level: advanced
1040: .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType1041: @*/
1042: PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)1043: {
1046: *method = bv->vmm;
1047: PetscFunctionReturn(0);
1048: }
1050: /*@
1051: BVGetColumn - Returns a Vec object that contains the entries of the
1052: requested column of the basis vectors object.
1054: Logically Collective on bv
1056: Input Parameters:
1057: + bv - the basis vectors context
1058: - j - the index of the requested column
1060: Output Parameter:
1061: . v - vector containing the jth column
1063: Notes:
1064: The returned Vec must be seen as a reference (not a copy) of the BV1065: column, that is, modifying the Vec will change the BV entries as well.
1067: The returned Vec must not be destroyed. BVRestoreColumn() must be
1068: called when it is no longer needed. At most, two columns can be fetched,
1069: that is, this function can only be called twice before the corresponding
1070: BVRestoreColumn() is invoked.
1072: A negative index j selects the i-th constraint, where i=-j. Constraints
1073: should not be modified.
1075: Level: beginner
1077: .seealso: BVRestoreColumn(), BVInsertConstraints()
1078: @*/
1079: PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)1080: {
1081: PetscInt l;
1085: BVCheckSizes(bv,1);
1086: BVCheckOp(bv,1,getcolumn);
1091: l = BVAvailableVec;
1093: (*bv->ops->getcolumn)(bv,j,v);
1094: bv->ci[l] = j;
1095: PetscObjectStateGet((PetscObject)bv->cv[l],&bv->st[l]);
1096: PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]);
1097: *v = bv->cv[l];
1098: PetscFunctionReturn(0);
1099: }
1101: /*@
1102: BVRestoreColumn - Restore a column obtained with BVGetColumn().
1104: Logically Collective on bv
1106: Input Parameters:
1107: + bv - the basis vectors context
1108: . j - the index of the column
1109: - v - vector obtained with BVGetColumn()
1111: Note:
1112: The arguments must match the corresponding call to BVGetColumn().
1114: Level: beginner
1116: .seealso: BVGetColumn()
1117: @*/
1118: PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)1119: {
1120: PetscObjectId id;
1121: PetscObjectState st;
1122: PetscInt l;
1126: BVCheckSizes(bv,1);
1133: l = (j==bv->ci[0])? 0: 1;
1134: PetscObjectGetId((PetscObject)*v,&id);
1136: PetscObjectStateGet((PetscObject)*v,&st);
1137: if (st!=bv->st[l]) PetscObjectStateIncrease((PetscObject)bv);
1138: if (bv->ops->restorecolumn) (*bv->ops->restorecolumn)(bv,j,v);
1139: else bv->cv[l] = NULL;
1140: bv->ci[l] = -bv->nc-1;
1141: bv->st[l] = -1;
1142: bv->id[l] = 0;
1143: *v = NULL;
1144: PetscFunctionReturn(0);
1145: }
1147: /*@C
1148: BVGetArray - Returns a pointer to a contiguous array that contains this
1149: processor's portion of the BV data.
1151: Logically Collective on bv
1153: Input Parameters:
1154: . bv - the basis vectors context
1156: Output Parameter:
1157: . a - location to put pointer to the array
1159: Notes:
1160: BVRestoreArray() must be called when access to the array is no longer needed.
1161: This operation may imply a data copy, for BV types that do not store
1162: data contiguously in memory.
1164: The pointer will normally point to the first entry of the first column,
1165: but if the BV has constraints then these go before the regular columns.
1167: Level: advanced
1169: .seealso: BVRestoreArray(), BVInsertConstraints()
1170: @*/
1171: PetscErrorCode BVGetArray(BV bv,PetscScalar **a)1172: {
1175: BVCheckSizes(bv,1);
1176: BVCheckOp(bv,1,getarray);
1177: (*bv->ops->getarray)(bv,a);
1178: PetscFunctionReturn(0);
1179: }
1181: /*@C
1182: BVRestoreArray - Restore the BV object after BVGetArray() has been called.
1184: Logically Collective on bv
1186: Input Parameters:
1187: + bv - the basis vectors context
1188: - a - location of pointer to array obtained from BVGetArray()
1190: Note:
1191: This operation may imply a data copy, for BV types that do not store
1192: data contiguously in memory.
1194: Level: advanced
1196: .seealso: BVGetColumn()
1197: @*/
1198: PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)1199: {
1202: BVCheckSizes(bv,1);
1203: if (bv->ops->restorearray) (*bv->ops->restorearray)(bv,a);
1204: if (a) *a = NULL;
1205: PetscObjectStateIncrease((PetscObject)bv);
1206: PetscFunctionReturn(0);
1207: }
1209: /*@C
1210: BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1211: contains this processor's portion of the BV data.
1213: Not Collective
1215: Input Parameters:
1216: . bv - the basis vectors context
1218: Output Parameter:
1219: . a - location to put pointer to the array
1221: Notes:
1222: BVRestoreArrayRead() must be called when access to the array is no
1223: longer needed. This operation may imply a data copy, for BV types that
1224: do not store data contiguously in memory.
1226: The pointer will normally point to the first entry of the first column,
1227: but if the BV has constraints then these go before the regular columns.
1229: Level: advanced
1231: .seealso: BVRestoreArray(), BVInsertConstraints()
1232: @*/
1233: PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)1234: {
1237: BVCheckSizes(bv,1);
1238: BVCheckOp(bv,1,getarrayread);
1239: (*bv->ops->getarrayread)(bv,a);
1240: PetscFunctionReturn(0);
1241: }
1243: /*@C
1244: BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1245: been called.
1247: Not Collective
1249: Input Parameters:
1250: + bv - the basis vectors context
1251: - a - location of pointer to array obtained from BVGetArrayRead()
1253: Level: advanced
1255: .seealso: BVGetColumn()
1256: @*/
1257: PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)1258: {
1261: BVCheckSizes(bv,1);
1262: if (bv->ops->restorearrayread) (*bv->ops->restorearrayread)(bv,a);
1263: if (a) *a = NULL;
1264: PetscFunctionReturn(0);
1265: }
1267: /*@
1268: BVCreateVec - Creates a new Vec object with the same type and dimensions
1269: as the columns of the basis vectors object.
1271: Collective on bv
1273: Input Parameter:
1274: . bv - the basis vectors context
1276: Output Parameter:
1277: . v - the new vector
1279: Note:
1280: The user is responsible of destroying the returned vector.
1282: Level: beginner
1284: .seealso: BVCreateMat()
1285: @*/
1286: PetscErrorCode BVCreateVec(BV bv,Vec *v)1287: {
1289: BVCheckSizes(bv,1);
1291: VecDuplicate(bv->t,v);
1292: PetscFunctionReturn(0);
1293: }
1295: /*@
1296: BVCreateMat - Creates a new Mat object of dense type and copies the contents
1297: of the BV object.
1299: Collective on bv
1301: Input Parameter:
1302: . bv - the basis vectors context
1304: Output Parameter:
1305: . A - the new matrix
1307: Notes:
1308: The user is responsible of destroying the returned matrix.
1310: The matrix contains all columns of the BV, not just the active columns.
1312: Level: intermediate
1314: .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1315: @*/
1316: PetscErrorCode BVCreateMat(BV bv,Mat *A)1317: {
1318: PetscScalar *aa;
1319: const PetscScalar *vv;
1322: BVCheckSizes(bv,1);
1325: MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,bv->N,bv->m,NULL,A);
1326: MatDenseGetArrayWrite(*A,&aa);
1327: BVGetArrayRead(bv,&vv);
1328: PetscArraycpy(aa,vv,bv->m*bv->n);
1329: BVRestoreArrayRead(bv,&vv);
1330: MatDenseRestoreArrayWrite(*A,&aa);
1331: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1332: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1333: PetscFunctionReturn(0);
1334: }
1336: /*@
1337: BVGetMat - Returns a Mat object of dense type that shares the memory of
1338: the BV object.
1340: Collective on bv
1342: Input Parameter:
1343: . bv - the basis vectors context
1345: Output Parameter:
1346: . A - the matrix
1348: Notes:
1349: The returned matrix contains only the active columns. If the content of
1350: the Mat is modified, these changes are also done in the BV object. The
1351: user must call BVRestoreMat() when no longer needed.
1353: This operation implies a call to BVGetArray(), which may result in data
1354: copies.
1356: Level: advanced
1358: .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1359: @*/
1360: PetscErrorCode BVGetMat(BV bv,Mat *A)1361: {
1362: PetscScalar *vv,*aa;
1363: PetscBool create=PETSC_FALSE;
1364: PetscInt m,cols;
1367: BVCheckSizes(bv,1);
1369: if (bv->ops->getmat) (*bv->ops->getmat)(bv,A);
1370: else {
1371: m = bv->k-bv->l;
1372: if (!bv->Aget) create=PETSC_TRUE;
1373: else {
1374: MatDenseGetArray(bv->Aget,&aa);
1376: MatGetSize(bv->Aget,NULL,&cols);
1377: if (cols!=m) {
1378: MatDestroy(&bv->Aget);
1379: create=PETSC_TRUE;
1380: }
1381: }
1382: BVGetArray(bv,&vv);
1383: if (create) {
1384: MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,vv,&bv->Aget); /* pass a pointer to avoid allocation of storage */
1385: MatDenseReplaceArray(bv->Aget,NULL); /* replace with a null pointer, the value after BVRestoreMat */
1386: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Aget);
1387: }
1388: MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->n); /* set the actual pointer */
1389: *A = bv->Aget;
1390: }
1391: PetscFunctionReturn(0);
1392: }
1394: /*@
1395: BVRestoreMat - Restores the Mat obtained with BVGetMat().
1397: Logically Collective on bv
1399: Input Parameters:
1400: + bv - the basis vectors context
1401: - A - the fetched matrix
1403: Note:
1404: A call to this function must match a previous call of BVGetMat().
1405: The effect is that the contents of the Mat are copied back to the
1406: BV internal data structures.
1408: Level: advanced
1410: .seealso: BVGetMat(), BVRestoreArray()
1411: @*/
1412: PetscErrorCode BVRestoreMat(BV bv,Mat *A)1413: {
1414: PetscScalar *vv,*aa;
1417: BVCheckSizes(bv,1);
1421: if (bv->ops->restoremat) (*bv->ops->restoremat)(bv,A);
1422: else {
1423: MatDenseGetArray(bv->Aget,&aa);
1424: vv = aa-(bv->nc+bv->l)*bv->n;
1425: MatDenseResetArray(bv->Aget);
1426: BVRestoreArray(bv,&vv);
1427: *A = NULL;
1428: }
1429: PetscFunctionReturn(0);
1430: }
1432: /*
1433: Copy all user-provided attributes of V to another BV object W
1434: */
1435: static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)1436: {
1437: BVSetType(W,((PetscObject)V)->type_name);
1438: W->orthog_type = V->orthog_type;
1439: W->orthog_ref = V->orthog_ref;
1440: W->orthog_eta = V->orthog_eta;
1441: W->orthog_block = V->orthog_block;
1442: if (V->matrix) PetscObjectReference((PetscObject)V->matrix);
1443: W->matrix = V->matrix;
1444: W->indef = V->indef;
1445: W->vmm = V->vmm;
1446: W->rrandom = V->rrandom;
1447: W->deftol = V->deftol;
1448: if (V->rand) PetscObjectReference((PetscObject)V->rand);
1449: W->rand = V->rand;
1450: W->sfocalled = V->sfocalled;
1451: if (V->ops->duplicate) (*V->ops->duplicate)(V,W);
1452: PetscObjectStateIncrease((PetscObject)W);
1453: PetscFunctionReturn(0);
1454: }
1456: /*@
1457: BVDuplicate - Creates a new basis vector object of the same type and
1458: dimensions as an existing one.
1460: Collective on V
1462: Input Parameter:
1463: . V - basis vectors context
1465: Output Parameter:
1466: . W - location to put the new BV1468: Notes:
1469: The new BV has the same type and dimensions as V, and it shares the same
1470: template vector. Also, the inner product matrix and orthogonalization
1471: options are copied.
1473: BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1474: for the new basis vectors. Use BVCopy() to copy the contents.
1476: Level: intermediate
1478: .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1479: @*/
1480: PetscErrorCode BVDuplicate(BV V,BV *W)1481: {
1484: BVCheckSizes(V,1);
1486: BVCreate(PetscObjectComm((PetscObject)V),W);
1487: BVSetSizesFromVec(*W,V->t,V->m);
1488: BVDuplicate_Private(V,*W);
1489: PetscFunctionReturn(0);
1490: }
1492: /*@
1493: BVDuplicateResize - Creates a new basis vector object of the same type and
1494: dimensions as an existing one, but with possibly different number of columns.
1496: Collective on V
1498: Input Parameters:
1499: + V - basis vectors context
1500: - m - the new number of columns
1502: Output Parameter:
1503: . W - location to put the new BV1505: Note:
1506: This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1507: contents of V are not copied to W.
1509: Level: intermediate
1511: .seealso: BVDuplicate(), BVResize()
1512: @*/
1513: PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)1514: {
1517: BVCheckSizes(V,1);
1520: BVCreate(PetscObjectComm((PetscObject)V),W);
1521: BVSetSizesFromVec(*W,V->t,m);
1522: BVDuplicate_Private(V,*W);
1523: PetscFunctionReturn(0);
1524: }
1526: /*@
1527: BVGetCachedBV - Returns a BV object stored internally that holds the
1528: result of B*X after a call to BVApplyMatrixBV().
1530: Not collective
1532: Input Parameter:
1533: . bv - the basis vectors context
1535: Output Parameter:
1536: . cached - the cached BV1538: Note:
1539: The cached BV is created if not available previously.
1541: Level: developer
1543: .seealso: BVApplyMatrixBV()
1544: @*/
1545: PetscErrorCode BVGetCachedBV(BV bv,BV *cached)1546: {
1549: BVCheckSizes(bv,1);
1550: if (!bv->cached) {
1551: BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached);
1552: BVSetSizesFromVec(bv->cached,bv->t,bv->m);
1553: BVDuplicate_Private(bv,bv->cached);
1554: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->cached);
1555: }
1556: *cached = bv->cached;
1557: PetscFunctionReturn(0);
1558: }
1560: /*@
1561: BVCopy - Copies a basis vector object into another one, W <- V.
1563: Logically Collective on V
1565: Input Parameter:
1566: . V - basis vectors context
1568: Output Parameter:
1569: . W - the copy
1571: Note:
1572: Both V and W must be distributed in the same manner; local copies are
1573: done. Only active columns (excluding the leading ones) are copied.
1574: In the destination W, columns are overwritten starting from the leading ones.
1575: Constraints are not copied.
1577: Level: beginner
1579: .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1580: @*/
1581: PetscErrorCode BVCopy(BV V,BV W)1582: {
1583: PetscScalar *womega;
1584: const PetscScalar *vomega;
1588: BVCheckSizes(V,1);
1589: BVCheckOp(V,1,copy);
1592: BVCheckSizes(W,2);
1596: if (V==W || !V->n) PetscFunctionReturn(0);
1598: PetscLogEventBegin(BV_Copy,V,W,0,0);
1599: if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1600: /* copy signature */
1601: BV_AllocateSignature(W);
1602: VecGetArrayRead(V->omega,&vomega);
1603: VecGetArray(W->omega,&womega);
1604: PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l);
1605: VecRestoreArray(W->omega,&womega);
1606: VecRestoreArrayRead(V->omega,&vomega);
1607: }
1608: (*V->ops->copy)(V,W);
1609: PetscLogEventEnd(BV_Copy,V,W,0,0);
1610: PetscObjectStateIncrease((PetscObject)W);
1611: PetscFunctionReturn(0);
1612: }
1614: /*@
1615: BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1617: Logically Collective on V
1619: Input Parameters:
1620: + V - basis vectors context
1621: - j - the column number to be copied
1623: Output Parameter:
1624: . w - the copied column
1626: Note:
1627: Both V and w must be distributed in the same manner; local copies are done.
1629: Level: beginner
1631: .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1632: @*/
1633: PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)1634: {
1635: PetscInt n,N;
1636: Vec z;
1640: BVCheckSizes(V,1);
1645: VecGetSize(w,&N);
1646: VecGetLocalSize(w,&n);
1649: PetscLogEventBegin(BV_Copy,V,w,0,0);
1650: BVGetColumn(V,j,&z);
1651: VecCopy(z,w);
1652: BVRestoreColumn(V,j,&z);
1653: PetscLogEventEnd(BV_Copy,V,w,0,0);
1654: PetscFunctionReturn(0);
1655: }
1657: /*@
1658: BVCopyColumn - Copies the values from one of the columns to another one.
1660: Logically Collective on V
1662: Input Parameters:
1663: + V - basis vectors context
1664: . j - the number of the source column
1665: - i - the number of the destination column
1667: Level: beginner
1669: .seealso: BVCopy(), BVCopyVec()
1670: @*/
1671: PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)1672: {
1673: Vec z,w;
1674: PetscScalar *omega;
1678: BVCheckSizes(V,1);
1681: if (j==i) PetscFunctionReturn(0);
1683: PetscLogEventBegin(BV_Copy,V,0,0,0);
1684: if (V->omega) {
1685: VecGetArray(V->omega,&omega);
1686: omega[i] = omega[j];
1687: VecRestoreArray(V->omega,&omega);
1688: }
1689: if (V->ops->copycolumn) (*V->ops->copycolumn)(V,j,i);
1690: else {
1691: BVGetColumn(V,j,&z);
1692: BVGetColumn(V,i,&w);
1693: VecCopy(z,w);
1694: BVRestoreColumn(V,j,&z);
1695: BVRestoreColumn(V,i,&w);
1696: }
1697: PetscLogEventEnd(BV_Copy,V,0,0,0);
1698: PetscObjectStateIncrease((PetscObject)V);
1699: PetscFunctionReturn(0);
1700: }
1702: static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)1703: {
1704: PetscInt ncols;
1706: ncols = left? bv->nc+bv->l: bv->m-bv->l;
1707: if (*split && ncols!=(*split)->m) BVDestroy(split);
1708: if (!*split) {
1709: BVCreate(PetscObjectComm((PetscObject)bv),split);
1710: PetscLogObjectParent((PetscObject)bv,(PetscObject)*split);
1711: (*split)->issplit = left? 1: 2;
1712: (*split)->splitparent = bv;
1713: BVSetSizesFromVec(*split,bv->t,ncols);
1714: BVDuplicate_Private(bv,*split);
1715: }
1716: (*split)->l = 0;
1717: (*split)->k = left? bv->l: bv->k-bv->l;
1718: (*split)->nc = left? bv->nc: 0;
1719: (*split)->m = ncols-(*split)->nc;
1720: if ((*split)->nc) {
1721: (*split)->ci[0] = -(*split)->nc-1;
1722: (*split)->ci[1] = -(*split)->nc-1;
1723: }
1724: if (left) PetscObjectStateGet((PetscObject)*split,&bv->lstate);
1725: else PetscObjectStateGet((PetscObject)*split,&bv->rstate);
1726: PetscFunctionReturn(0);
1727: }
1729: /*@
1730: BVGetSplit - Splits the BV object into two BV objects that share the
1731: internal data, one of them containing the leading columns and the other
1732: one containing the remaining columns.
1734: Logically Collective on bv
1736: Input Parameter:
1737: . bv - the basis vectors context
1739: Output Parameters:
1740: + L - left BV containing leading columns (can be NULL)
1741: - R - right BV containing remaining columns (can be NULL)
1743: Notes:
1744: The columns are split in two sets. The leading columns (including the
1745: constraints) are assigned to the left BV and the remaining columns
1746: are assigned to the right BV. The number of leading columns, as
1747: specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1748: guarantee that both L and R have at least one column).
1750: The returned BV's must be seen as references (not copies) of the input
1751: BV, that is, modifying them will change the entries of bv as well.
1752: The returned BV's must not be destroyed. BVRestoreSplit() must be called
1753: when they are no longer needed.
1755: Pass NULL for any of the output BV's that is not needed.
1757: Level: advanced
1759: .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints()
1760: @*/
1761: PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)1762: {
1765: BVCheckSizes(bv,1);
1768: bv->lsplit = bv->nc+bv->l;
1769: BVGetSplit_Private(bv,PETSC_TRUE,&bv->L);
1770: BVGetSplit_Private(bv,PETSC_FALSE,&bv->R);
1771: if (L) *L = bv->L;
1772: if (R) *R = bv->R;
1773: PetscFunctionReturn(0);
1774: }
1776: /*@
1777: BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1779: Logically Collective on bv
1781: Input Parameters:
1782: + bv - the basis vectors context
1783: . L - left BV obtained with BVGetSplit()
1784: - R - right BV obtained with BVGetSplit()
1786: Note:
1787: The arguments must match the corresponding call to BVGetSplit().
1789: Level: advanced
1791: .seealso: BVGetSplit()
1792: @*/
1793: PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)1794: {
1797: BVCheckSizes(bv,1);
1804: if (bv->ops->restoresplit) (*bv->ops->restoresplit)(bv,L,R);
1805: bv->lsplit = 0;
1806: if (L) *L = NULL;
1807: if (R) *R = NULL;
1808: PetscFunctionReturn(0);
1809: }
1811: /*@
1812: BVSetDefiniteTolerance - Set the tolerance to be used when checking a
1813: definite inner product.
1815: Logically Collective on bv
1817: Input Parameters:
1818: + bv - basis vectors
1819: - deftol - the tolerance
1821: Options Database Key:
1822: . -bv_definite_tol <deftol> - the tolerance
1824: Notes:
1825: When using a non-standard inner product, see BVSetMatrix(), the solver needs
1826: to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
1827: been declared indefinite, the value z'*B*z must be positive, but due to
1828: rounding error a tiny value may become negative. A tolerance is used to
1829: detect this situation. Likewise, in complex arithmetic z'*B*z should be
1830: real, and we use the same tolerance to check whether a nonzero imaginary part
1831: can be considered negligible.
1833: This function sets this tolerance, which defaults to 10*eps, where eps is
1834: the machine epsilon. The default value should be good for most applications.
1836: Level: advanced
1838: .seealso: BVSetMatrix()
1839: @*/
1840: PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)1841: {
1844: if (deftol == PETSC_DEFAULT) bv->deftol = 10*PETSC_MACHINE_EPSILON;
1845: else {
1847: bv->deftol = deftol;
1848: }
1849: PetscFunctionReturn(0);
1850: }
1852: /*@
1853: BVGetDefiniteTolerance - Returns the tolerance for checking a definite
1854: inner product.
1856: Not Collective
1858: Input Parameter:
1859: . bv - the basis vectors
1861: Output Parameters:
1862: . deftol - the tolerance
1864: Level: advanced
1866: .seealso: BVSetDefiniteTolerance()
1867: @*/
1868: PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)1869: {
1872: *deftol = bv->deftol;
1873: PetscFunctionReturn(0);
1874: }