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: BV operations involving global communication
12: */
14: #include <slepc/private/bvimpl.h> 16: /*
17: BVDot for the particular case of non-standard inner product with
18: matrix B, which is assumed to be symmetric (or complex Hermitian)
19: */
20: static inline PetscErrorCode BVDot_Private(BV X,BV Y,Mat M) 21: {
22: PetscObjectId idx,idy;
23: PetscInt i,j,jend,m;
24: PetscScalar *marray;
25: PetscBool symm=PETSC_FALSE;
26: Mat B;
27: Vec z;
29: BVCheckOp(Y,1,dotvec);
30: MatGetSize(M,&m,NULL);
31: MatDenseGetArray(M,&marray);
32: PetscObjectGetId((PetscObject)X,&idx);
33: PetscObjectGetId((PetscObject)Y,&idy);
34: B = Y->matrix;
35: Y->matrix = NULL;
36: if (idx==idy) symm=PETSC_TRUE; /* M=X'BX is symmetric */
37: jend = X->k;
38: for (j=X->l;j<jend;j++) {
39: if (symm) Y->k = j+1;
40: BVGetColumn(X->cached,j,&z);
41: (*Y->ops->dotvec)(Y,z,marray+j*m+Y->l);
42: BVRestoreColumn(X->cached,j,&z);
43: if (symm) {
44: for (i=X->l;i<j;i++)
45: marray[j+i*m] = PetscConj(marray[i+j*m]);
46: }
47: }
48: MatDenseRestoreArray(M,&marray);
49: Y->matrix = B;
50: PetscFunctionReturn(0);
51: }
53: /*@
54: BVDot - Computes the 'block-dot' product of two basis vectors objects.
56: Collective on X
58: Input Parameters:
59: + X - first basis vectors
60: - Y - second basis vectors
62: Output Parameter:
63: . M - the resulting matrix
65: Notes:
66: This is the generalization of VecDot() for a collection of vectors, M = Y^H*X.
67: The result is a matrix M whose entry m_ij is equal to y_i^H x_j (where y_i^H
68: denotes the conjugate transpose of y_i).
70: If a non-standard inner product has been specified with BVSetMatrix(),
71: then the result is M = Y^H*B*X. In this case, both X and Y must have
72: the same associated matrix.
74: On entry, M must be a sequential dense Mat with dimensions m,n at least, where
75: m is the number of active columns of Y and n is the number of active columns of X.
76: Only rows (resp. columns) of M starting from ly (resp. lx) are computed,
77: where ly (resp. lx) is the number of leading columns of Y (resp. X).
79: X and Y need not be different objects.
81: Level: intermediate
83: .seealso: BVDotVec(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
84: @*/
85: PetscErrorCode BVDot(BV X,BV Y,Mat M) 86: {
87: PetscInt m,n;
93: BVCheckSizes(X,1);
95: BVCheckSizes(Y,2);
100: MatGetSize(M,&m,&n);
105: if (X->l==X->k || Y->l==Y->k) PetscFunctionReturn(0);
107: PetscLogEventBegin(BV_Dot,X,Y,0,0);
108: if (X->matrix) { /* non-standard inner product */
109: /* compute BX first */
110: BV_IPMatMultBV(X);
111: if (X->vmm==BV_MATMULT_VECS) {
112: /* perform computation column by column */
113: BVDot_Private(X,Y,M);
114: } else (*X->ops->dot)(X->cached,Y,M);
115: } else (*X->ops->dot)(X,Y,M);
116: PetscLogEventEnd(BV_Dot,X,Y,0,0);
117: PetscFunctionReturn(0);
118: }
120: /*@
121: BVDotVec - Computes multiple dot products of a vector against all the
122: column vectors of a BV.
124: Collective on X
126: Input Parameters:
127: + X - basis vectors
128: - y - a vector
130: Output Parameter:
131: . m - an array where the result must be placed
133: Notes:
134: This is analogue to VecMDot(), but using BV to represent a collection
135: of vectors. The result is m = X^H*y, so m_i is equal to x_j^H y. Note
136: that here X is transposed as opposed to BVDot().
138: If a non-standard inner product has been specified with BVSetMatrix(),
139: then the result is m = X^H*B*y.
141: The length of array m must be equal to the number of active columns of X
142: minus the number of leading columns, i.e. the first entry of m is the
143: product of the first non-leading column with y.
145: Level: intermediate
147: .seealso: BVDot(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
148: @*/
149: PetscErrorCode BVDotVec(BV X,Vec y,PetscScalar m[])150: {
151: PetscInt n;
156: BVCheckSizes(X,1);
157: BVCheckOp(X,1,dotvec);
161: VecGetLocalSize(y,&n);
164: PetscLogEventBegin(BV_DotVec,X,y,0,0);
165: (*X->ops->dotvec)(X,y,m);
166: PetscLogEventEnd(BV_DotVec,X,y,0,0);
167: PetscFunctionReturn(0);
168: }
170: /*@
171: BVDotVecBegin - Starts a split phase dot product computation.
173: Input Parameters:
174: + X - basis vectors
175: . y - a vector
176: - m - an array where the result will go (can be NULL)
178: Note:
179: Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
181: Level: advanced
183: .seealso: BVDotVecEnd(), BVDotVec()
184: @*/
185: PetscErrorCode BVDotVecBegin(BV X,Vec y,PetscScalar *m)186: {
187: PetscInt i,n,nv;
188: PetscSplitReduction *sr;
189: MPI_Comm comm;
194: BVCheckSizes(X,1);
198: VecGetLocalSize(y,&n);
201: if (X->ops->dotvec_begin) (*X->ops->dotvec_begin)(X,y,m);
202: else {
203: BVCheckOp(X,1,dotvec_local);
204: nv = X->k-X->l;
205: PetscObjectGetComm((PetscObject)X,&comm);
206: PetscSplitReductionGet(comm,&sr);
208: for (i=0;i<nv;i++) {
209: if (sr->numopsbegin+i >= sr->maxops) PetscSplitReductionExtend(sr);
210: sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
211: sr->invecs[sr->numopsbegin+i] = (void*)X;
212: }
213: PetscLogEventBegin(BV_DotVec,X,y,0,0);
214: (*X->ops->dotvec_local)(X,y,sr->lvalues+sr->numopsbegin);
215: sr->numopsbegin += nv;
216: PetscLogEventEnd(BV_DotVec,X,y,0,0);
217: }
218: PetscFunctionReturn(0);
219: }
221: /*@
222: BVDotVecEnd - Ends a split phase dot product computation.
224: Input Parameters:
225: + X - basis vectors
226: . y - a vector
227: - m - an array where the result will go
229: Note:
230: Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
232: Level: advanced
234: .seealso: BVDotVecBegin(), BVDotVec()
235: @*/
236: PetscErrorCode BVDotVecEnd(BV X,Vec y,PetscScalar *m)237: {
238: PetscInt i,nv;
239: PetscSplitReduction *sr;
240: MPI_Comm comm;
244: BVCheckSizes(X,1);
246: if (X->ops->dotvec_end) (*X->ops->dotvec_end)(X,y,m);
247: else {
248: nv = X->k-X->l;
249: PetscObjectGetComm((PetscObject)X,&comm);
250: PetscSplitReductionGet(comm,&sr);
251: PetscSplitReductionEnd(sr);
256: for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
258: /* Finished getting all the results so reset to no outstanding requests */
259: if (sr->numopsend == sr->numopsbegin) {
260: sr->state = STATE_BEGIN;
261: sr->numopsend = 0;
262: sr->numopsbegin = 0;
263: }
264: }
265: PetscFunctionReturn(0);
266: }
268: /*@
269: BVDotColumn - Computes multiple dot products of a column against all the
270: previous columns of a BV.
272: Collective on X
274: Input Parameters:
275: + X - basis vectors
276: - j - the column index
278: Output Parameter:
279: . q - an array where the result must be placed
281: Notes:
282: This operation is equivalent to BVDotVec() but it uses column j of X
283: rather than taking a Vec as an argument. The number of active columns of
284: X is set to j before the computation, and restored afterwards.
285: If X has leading columns specified, then these columns do not participate
286: in the computation. Therefore, the length of array q must be equal to j
287: minus the number of leading columns.
289: Developer Notes:
290: If q is NULL, then the result is written in position nc+l of the internal
291: buffer vector, see BVGetBufferVec().
293: Level: advanced
295: .seealso: BVDot(), BVDotVec(), BVSetActiveColumns(), BVSetMatrix()
296: @*/
297: PetscErrorCode BVDotColumn(BV X,PetscInt j,PetscScalar *q)298: {
299: PetscInt ksave;
300: Vec y;
305: BVCheckSizes(X,1);
306: BVCheckOp(X,1,dotvec);
311: PetscLogEventBegin(BV_DotVec,X,0,0,0);
312: ksave = X->k;
313: X->k = j;
314: if (!q && !X->buffer) BVGetBufferVec(X,&X->buffer);
315: BVGetColumn(X,j,&y);
316: (*X->ops->dotvec)(X,y,q);
317: BVRestoreColumn(X,j,&y);
318: X->k = ksave;
319: PetscLogEventEnd(BV_DotVec,X,0,0,0);
320: PetscFunctionReturn(0);
321: }
323: /*@
324: BVDotColumnBegin - Starts a split phase dot product computation.
326: Input Parameters:
327: + X - basis vectors
328: . j - the column index
329: - m - an array where the result will go (can be NULL)
331: Note:
332: Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
334: Level: advanced
336: .seealso: BVDotColumnEnd(), BVDotColumn()
337: @*/
338: PetscErrorCode BVDotColumnBegin(BV X,PetscInt j,PetscScalar *m)339: {
340: PetscInt i,nv,ksave;
341: PetscSplitReduction *sr;
342: MPI_Comm comm;
343: Vec y;
348: BVCheckSizes(X,1);
352: ksave = X->k;
353: X->k = j;
354: BVGetColumn(X,j,&y);
356: if (X->ops->dotvec_begin) (*X->ops->dotvec_begin)(X,y,m);
357: else {
358: BVCheckOp(X,1,dotvec_local);
359: nv = X->k-X->l;
360: PetscObjectGetComm((PetscObject)X,&comm);
361: PetscSplitReductionGet(comm,&sr);
363: for (i=0;i<nv;i++) {
364: if (sr->numopsbegin+i >= sr->maxops) PetscSplitReductionExtend(sr);
365: sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
366: sr->invecs[sr->numopsbegin+i] = (void*)X;
367: }
368: PetscLogEventBegin(BV_DotVec,X,0,0,0);
369: (*X->ops->dotvec_local)(X,y,sr->lvalues+sr->numopsbegin);
370: sr->numopsbegin += nv;
371: PetscLogEventEnd(BV_DotVec,X,0,0,0);
372: }
373: BVRestoreColumn(X,j,&y);
374: X->k = ksave;
375: PetscFunctionReturn(0);
376: }
378: /*@
379: BVDotColumnEnd - Ends a split phase dot product computation.
381: Input Parameters:
382: + X - basis vectors
383: . j - the column index
384: - m - an array where the result will go
386: Notes:
387: Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
389: Level: advanced
391: .seealso: BVDotColumnBegin(), BVDotColumn()
392: @*/
393: PetscErrorCode BVDotColumnEnd(BV X,PetscInt j,PetscScalar *m)394: {
395: PetscInt i,nv,ksave;
396: PetscSplitReduction *sr;
397: MPI_Comm comm;
398: Vec y;
403: BVCheckSizes(X,1);
407: ksave = X->k;
408: X->k = j;
410: if (X->ops->dotvec_end) {
411: BVGetColumn(X,j,&y);
412: (*X->ops->dotvec_end)(X,y,m);
413: BVRestoreColumn(X,j,&y);
414: } else {
415: nv = X->k-X->l;
416: PetscObjectGetComm((PetscObject)X,&comm);
417: PetscSplitReductionGet(comm,&sr);
418: PetscSplitReductionEnd(sr);
423: for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
425: /* Finished getting all the results so reset to no outstanding requests */
426: if (sr->numopsend == sr->numopsbegin) {
427: sr->state = STATE_BEGIN;
428: sr->numopsend = 0;
429: sr->numopsbegin = 0;
430: }
431: }
432: X->k = ksave;
433: PetscFunctionReturn(0);
434: }
436: static inline PetscErrorCode BVNorm_Private(BV bv,Vec z,NormType type,PetscReal *val)437: {
438: PetscScalar p;
440: BV_IPMatMult(bv,z);
441: VecDot(bv->Bx,z,&p);
442: BV_SafeSqrt(bv,p,val);
443: PetscFunctionReturn(0);
444: }
446: static inline PetscErrorCode BVNorm_Begin_Private(BV bv,Vec z,NormType type,PetscReal *val)447: {
448: PetscScalar p;
450: BV_IPMatMult(bv,z);
451: VecDotBegin(bv->Bx,z,&p);
452: PetscFunctionReturn(0);
453: }
455: static inline PetscErrorCode BVNorm_End_Private(BV bv,Vec z,NormType type,PetscReal *val)456: {
457: PetscScalar p;
459: VecDotEnd(bv->Bx,z,&p);
460: BV_SafeSqrt(bv,p,val);
461: PetscFunctionReturn(0);
462: }
464: /*@
465: BVNorm - Computes the matrix norm of the BV.
467: Collective on bv
469: Input Parameters:
470: + bv - basis vectors
471: - type - the norm type
473: Output Parameter:
474: . val - the norm
476: Notes:
477: All active columns (except the leading ones) are considered as a matrix.
478: The allowed norms are NORM_1, NORM_FROBENIUS, and NORM_INFINITY.
480: This operation fails if a non-standard inner product has been
481: specified with BVSetMatrix().
483: Level: intermediate
485: .seealso: BVNormVec(), BVNormColumn(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
486: @*/
487: PetscErrorCode BVNorm(BV bv,NormType type,PetscReal *val)488: {
493: BVCheckSizes(bv,1);
498: PetscLogEventBegin(BV_Norm,bv,0,0,0);
499: (*bv->ops->norm)(bv,-1,type,val);
500: PetscLogEventEnd(BV_Norm,bv,0,0,0);
501: PetscFunctionReturn(0);
502: }
504: /*@
505: BVNormVec - Computes the norm of a given vector.
507: Collective on bv
509: Input Parameters:
510: + bv - basis vectors
511: . v - the vector
512: - type - the norm type
514: Output Parameter:
515: . val - the norm
517: Notes:
518: This is the analogue of BVNormColumn() but for a vector that is not in the BV.
519: If a non-standard inner product has been specified with BVSetMatrix(),
520: then the returned value is sqrt(v'*B*v), where B is the inner product
521: matrix (argument 'type' is ignored). Otherwise, VecNorm() is called.
523: Level: developer
525: .seealso: BVNorm(), BVNormColumn(), BVSetMatrix()
526: @*/
527: PetscErrorCode BVNormVec(BV bv,Vec v,NormType type,PetscReal *val)528: {
529: PetscInt n;
536: BVCheckSizes(bv,1);
542: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
543: if (bv->matrix) { /* non-standard inner product */
544: VecGetLocalSize(v,&n);
546: BVNorm_Private(bv,v,type,val);
547: } else VecNorm(v,type,val);
548: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
549: PetscFunctionReturn(0);
550: }
552: /*@
553: BVNormVecBegin - Starts a split phase norm computation.
555: Input Parameters:
556: + bv - basis vectors
557: . v - the vector
558: . type - the norm type
559: - val - the norm
561: Note:
562: Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
564: Level: advanced
566: .seealso: BVNormVecEnd(), BVNormVec()
567: @*/
568: PetscErrorCode BVNormVecBegin(BV bv,Vec v,NormType type,PetscReal *val)569: {
570: PetscInt n;
577: BVCheckSizes(bv,1);
583: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
584: if (bv->matrix) { /* non-standard inner product */
585: VecGetLocalSize(v,&n);
587: BVNorm_Begin_Private(bv,v,type,val);
588: } else VecNormBegin(v,type,val);
589: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
590: PetscFunctionReturn(0);
591: }
593: /*@
594: BVNormVecEnd - Ends a split phase norm computation.
596: Input Parameters:
597: + bv - basis vectors
598: . v - the vector
599: . type - the norm type
600: - val - the norm
602: Note:
603: Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
605: Level: advanced
607: .seealso: BVNormVecBegin(), BVNormVec()
608: @*/
609: PetscErrorCode BVNormVecEnd(BV bv,Vec v,NormType type,PetscReal *val)610: {
615: BVCheckSizes(bv,1);
619: if (bv->matrix) BVNorm_End_Private(bv,v,type,val); /* non-standard inner product */
620: else VecNormEnd(v,type,val);
621: PetscFunctionReturn(0);
622: }
624: /*@
625: BVNormColumn - Computes the vector norm of a selected column.
627: Collective on bv
629: Input Parameters:
630: + bv - basis vectors
631: . j - column number to be used
632: - type - the norm type
634: Output Parameter:
635: . val - the norm
637: Notes:
638: The norm of V[j] is computed (NORM_1, NORM_2, or NORM_INFINITY).
639: If a non-standard inner product has been specified with BVSetMatrix(),
640: then the returned value is sqrt(V[j]'*B*V[j]),
641: where B is the inner product matrix (argument 'type' is ignored).
643: Level: intermediate
645: .seealso: BVNorm(), BVNormVec(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
646: @*/
647: PetscErrorCode BVNormColumn(BV bv,PetscInt j,NormType type,PetscReal *val)648: {
649: Vec z;
656: BVCheckSizes(bv,1);
661: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
662: if (bv->matrix) { /* non-standard inner product */
663: BVGetColumn(bv,j,&z);
664: BVNorm_Private(bv,z,type,val);
665: BVRestoreColumn(bv,j,&z);
666: } else (*bv->ops->norm)(bv,j,type,val);
667: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
668: PetscFunctionReturn(0);
669: }
671: /*@
672: BVNormColumnBegin - Starts a split phase norm computation.
674: Input Parameters:
675: + bv - basis vectors
676: . j - column number to be used
677: . type - the norm type
678: - val - the norm
680: Note:
681: Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
683: Level: advanced
685: .seealso: BVNormColumnEnd(), BVNormColumn()
686: @*/
687: PetscErrorCode BVNormColumnBegin(BV bv,PetscInt j,NormType type,PetscReal *val)688: {
689: PetscSplitReduction *sr;
690: PetscReal lresult;
691: MPI_Comm comm;
692: Vec z;
699: BVCheckSizes(bv,1);
704: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
705: BVGetColumn(bv,j,&z);
706: if (bv->matrix) BVNorm_Begin_Private(bv,z,type,val); /* non-standard inner product */
707: else if (bv->ops->norm_begin) (*bv->ops->norm_begin)(bv,j,type,val);
708: else {
709: BVCheckOp(bv,1,norm_local);
710: PetscObjectGetComm((PetscObject)z,&comm);
711: PetscSplitReductionGet(comm,&sr);
713: if (sr->numopsbegin >= sr->maxops) PetscSplitReductionExtend(sr);
714: sr->invecs[sr->numopsbegin] = (void*)bv;
715: (*bv->ops->norm_local)(bv,j,type,&lresult);
716: if (type == NORM_2) lresult = lresult*lresult;
717: if (type == NORM_MAX) sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_MAX;
718: else sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_SUM;
719: sr->lvalues[sr->numopsbegin++] = lresult;
720: }
721: BVRestoreColumn(bv,j,&z);
722: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
723: PetscFunctionReturn(0);
724: }
726: /*@
727: BVNormColumnEnd - Ends a split phase norm computation.
729: Input Parameters:
730: + bv - basis vectors
731: . j - column number to be used
732: . type - the norm type
733: - val - the norm
735: Note:
736: Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
738: Level: advanced
740: .seealso: BVNormColumnBegin(), BVNormColumn()
741: @*/
742: PetscErrorCode BVNormColumnEnd(BV bv,PetscInt j,NormType type,PetscReal *val)743: {
744: PetscSplitReduction *sr;
745: MPI_Comm comm;
746: Vec z;
753: BVCheckSizes(bv,1);
757: BVGetColumn(bv,j,&z);
758: if (bv->matrix) BVNorm_End_Private(bv,z,type,val); /* non-standard inner product */
759: else if (bv->ops->norm_end) (*bv->ops->norm_end)(bv,j,type,val);
760: else {
761: PetscObjectGetComm((PetscObject)z,&comm);
762: PetscSplitReductionGet(comm,&sr);
763: PetscSplitReductionEnd(sr);
768: *val = PetscRealPart(sr->gvalues[sr->numopsend++]);
769: if (type == NORM_2) *val = PetscSqrtReal(*val);
770: if (sr->numopsend == sr->numopsbegin) {
771: sr->state = STATE_BEGIN;
772: sr->numopsend = 0;
773: sr->numopsbegin = 0;
774: }
775: }
776: BVRestoreColumn(bv,j,&z);
777: PetscFunctionReturn(0);
778: }
780: /*@
781: BVNormalize - Normalize all columns (starting from the leading ones).
783: Collective on bv
785: Input Parameters:
786: + bv - basis vectors
787: - eigi - (optional) imaginary parts of eigenvalues
789: Notes:
790: On output, all columns will have unit norm. The normalization is done with
791: respect to the 2-norm, or to the B-norm if a non-standard inner product has
792: been specified with BVSetMatrix(), see BVNormColumn().
794: If the optional argument eigi is passed (taken into account only in real
795: scalars) it is interpreted as the imaginary parts of the eigenvalues and
796: the BV is supposed to contain the corresponding eigenvectors. Suppose the
797: first three values are eigi = { 0, alpha, -alpha }, then the first column
798: is normalized as usual, but the second and third ones are normalized assuming
799: that they contain the real and imaginary parts of a complex conjugate pair of
800: eigenvectors.
802: If eigi is passed, the inner-product matrix is ignored.
804: If there are leading columns, they are not modified (are assumed to be already
805: normalized).
807: Level: intermediate
809: .seealso: BVNormColumn()
810: @*/
811: PetscErrorCode BVNormalize(BV bv,PetscScalar *eigi)812: {
813: PetscReal norm;
814: PetscInt i;
815: Vec v;
816: #if !defined(PETSC_USE_COMPLEX)
817: Vec v1;
818: #endif
822: BVCheckSizes(bv,1);
824: PetscLogEventBegin(BV_Normalize,bv,0,0,0);
825: if (bv->matrix && !eigi) {
826: for (i=bv->l;i<bv->k;i++) {
827: BVNormColumn(bv,i,NORM_2,&norm);
828: BVScaleColumn(bv,i,1.0/norm);
829: }
830: } else if (bv->ops->normalize) (*bv->ops->normalize)(bv,eigi);
831: else {
832: for (i=bv->l;i<bv->k;i++) {
833: #if !defined(PETSC_USE_COMPLEX)
834: if (eigi && eigi[i] != 0.0) {
835: BVGetColumn(bv,i,&v);
836: BVGetColumn(bv,i+1,&v1);
837: VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
838: BVRestoreColumn(bv,i,&v);
839: BVRestoreColumn(bv,i+1,&v1);
840: i++;
841: } else
842: #endif
843: {
844: BVGetColumn(bv,i,&v);
845: VecNormalize(v,NULL);
846: BVRestoreColumn(bv,i,&v);
847: }
848: }
849: }
850: PetscLogEventEnd(BV_Normalize,bv,0,0,0);
851: PetscObjectStateIncrease((PetscObject)bv);
852: PetscFunctionReturn(0);
853: }
855: /*
856: Compute Y^H*A*X: right part column by column (with MatMult) and bottom
857: part row by row (with MatMultHermitianTranspose); result placed in marray[*,ldm]
858: */
859: static inline PetscErrorCode BVMatProject_Vec(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)860: {
861: PetscInt i,j,lx,ly,kx,ky,ulim;
862: Vec z,f;
864: lx = X->l; kx = X->k;
865: ly = Y->l; ky = Y->k;
866: BVCreateVec(X,&f);
867: BVCheckOp(Y,3,dotvec);
868: for (j=lx;j<kx;j++) {
869: BVGetColumn(X,j,&z);
870: MatMult(A,z,f);
871: BVRestoreColumn(X,j,&z);
872: ulim = PetscMin(ly+(j-lx)+1,ky);
873: Y->l = 0; Y->k = ulim;
874: (*Y->ops->dotvec)(Y,f,marray+j*ldm);
875: if (symm) {
876: for (i=0;i<j;i++) marray[j+i*ldm] = PetscConj(marray[i+j*ldm]);
877: }
878: }
879: if (!symm) {
880: BVCheckOp(X,1,dotvec);
881: BV_AllocateCoeffs(Y);
882: for (j=ly;j<ky;j++) {
883: BVGetColumn(Y,j,&z);
884: MatMultHermitianTranspose(A,z,f);
885: BVRestoreColumn(Y,j,&z);
886: ulim = PetscMin(lx+(j-ly),kx);
887: X->l = 0; X->k = ulim;
888: (*X->ops->dotvec)(X,f,Y->h);
889: for (i=0;i<ulim;i++) marray[j+i*ldm] = PetscConj(Y->h[i]);
890: }
891: }
892: VecDestroy(&f);
893: X->l = lx; X->k = kx;
894: Y->l = ly; Y->k = ky;
895: PetscFunctionReturn(0);
896: }
898: /*
899: Compute Y^H*A*X= [ -- | Y0'*W1 ]
900: [ Y1'*W0 | Y1'*W1 ]
901: Allocates auxiliary BV to store the result of A*X, then one BVDot902: call for top-right part and another one for bottom part;
903: result placed in marray[*,ldm]
904: */
905: static inline PetscErrorCode BVMatProject_MatMult(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm)906: {
907: PetscInt j,lx,ly,kx,ky;
908: const PetscScalar *harray;
909: Mat H;
910: BV W;
912: lx = X->l; kx = X->k;
913: ly = Y->l; ky = Y->k;
914: BVDuplicate(X,&W);
915: X->l = 0; X->k = kx;
916: W->l = 0; W->k = kx;
917: BVMatMult(X,A,W);
919: /* top-right part, Y0'*AX1 */
920: if (ly>0 && lx<kx) {
921: MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H);
922: W->l = lx; W->k = kx;
923: Y->l = 0; Y->k = ly;
924: BVDot(W,Y,H);
925: MatDenseGetArrayRead(H,&harray);
926: for (j=lx;j<kx;j++) PetscArraycpy(marray+j*ldm,harray+j*ly,ly);
927: MatDenseRestoreArrayRead(H,&harray);
928: MatDestroy(&H);
929: }
931: /* bottom part, Y1'*AX */
932: if (kx>0 && ly<ky) {
933: MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H);
934: W->l = 0; W->k = kx;
935: Y->l = ly; Y->k = ky;
936: BVDot(W,Y,H);
937: MatDenseGetArrayRead(H,&harray);
938: for (j=0;j<kx;j++) PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly);
939: MatDenseRestoreArrayRead(H,&harray);
940: MatDestroy(&H);
941: }
942: BVDestroy(&W);
943: X->l = lx; X->k = kx;
944: Y->l = ly; Y->k = ky;
945: PetscFunctionReturn(0);
946: }
948: /*
949: Compute Y^H*A*X= [ -- | Y0'*W1 ]
950: [ Y1'*W0 | Y1'*W1 ]
951: First stage: allocate auxiliary BV to store A*X1, one BVDot for right part;
952: Second stage: resize BV to accommodate A'*Y1, then call BVDot for transpose of
953: bottom-left part; result placed in marray[*,ldm]
954: */
955: static inline PetscErrorCode BVMatProject_MatMult_2(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)956: {
957: PetscInt i,j,lx,ly,kx,ky;
958: const PetscScalar *harray;
959: Mat H;
960: BV W;
962: lx = X->l; kx = X->k;
963: ly = Y->l; ky = Y->k;
965: /* right part, Y'*AX1 */
966: BVDuplicateResize(X,kx-lx,&W);
967: if (ky>0 && lx<kx) {
968: BVMatMult(X,A,W);
969: MatCreateSeqDense(PETSC_COMM_SELF,ky,kx-lx,NULL,&H);
970: Y->l = 0; Y->k = ky;
971: BVDot(W,Y,H);
972: MatDenseGetArrayRead(H,&harray);
973: for (j=lx;j<kx;j++) PetscArraycpy(marray+j*ldm,harray+(j-lx)*ky,ky);
974: MatDenseRestoreArrayRead(H,&harray);
975: MatDestroy(&H);
976: }
978: /* bottom-left part, Y1'*AX0 */
979: if (lx>0 && ly<ky) {
980: if (symm) {
981: /* do not compute, just copy symmetric elements */
982: for (i=ly;i<ky;i++) {
983: for (j=0;j<lx;j++) marray[i+j*ldm] = PetscConj(marray[j+i*ldm]);
984: }
985: } else {
986: BVResize(W,ky-ly,PETSC_FALSE);
987: Y->l = ly; Y->k = ky;
988: BVMatMultHermitianTranspose(Y,A,W);
989: MatCreateSeqDense(PETSC_COMM_SELF,lx,ky-ly,NULL,&H);
990: X->l = 0; X->k = lx;
991: BVDot(W,X,H);
992: MatDenseGetArrayRead(H,&harray);
993: for (i=0;i<ky-ly;i++) {
994: for (j=0;j<lx;j++) {
995: marray[i+j*ldm+ly] = PetscConj(harray[j+i*(ky-ly)]);
996: }
997: }
998: MatDenseRestoreArrayRead(H,&harray);
999: MatDestroy(&H);
1000: }
1001: }
1002: BVDestroy(&W);
1003: X->l = lx; X->k = kx;
1004: Y->l = ly; Y->k = ky;
1005: PetscFunctionReturn(0);
1006: }
1008: /*
1009: Compute Y^H*X = [ -- | Y0'*X1 ] (X contains A*X):
1010: [ Y1'*X0 | Y1'*X1 ]
1011: one BVDot call for top-right part and another one for bottom part;
1012: result placed in marray[*,ldm]
1013: */
1014: static inline PetscErrorCode BVMatProject_Dot(BV X,BV Y,PetscScalar *marray,PetscInt ldm)1015: {
1016: PetscInt j,lx,ly,kx,ky;
1017: const PetscScalar *harray;
1018: Mat H;
1020: lx = X->l; kx = X->k;
1021: ly = Y->l; ky = Y->k;
1023: /* top-right part, Y0'*X1 */
1024: if (ly>0 && lx<kx) {
1025: MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H);
1026: X->l = lx; X->k = kx;
1027: Y->l = 0; Y->k = ly;
1028: BVDot(X,Y,H);
1029: MatDenseGetArrayRead(H,&harray);
1030: for (j=lx;j<kx;j++) PetscArraycpy(marray+j*ldm,harray+j*ly,ly);
1031: MatDenseRestoreArrayRead(H,&harray);
1032: MatDestroy(&H);
1033: }
1035: /* bottom part, Y1'*X */
1036: if (kx>0 && ly<ky) {
1037: MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H);
1038: X->l = 0; X->k = kx;
1039: Y->l = ly; Y->k = ky;
1040: BVDot(X,Y,H);
1041: MatDenseGetArrayRead(H,&harray);
1042: for (j=0;j<kx;j++) PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly);
1043: MatDenseRestoreArrayRead(H,&harray);
1044: MatDestroy(&H);
1045: }
1046: X->l = lx; X->k = kx;
1047: Y->l = ly; Y->k = ky;
1048: PetscFunctionReturn(0);
1049: }
1051: /*@
1052: BVMatProject - Computes the projection of a matrix onto a subspace.
1054: Collective on X
1056: Input Parameters:
1057: + X - basis vectors
1058: . A - (optional) matrix to be projected
1059: . Y - left basis vectors, can be equal to X
1060: - M - Mat object where the result must be placed
1062: Output Parameter:
1063: . M - the resulting matrix
1065: Notes:
1066: If A=NULL, then it is assumed that X already contains A*X.
1068: This operation is similar to BVDot(), with important differences.
1069: The goal is to compute the matrix resulting from the orthogonal projection
1070: of A onto the subspace spanned by the columns of X, M = X^H*A*X, or the
1071: oblique projection onto X along Y, M = Y^H*A*X.
1073: A difference with respect to BVDot() is that the standard inner product
1074: is always used, regardless of a non-standard inner product being specified
1075: with BVSetMatrix().
1077: On entry, M must be a sequential dense Mat with dimensions ky,kx at least,
1078: where ky (resp. kx) is the number of active columns of Y (resp. X).
1079: Another difference with respect to BVDot() is that all entries of M are
1080: computed except the leading ly,lx part, where ly (resp. lx) is the
1081: number of leading columns of Y (resp. X). Hence, the leading columns of
1082: X and Y participate in the computation, as opposed to BVDot().
1083: The leading part of M is assumed to be already available from previous
1084: computations.
1086: In the orthogonal projection case, Y=X, some computation can be saved if
1087: A is real symmetric (or complex Hermitian). In order to exploit this
1088: property, the symmetry flag of A must be set with MatSetOption().
1090: Level: intermediate
1092: .seealso: BVDot(), BVSetActiveColumns(), BVSetMatrix()
1093: @*/
1094: PetscErrorCode BVMatProject(BV X,Mat A,BV Y,Mat M)1095: {
1096: PetscBool set,flg,symm=PETSC_FALSE;
1097: PetscInt m,n;
1098: PetscScalar *marray;
1099: Mat Xmatrix,Ymatrix;
1100: PetscObjectId idx,idy;
1107: BVCheckSizes(X,1);
1108: if (A) {
1111: }
1113: BVCheckSizes(Y,3);
1118: MatGetSize(M,&m,&n);
1123: PetscLogEventBegin(BV_MatProject,X,A,Y,0);
1124: /* temporarily set standard inner product */
1125: Xmatrix = X->matrix;
1126: Ymatrix = Y->matrix;
1127: X->matrix = Y->matrix = NULL;
1129: PetscObjectGetId((PetscObject)X,&idx);
1130: PetscObjectGetId((PetscObject)Y,&idy);
1131: if (A && idx==idy) { /* check symmetry of M=X'AX */
1132: MatIsHermitianKnown(A,&set,&flg);
1133: symm = set? flg: PETSC_FALSE;
1134: }
1136: MatDenseGetArray(M,&marray);
1138: if (A) {
1139: if (X->vmm==BV_MATMULT_VECS) {
1140: /* perform computation column by column */
1141: BVMatProject_Vec(X,A,Y,marray,m,symm);
1142: } else {
1143: /* use BVMatMult, then BVDot */
1144: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&flg);
1145: if (symm || (flg && X->l>=X->k/2 && Y->l>=Y->k/2)) BVMatProject_MatMult_2(X,A,Y,marray,m,symm);
1146: else BVMatProject_MatMult(X,A,Y,marray,m);
1147: }
1148: } else {
1149: /* use BVDot on subblocks */
1150: BVMatProject_Dot(X,Y,marray,m);
1151: }
1153: MatDenseRestoreArray(M,&marray);
1154: PetscLogEventEnd(BV_MatProject,X,A,Y,0);
1155: /* restore non-standard inner product */
1156: X->matrix = Xmatrix;
1157: Y->matrix = Ymatrix;
1158: PetscFunctionReturn(0);
1159: }