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: BV operations involving global communication
12: */
14: #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
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: PETSC_STATIC_INLINE PetscErrorCode BVDot_Private(BV X,BV Y,Mat M) 21: {
23: PetscObjectId idx,idy;
24: PetscInt i,j,jend,m;
25: PetscScalar *marray;
26: PetscBool symm=PETSC_FALSE;
27: Mat B;
28: Vec z;
31: MatGetSize(M,&m,NULL);
32: MatDenseGetArray(M,&marray);
33: PetscObjectGetId((PetscObject)X,&idx);
34: PetscObjectGetId((PetscObject)Y,&idy);
35: B = Y->matrix;
36: Y->matrix = NULL;
37: if (idx==idy) symm=PETSC_TRUE; /* M=X'BX is symmetric */
38: jend = X->k;
39: for (j=X->l;j<jend;j++) {
40: if (symm) Y->k = j+1;
41: BVGetColumn(X->cached,j,&z);
42: (*Y->ops->dotvec)(Y,z,marray+j*m+Y->l);
43: BVRestoreColumn(X->cached,j,&z);
44: if (symm) {
45: for (i=X->l;i<j;i++)
46: marray[j+i*m] = PetscConj(marray[i+j*m]);
47: }
48: }
49: MatDenseRestoreArray(M,&marray);
50: Y->matrix = B;
51: return(0);
52: }
54: /*@
55: BVDot - Computes the 'block-dot' product of two basis vectors objects.
57: Collective on BV 59: Input Parameters:
60: + X, Y - basis vectors
61: - M - Mat object where the result must be placed
63: Output Parameter:
64: . M - the resulting matrix
66: Notes:
67: This is the generalization of VecDot() for a collection of vectors, M = Y^H*X.
68: The result is a matrix M whose entry m_ij is equal to y_i^H x_j (where y_i^H
69: denotes the conjugate transpose of y_i).
71: If a non-standard inner product has been specified with BVSetMatrix(),
72: then the result is M = Y^H*B*X. In this case, both X and Y must have
73: the same associated matrix.
75: On entry, M must be a sequential dense Mat with dimensions m,n at least, where
76: m is the number of active columns of Y and n is the number of active columns of X.
77: Only rows (resp. columns) of M starting from ly (resp. lx) are computed,
78: where ly (resp. lx) is the number of leading columns of Y (resp. X).
80: X and Y need not be different objects.
82: Level: intermediate
84: .seealso: BVDotVec(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
85: @*/
86: PetscErrorCode BVDot(BV X,BV Y,Mat M) 87: {
89: PetscBool match;
90: PetscInt m,n;
97: BVCheckSizes(X,1);
99: BVCheckSizes(Y,2);
102: PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&match);
103: if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Mat argument must be of type seqdense");
105: MatGetSize(M,&m,&n);
106: if (m<Y->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %D rows, should have at least %D",m,Y->k);
107: if (n<X->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %D columns, should have at least %D",n,X->k);
108: if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);
109: if (X->matrix!=Y->matrix) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"X and Y must have the same inner product matrix");
110: if (X->l==X->k || Y->l==Y->k) return(0);
112: PetscLogEventBegin(BV_Dot,X,Y,0,0);
113: if (X->matrix) { /* non-standard inner product */
114: /* compute BX first */
115: BV_IPMatMultBV(X);
116: if (X->vmm==BV_MATMULT_VECS) {
117: /* perform computation column by column */
118: BVDot_Private(X,Y,M);
119: } else {
120: (*X->ops->dot)(X->cached,Y,M);
121: }
122: } else {
123: (*X->ops->dot)(X,Y,M);
124: }
125: PetscLogEventEnd(BV_Dot,X,Y,0,0);
126: return(0);
127: }
129: /*@
130: BVDotVec - Computes multiple dot products of a vector against all the
131: column vectors of a BV.
133: Collective on BV and Vec
135: Input Parameters:
136: + X - basis vectors
137: - y - a vector
139: Output Parameter:
140: . m - an array where the result must be placed
142: Notes:
143: This is analogue to VecMDot(), but using BV to represent a collection
144: of vectors. The result is m = X^H*y, so m_i is equal to x_j^H y. Note
145: that here X is transposed as opposed to BVDot().
147: If a non-standard inner product has been specified with BVSetMatrix(),
148: then the result is m = X^H*B*y.
150: The length of array m must be equal to the number of active columns of X
151: minus the number of leading columns, i.e. the first entry of m is the
152: product of the first non-leading column with y.
154: Level: intermediate
156: .seealso: BVDot(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
157: @*/
158: PetscErrorCode BVDotVec(BV X,Vec y,PetscScalar m[])159: {
161: PetscInt n;
167: BVCheckSizes(X,1);
171: VecGetLocalSize(y,&n);
172: if (X->n!=n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, y %D",X->n,n);
174: PetscLogEventBegin(BV_DotVec,X,y,0,0);
175: (*X->ops->dotvec)(X,y,m);
176: PetscLogEventEnd(BV_DotVec,X,y,0,0);
177: return(0);
178: }
180: /*@
181: BVDotVecBegin - Starts a split phase dot product computation.
183: Input Parameters:
184: + X - basis vectors
185: . y - a vector
186: - m - an array where the result will go (can be NULL)
188: Note:
189: Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
191: Level: advanced
193: .seealso: BVDotVecEnd(), BVDotVec()
194: @*/
195: PetscErrorCode BVDotVecBegin(BV X,Vec y,PetscScalar *m)196: {
197: PetscErrorCode ierr;
198: PetscInt i,n,nv;
199: PetscSplitReduction *sr;
200: MPI_Comm comm;
206: BVCheckSizes(X,1);
210: VecGetLocalSize(y,&n);
211: if (X->n!=n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, y %D",X->n,n);
213: if (X->ops->dotvec_begin) {
214: (*X->ops->dotvec_begin)(X,y,m);
215: } else {
216: nv = X->k-X->l;
217: PetscObjectGetComm((PetscObject)X,&comm);
218: PetscSplitReductionGet(comm,&sr);
219: if (sr->state != STATE_BEGIN) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
220: for (i=0;i<nv;i++) {
221: if (sr->numopsbegin+i >= sr->maxops) {
222: PetscSplitReductionExtend(sr);
223: }
224: sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
225: sr->invecs[sr->numopsbegin+i] = (void*)X;
226: }
227: PetscLogEventBegin(BV_DotVec,X,y,0,0);
228: (*X->ops->dotvec_local)(X,y,sr->lvalues+sr->numopsbegin);
229: sr->numopsbegin += nv;
230: PetscLogEventEnd(BV_DotVec,X,y,0,0);
231: }
232: return(0);
233: }
235: /*@
236: BVDotVecEnd - Ends a split phase dot product computation.
238: Input Parameters:
239: + X - basis vectors
240: . y - a vector
241: - m - an array where the result will go
243: Note:
244: Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
246: Level: advanced
248: .seealso: BVDotVecBegin(), BVDotVec()
249: @*/
250: PetscErrorCode BVDotVecEnd(BV X,Vec y,PetscScalar *m)251: {
252: PetscErrorCode ierr;
253: PetscInt i,nv;
254: PetscSplitReduction *sr;
255: MPI_Comm comm;
260: BVCheckSizes(X,1);
262: if (X->ops->dotvec_end) {
263: (*X->ops->dotvec_end)(X,y,m);
264: } else {
265: nv = X->k-X->l;
266: PetscObjectGetComm((PetscObject)X,&comm);
267: PetscSplitReductionGet(comm,&sr);
268: PetscSplitReductionEnd(sr);
270: if (sr->numopsend >= sr->numopsbegin) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
271: if ((void*)X != sr->invecs[sr->numopsend]) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
272: if (sr->reducetype[sr->numopsend] != PETSC_SR_REDUCE_SUM) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
273: for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
275: /* Finished getting all the results so reset to no outstanding requests */
276: if (sr->numopsend == sr->numopsbegin) {
277: sr->state = STATE_BEGIN;
278: sr->numopsend = 0;
279: sr->numopsbegin = 0;
280: }
281: }
282: return(0);
283: }
285: /*@
286: BVDotColumn - Computes multiple dot products of a column against all the
287: previous columns of a BV.
289: Collective on BV291: Input Parameters:
292: + X - basis vectors
293: - j - the column index
295: Output Parameter:
296: . q - an array where the result must be placed
298: Notes:
299: This operation is equivalent to BVDotVec() but it uses column j of X
300: rather than taking a Vec as an argument. The number of active columns of
301: X is set to j before the computation, and restored afterwards.
302: If X has leading columns specified, then these columns do not participate
303: in the computation. Therefore, the length of array q must be equal to j
304: minus the number of leading columns.
306: Developer Notes:
307: If q is NULL, then the result is written in position nc+l of the internal
308: buffer vector, see BVGetBufferVec().
310: Level: advanced
312: .seealso: BVDot(), BVDotVec(), BVSetActiveColumns(), BVSetMatrix()
313: @*/
314: PetscErrorCode BVDotColumn(BV X,PetscInt j,PetscScalar *q)315: {
317: PetscInt ksave;
318: Vec y;
324: BVCheckSizes(X,1);
326: if (j<0) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
327: if (j>=X->m) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,X->m);
329: PetscLogEventBegin(BV_DotVec,X,0,0,0);
330: ksave = X->k;
331: X->k = j;
332: if (!q && !X->buffer) { BVGetBufferVec(X,&X->buffer); }
333: BVGetColumn(X,j,&y);
334: (*X->ops->dotvec)(X,y,q);
335: BVRestoreColumn(X,j,&y);
336: X->k = ksave;
337: PetscLogEventEnd(BV_DotVec,X,0,0,0);
338: return(0);
339: }
341: /*@
342: BVDotColumnBegin - Starts a split phase dot product computation.
344: Input Parameters:
345: + X - basis vectors
346: - j - the column index
347: - m - an array where the result will go (can be NULL)
349: Note:
350: Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
352: Level: advanced
354: .seealso: BVDotColumnEnd(), BVDotColumn()
355: @*/
356: PetscErrorCode BVDotColumnBegin(BV X,PetscInt j,PetscScalar *m)357: {
358: PetscErrorCode ierr;
359: PetscInt i,nv,ksave;
360: PetscSplitReduction *sr;
361: MPI_Comm comm;
362: Vec y;
368: BVCheckSizes(X,1);
370: if (j<0) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
371: if (j>=X->m) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,X->m);
372: ksave = X->k;
373: X->k = j;
374: BVGetColumn(X,j,&y);
376: if (X->ops->dotvec_begin) {
377: (*X->ops->dotvec_begin)(X,y,m);
378: } else {
379: nv = X->k-X->l;
380: PetscObjectGetComm((PetscObject)X,&comm);
381: PetscSplitReductionGet(comm,&sr);
382: if (sr->state != STATE_BEGIN) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
383: for (i=0;i<nv;i++) {
384: if (sr->numopsbegin+i >= sr->maxops) {
385: PetscSplitReductionExtend(sr);
386: }
387: sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
388: sr->invecs[sr->numopsbegin+i] = (void*)X;
389: }
390: PetscLogEventBegin(BV_DotVec,X,0,0,0);
391: (*X->ops->dotvec_local)(X,y,sr->lvalues+sr->numopsbegin);
392: sr->numopsbegin += nv;
393: PetscLogEventEnd(BV_DotVec,X,0,0,0);
394: }
395: BVRestoreColumn(X,j,&y);
396: X->k = ksave;
397: return(0);
398: }
400: /*@
401: BVDotColumnEnd - Ends a split phase dot product computation.
403: Input Parameters:
404: + X - basis vectors
405: . j - the column index
406: - m - an array where the result will go
408: Notes:
409: Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
411: Level: advanced
413: .seealso: BVDotColumnBegin(), BVDotColumn()
414: @*/
415: PetscErrorCode BVDotColumnEnd(BV X,PetscInt j,PetscScalar *m)416: {
417: PetscErrorCode ierr;
418: PetscInt i,nv,ksave;
419: PetscSplitReduction *sr;
420: MPI_Comm comm;
421: Vec y;
427: BVCheckSizes(X,1);
429: if (j<0) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
430: if (j>=X->m) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,X->m);
431: ksave = X->k;
432: X->k = j;
434: if (X->ops->dotvec_end) {
435: BVGetColumn(X,j,&y);
436: (*X->ops->dotvec_end)(X,y,m);
437: BVRestoreColumn(X,j,&y);
438: } else {
439: nv = X->k-X->l;
440: PetscObjectGetComm((PetscObject)X,&comm);
441: PetscSplitReductionGet(comm,&sr);
442: PetscSplitReductionEnd(sr);
444: if (sr->numopsend >= sr->numopsbegin) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
445: if ((void*)X != sr->invecs[sr->numopsend]) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
446: if (sr->reducetype[sr->numopsend] != PETSC_SR_REDUCE_SUM) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
447: for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
449: /* Finished getting all the results so reset to no outstanding requests */
450: if (sr->numopsend == sr->numopsbegin) {
451: sr->state = STATE_BEGIN;
452: sr->numopsend = 0;
453: sr->numopsbegin = 0;
454: }
455: }
456: X->k = ksave;
457: return(0);
458: }
460: PETSC_STATIC_INLINE PetscErrorCode BVNorm_Private(BV bv,Vec z,NormType type,PetscReal *val)461: {
463: PetscScalar p;
466: BV_IPMatMult(bv,z);
467: VecDot(bv->Bx,z,&p);
468: BV_SafeSqrt(bv,p,val);
469: return(0);
470: }
472: PETSC_STATIC_INLINE PetscErrorCode BVNorm_Begin_Private(BV bv,Vec z,NormType type,PetscReal *val)473: {
475: PetscScalar p;
478: BV_IPMatMult(bv,z);
479: VecDotBegin(bv->Bx,z,&p);
480: return(0);
481: }
483: PETSC_STATIC_INLINE PetscErrorCode BVNorm_End_Private(BV bv,Vec z,NormType type,PetscReal *val)484: {
486: PetscScalar p;
489: VecDotEnd(bv->Bx,z,&p);
490: BV_SafeSqrt(bv,p,val);
491: return(0);
492: }
494: /*@
495: BVNorm - Computes the matrix norm of the BV.
497: Collective on BV499: Input Parameters:
500: + bv - basis vectors
501: - type - the norm type
503: Output Parameter:
504: . val - the norm
506: Notes:
507: All active columns (except the leading ones) are considered as a matrix.
508: The allowed norms are NORM_1, NORM_FROBENIUS, and NORM_INFINITY.
510: This operation fails if a non-standard inner product has been
511: specified with BVSetMatrix().
513: Level: intermediate
515: .seealso: BVNormVec(), BVNormColumn(), BVSetActiveColumns(), BVSetMatrix()
516: @*/
517: PetscErrorCode BVNorm(BV bv,NormType type,PetscReal *val)518: {
526: BVCheckSizes(bv,1);
528: if (type==NORM_2 || type==NORM_1_AND_2) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
529: if (bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Matrix norm not available for non-standard inner product");
531: PetscLogEventBegin(BV_Norm,bv,0,0,0);
532: (*bv->ops->norm)(bv,-1,type,val);
533: PetscLogEventEnd(BV_Norm,bv,0,0,0);
534: return(0);
535: }
537: /*@
538: BVNormVec - Computes the norm of a given vector.
540: Collective on BV542: Input Parameters:
543: + bv - basis vectors
544: . v - the vector
545: - type - the norm type
547: Output Parameter:
548: . val - the norm
550: Notes:
551: This is the analogue of BVNormColumn() but for a vector that is not in the BV.
552: If a non-standard inner product has been specified with BVSetMatrix(),
553: then the returned value is sqrt(v'*B*v), where B is the inner product
554: matrix (argument 'type' is ignored). Otherwise, VecNorm() is called.
556: Level: developer
558: .seealso: BVNorm(), BVNormColumn(), BVSetMatrix()
559: @*/
560: PetscErrorCode BVNormVec(BV bv,Vec v,NormType type,PetscReal *val)561: {
563: PetscInt n;
571: BVCheckSizes(bv,1);
575: if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
577: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
578: if (bv->matrix) { /* non-standard inner product */
579: VecGetLocalSize(v,&n);
580: if (bv->n!=n) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %D, v %D",bv->n,n);
581: BVNorm_Private(bv,v,type,val);
582: } else {
583: VecNorm(v,type,val);
584: }
585: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
586: return(0);
587: }
589: /*@
590: BVNormVecBegin - Starts a split phase norm computation.
592: Input Parameters:
593: + bv - basis vectors
594: . v - the vector
595: . type - the norm type
596: - val - the norm
598: Note:
599: Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
601: Level: advanced
603: .seealso: BVNormVecEnd(), BVNormVec()
604: @*/
605: PetscErrorCode BVNormVecBegin(BV bv,Vec v,NormType type,PetscReal *val)606: {
608: PetscInt n;
616: BVCheckSizes(bv,1);
620: if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
622: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
623: if (bv->matrix) { /* non-standard inner product */
624: VecGetLocalSize(v,&n);
625: if (bv->n!=n) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %D, v %D",bv->n,n);
626: BVNorm_Begin_Private(bv,v,type,val);
627: } else {
628: VecNormBegin(v,type,val);
629: }
630: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
631: return(0);
632: }
634: /*@
635: BVNormVecEnd - Ends a split phase norm computation.
637: Input Parameters:
638: + bv - basis vectors
639: . v - the vector
640: . type - the norm type
641: - val - the norm
643: Note:
644: Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
646: Level: advanced
648: .seealso: BVNormVecBegin(), BVNormVec()
649: @*/
650: PetscErrorCode BVNormVecEnd(BV bv,Vec v,NormType type,PetscReal *val)651: {
659: BVCheckSizes(bv,1);
661: if (type==NORM_1_AND_2) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
663: if (bv->matrix) { /* non-standard inner product */
664: BVNorm_End_Private(bv,v,type,val);
665: } else {
666: VecNormEnd(v,type,val);
667: }
668: return(0);
669: }
671: /*@
672: BVNormColumn - Computes the vector norm of a selected column.
674: Collective on BV676: Input Parameters:
677: + bv - basis vectors
678: . j - column number to be used
679: - type - the norm type
681: Output Parameter:
682: . val - the norm
684: Notes:
685: The norm of V[j] is computed (NORM_1, NORM_2, or NORM_INFINITY).
686: If a non-standard inner product has been specified with BVSetMatrix(),
687: then the returned value is sqrt(V[j]'*B*V[j]),
688: where B is the inner product matrix (argument 'type' is ignored).
690: Level: intermediate
692: .seealso: BVNorm(), BVNormVec(), BVSetActiveColumns(), BVSetMatrix()
693: @*/
694: PetscErrorCode BVNormColumn(BV bv,PetscInt j,NormType type,PetscReal *val)695: {
697: Vec z;
705: BVCheckSizes(bv,1);
707: if (j<0 || j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %D, the number of columns is %D",j,bv->m);
708: if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
710: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
711: if (bv->matrix) { /* non-standard inner product */
712: BVGetColumn(bv,j,&z);
713: BVNorm_Private(bv,z,type,val);
714: BVRestoreColumn(bv,j,&z);
715: } else {
716: (*bv->ops->norm)(bv,j,type,val);
717: }
718: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
719: return(0);
720: }
722: /*@
723: BVNormColumnBegin - Starts a split phase norm computation.
725: Input Parameters:
726: + bv - basis vectors
727: . j - column number to be used
728: . type - the norm type
729: - val - the norm
731: Note:
732: Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
734: Level: advanced
736: .seealso: BVNormColumnEnd(), BVNormColumn()
737: @*/
738: PetscErrorCode BVNormColumnBegin(BV bv,PetscInt j,NormType type,PetscReal *val)739: {
740: PetscErrorCode ierr;
741: PetscSplitReduction *sr;
742: PetscReal lresult;
743: MPI_Comm comm;
744: Vec z;
752: BVCheckSizes(bv,1);
754: if (j<0 || j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %D, the number of columns is %D",j,bv->m);
755: if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
757: PetscLogEventBegin(BV_NormVec,bv,0,0,0);
758: BVGetColumn(bv,j,&z);
759: if (bv->matrix) { /* non-standard inner product */
760: BVNorm_Begin_Private(bv,z,type,val);
761: } else if (bv->ops->norm_begin) {
762: (*bv->ops->norm_begin)(bv,j,type,val);
763: } else {
764: PetscObjectGetComm((PetscObject)z,&comm);
765: PetscSplitReductionGet(comm,&sr);
766: if (sr->state != STATE_BEGIN) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
767: if (sr->numopsbegin >= sr->maxops) {
768: PetscSplitReductionExtend(sr);
769: }
770: sr->invecs[sr->numopsbegin] = (void*)bv;
771: (*bv->ops->norm_local)(bv,j,type,&lresult);
772: if (type == NORM_2) lresult = lresult*lresult;
773: if (type == NORM_MAX) sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_MAX;
774: else sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_SUM;
775: sr->lvalues[sr->numopsbegin++] = lresult;
776: }
777: BVRestoreColumn(bv,j,&z);
778: PetscLogEventEnd(BV_NormVec,bv,0,0,0);
779: return(0);
780: }
782: /*@
783: BVNormColumnEnd - Ends a split phase norm computation.
785: Input Parameters:
786: + bv - basis vectors
787: . j - column number to be used
788: . type - the norm type
789: - val - the norm
791: Note:
792: Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
794: Level: advanced
796: .seealso: BVNormColumnBegin(), BVNormColumn()
797: @*/
798: PetscErrorCode BVNormColumnEnd(BV bv,PetscInt j,NormType type,PetscReal *val)799: {
800: PetscErrorCode ierr;
801: PetscSplitReduction *sr;
802: MPI_Comm comm;
803: Vec z;
811: BVCheckSizes(bv,1);
813: if (type==NORM_1_AND_2) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
815: BVGetColumn(bv,j,&z);
816: if (bv->matrix) { /* non-standard inner product */
817: BVNorm_End_Private(bv,z,type,val);
818: } else if (bv->ops->norm_end) {
819: (*bv->ops->norm_end)(bv,j,type,val);
820: } else {
821: PetscObjectGetComm((PetscObject)z,&comm);
822: PetscSplitReductionGet(comm,&sr);
823: PetscSplitReductionEnd(sr);
825: if (sr->numopsend >= sr->numopsbegin) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times then VecxxxBegin()");
826: if ((void*)bv != sr->invecs[sr->numopsend]) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() in a different order or with a different vector than VecxxxBegin()");
827: if (sr->reducetype[sr->numopsend] != PETSC_SR_REDUCE_MAX && type == NORM_MAX) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called BVNormEnd(,NORM_MAX,) on a reduction started with VecDotBegin() or NORM_1 or NORM_2");
828: *val = PetscRealPart(sr->gvalues[sr->numopsend++]);
829: if (type == NORM_2) *val = PetscSqrtReal(*val);
830: if (sr->numopsend == sr->numopsbegin) {
831: sr->state = STATE_BEGIN;
832: sr->numopsend = 0;
833: sr->numopsbegin = 0;
834: }
835: }
836: BVRestoreColumn(bv,j,&z);
837: return(0);
838: }
840: /*
841: Compute Y^H*A*X: right part column by column (with MatMult) and bottom
842: part row by row (with MatMultHermitianTranspose); result placed in marray[*,ldm]
843: */
844: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_Vec(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)845: {
847: PetscInt i,j,lx,ly,kx,ky,ulim;
848: Vec z,f;
851: lx = X->l; kx = X->k;
852: ly = Y->l; ky = Y->k;
853: BVCreateVec(X,&f);
854: for (j=lx;j<kx;j++) {
855: BVGetColumn(X,j,&z);
856: MatMult(A,z,f);
857: BVRestoreColumn(X,j,&z);
858: ulim = PetscMin(ly+(j-lx)+1,ky);
859: Y->l = 0; Y->k = ulim;
860: (*Y->ops->dotvec)(Y,f,marray+j*ldm);
861: if (symm) {
862: for (i=0;i<j;i++) marray[j+i*ldm] = PetscConj(marray[i+j*ldm]);
863: }
864: }
865: if (!symm) {
866: BV_AllocateCoeffs(Y);
867: for (j=ly;j<ky;j++) {
868: BVGetColumn(Y,j,&z);
869: MatMultHermitianTranspose(A,z,f);
870: BVRestoreColumn(Y,j,&z);
871: ulim = PetscMin(lx+(j-ly),kx);
872: X->l = 0; X->k = ulim;
873: (*X->ops->dotvec)(X,f,Y->h);
874: for (i=0;i<ulim;i++) marray[j+i*ldm] = PetscConj(Y->h[i]);
875: }
876: }
877: VecDestroy(&f);
878: X->l = lx; X->k = kx;
879: Y->l = ly; Y->k = ky;
880: return(0);
881: }
883: /*
884: Compute Y^H*A*X= [ -- | Y0'*W1 ]
885: [ Y1'*W0 | Y1'*W1 ]
886: Allocates auxiliary BV to store the result of A*X, then one BVDot887: call for top-right part and another one for bottom part;
888: result placed in marray[*,ldm]
889: */
890: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_MatMult(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm)891: {
893: PetscInt j,lx,ly,kx,ky;
894: PetscScalar *harray;
895: Mat H;
896: BV W;
899: lx = X->l; kx = X->k;
900: ly = Y->l; ky = Y->k;
901: BVDuplicate(X,&W);
902: X->l = 0; X->k = kx;
903: W->l = 0; W->k = kx;
904: BVMatMult(X,A,W);
906: /* top-right part, Y0'*AX1 */
907: if (ly>0 && lx<kx) {
908: MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H);
909: W->l = lx; W->k = kx;
910: Y->l = 0; Y->k = ly;
911: BVDot(W,Y,H);
912: MatDenseGetArray(H,&harray);
913: for (j=lx;j<kx;j++) {
914: PetscMemcpy(marray+j*ldm,harray+j*ly,ly*sizeof(PetscScalar));
915: }
916: MatDenseRestoreArray(H,&harray);
917: MatDestroy(&H);
918: }
920: /* bottom part, Y1'*AX */
921: if (kx>0 && ly<ky) {
922: MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H);
923: W->l = 0; W->k = kx;
924: Y->l = ly; Y->k = ky;
925: BVDot(W,Y,H);
926: MatDenseGetArray(H,&harray);
927: for (j=0;j<kx;j++) {
928: PetscMemcpy(marray+j*ldm+ly,harray+j*ky+ly,(ky-ly)*sizeof(PetscScalar));
929: }
930: MatDenseRestoreArray(H,&harray);
931: MatDestroy(&H);
932: }
933: BVDestroy(&W);
934: X->l = lx; X->k = kx;
935: Y->l = ly; Y->k = ky;
936: return(0);
937: }
939: /*
940: Compute Y^H*A*X= [ -- | Y0'*W1 ]
941: [ Y1'*W0 | Y1'*W1 ]
942: First stage: allocate auxiliary BV to store A*X1, one BVDot for right part;
943: Second stage: resize BV to accomodate A'*Y1, then call BVDot for transpose of
944: bottom-left part; result placed in marray[*,ldm]
945: */
946: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_MatMult_2(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)947: {
949: PetscInt i,j,lx,ly,kx,ky;
950: PetscScalar *harray;
951: Mat H;
952: BV W;
955: lx = X->l; kx = X->k;
956: ly = Y->l; ky = Y->k;
958: /* right part, Y'*AX1 */
959: BVDuplicateResize(X,kx-lx,&W);
960: if (ky>0 && lx<kx) {
961: BVMatMult(X,A,W);
962: MatCreateSeqDense(PETSC_COMM_SELF,ky,kx-lx,NULL,&H);
963: Y->l = 0; Y->k = ky;
964: BVDot(W,Y,H);
965: MatDenseGetArray(H,&harray);
966: for (j=lx;j<kx;j++) {
967: PetscMemcpy(marray+j*ldm,harray+(j-lx)*ky,ky*sizeof(PetscScalar));
968: }
969: MatDenseRestoreArray(H,&harray);
970: MatDestroy(&H);
971: }
973: /* bottom-left part, Y1'*AX0 */
974: if (lx>0 && ly<ky) {
975: if (symm) {
976: /* do not compute, just copy symmetric elements */
977: for (i=ly;i<ky;i++) {
978: for (j=0;j<lx;j++) marray[i+j*ldm] = PetscConj(marray[j+i*ldm]);
979: }
980: } else {
981: BVResize(W,ky-ly,PETSC_FALSE);
982: Y->l = ly; Y->k = ky;
983: BVMatMultHermitianTranspose(Y,A,W);
984: MatCreateSeqDense(PETSC_COMM_SELF,lx,ky-ly,NULL,&H);
985: X->l = 0; X->k = lx;
986: BVDot(W,X,H);
987: MatDenseGetArray(H,&harray);
988: for (i=0;i<ky-ly;i++) {
989: for (j=0;j<lx;j++) {
990: marray[i+j*ldm+ly] = PetscConj(harray[j+i*(ky-ly)]);
991: }
992: }
993: MatDenseRestoreArray(H,&harray);
994: MatDestroy(&H);
995: }
996: }
997: BVDestroy(&W);
998: X->l = lx; X->k = kx;
999: Y->l = ly; Y->k = ky;
1000: return(0);
1001: }
1003: /*
1004: Compute Y^H*X = [ -- | Y0'*X1 ] (X contains A*X):
1005: [ Y1'*X0 | Y1'*X1 ]
1006: one BVDot call for top-right part and another one for bottom part;
1007: result placed in marray[*,ldm]
1008: */
1009: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_Dot(BV X,BV Y,PetscScalar *marray,PetscInt ldm)1010: {
1012: PetscInt j,lx,ly,kx,ky;
1013: PetscScalar *harray;
1014: Mat H;
1017: lx = X->l; kx = X->k;
1018: ly = Y->l; ky = Y->k;
1020: /* top-right part, Y0'*X1 */
1021: if (ly>0 && lx<kx) {
1022: MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H);
1023: X->l = lx; X->k = kx;
1024: Y->l = 0; Y->k = ly;
1025: BVDot(X,Y,H);
1026: MatDenseGetArray(H,&harray);
1027: for (j=lx;j<kx;j++) {
1028: PetscMemcpy(marray+j*ldm,harray+j*ly,ly*sizeof(PetscScalar));
1029: }
1030: MatDenseRestoreArray(H,&harray);
1031: MatDestroy(&H);
1032: }
1034: /* bottom part, Y1'*X */
1035: if (kx>0 && ly<ky) {
1036: MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H);
1037: X->l = 0; X->k = kx;
1038: Y->l = ly; Y->k = ky;
1039: BVDot(X,Y,H);
1040: MatDenseGetArray(H,&harray);
1041: for (j=0;j<kx;j++) {
1042: PetscMemcpy(marray+j*ldm+ly,harray+j*ky+ly,(ky-ly)*sizeof(PetscScalar));
1043: }
1044: MatDenseRestoreArray(H,&harray);
1045: MatDestroy(&H);
1046: }
1047: X->l = lx; X->k = kx;
1048: Y->l = ly; Y->k = ky;
1049: return(0);
1050: }
1052: /*@
1053: BVMatProject - Computes the projection of a matrix onto a subspace.
1055: Collective on BV1057: Input Parameters:
1058: + X - basis vectors
1059: . A - (optional) matrix to be projected
1060: . Y - left basis vectors, can be equal to X
1061: - M - Mat object where the result must be placed
1063: Output Parameter:
1064: . M - the resulting matrix
1066: Notes:
1067: If A=NULL, then it is assumed that X already contains A*X.
1069: This operation is similar to BVDot(), with important differences.
1070: The goal is to compute the matrix resulting from the orthogonal projection
1071: of A onto the subspace spanned by the columns of X, M = X^H*A*X, or the
1072: oblique projection onto X along Y, M = Y^H*A*X.
1074: A difference with respect to BVDot() is that the standard inner product
1075: is always used, regardless of a non-standard inner product being specified
1076: with BVSetMatrix().
1078: On entry, M must be a sequential dense Mat with dimensions ky,kx at least,
1079: where ky (resp. kx) is the number of active columns of Y (resp. X).
1080: Another difference with respect to BVDot() is that all entries of M are
1081: computed except the leading ly,lx part, where ly (resp. lx) is the
1082: number of leading columns of Y (resp. X). Hence, the leading columns of
1083: X and Y participate in the computation, as opposed to BVDot().
1084: The leading part of M is assumed to be already available from previous
1085: computations.
1087: In the orthogonal projection case, Y=X, some computation can be saved if
1088: A is real symmetric (or complex Hermitian). In order to exploit this
1089: property, the symmetry flag of A must be set with MatSetOption().
1091: Level: intermediate
1093: .seealso: BVDot(), BVSetActiveColumns(), BVSetMatrix()
1094: @*/
1095: PetscErrorCode BVMatProject(BV X,Mat A,BV Y,Mat M)1096: {
1098: PetscBool match,set,flg,symm=PETSC_FALSE;
1099: PetscInt m,n;
1100: PetscScalar *marray;
1101: Mat Xmatrix,Ymatrix;
1102: PetscObjectId idx,idy;
1110: BVCheckSizes(X,1);
1111: if (A) {
1114: }
1116: BVCheckSizes(Y,3);
1119: PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&match);
1120: if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Matrix M must be of type seqdense");
1122: MatGetSize(M,&m,&n);
1123: if (m<Y->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %D rows, should have at least %D",m,Y->k);
1124: if (n<X->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %D columns, should have at least %D",n,X->k);
1125: if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);
1127: PetscLogEventBegin(BV_MatProject,X,A,Y,0);
1128: /* temporarily set standard inner product */
1129: Xmatrix = X->matrix;
1130: Ymatrix = Y->matrix;
1131: X->matrix = Y->matrix = NULL;
1133: PetscObjectGetId((PetscObject)X,&idx);
1134: PetscObjectGetId((PetscObject)Y,&idy);
1135: if (A && idx==idy) { /* check symmetry of M=X'AX */
1136: MatIsHermitianKnown(A,&set,&flg);
1137: symm = set? flg: PETSC_FALSE;
1138: }
1140: MatDenseGetArray(M,&marray);
1142: if (A) {
1143: if (X->vmm==BV_MATMULT_VECS) {
1144: /* perform computation column by column */
1145: BVMatProject_Vec(X,A,Y,marray,m,symm);
1146: } else {
1147: /* use BVMatMult, then BVDot */
1148: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&flg);
1149: if (symm || (flg && X->l>=X->k/2 && Y->l>=Y->k/2)) {
1150: BVMatProject_MatMult_2(X,A,Y,marray,m,symm);
1151: } else {
1152: BVMatProject_MatMult(X,A,Y,marray,m);
1153: }
1154: }
1155: } else {
1156: /* use BVDot on subblocks */
1157: BVMatProject_Dot(X,Y,marray,m);
1158: }
1160: MatDenseRestoreArray(M,&marray);
1161: PetscLogEventEnd(BV_MatProject,X,A,Y,0);
1162: /* restore non-standard inner product */
1163: X->matrix = Xmatrix;
1164: Y->matrix = Ymatrix;
1165: return(0);
1166: }