Actual source code: bvglobal.c
slepc-main 2024-11-15
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: PetscFunctionBegin;
30: BVCheckOp(Y,1,dotvec);
31: PetscCall(MatGetSize(M,&m,NULL));
32: PetscCall(MatDenseGetArray(M,&marray));
33: PetscCall(PetscObjectGetId((PetscObject)X,&idx));
34: PetscCall(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: PetscCall(BVGetColumn(X->cached,j,&z));
42: PetscUseTypeMethod(Y,dotvec,z,marray+j*m+Y->l);
43: PetscCall(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: PetscCall(MatDenseRestoreArray(M,&marray));
50: Y->matrix = B;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@
55: BVDot - Computes the 'block-dot' product of two basis vectors objects.
57: Collective
59: Input Parameters:
60: + X - first basis vectors
61: - Y - second basis vectors
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: {
88: PetscInt m,n;
90: PetscFunctionBegin;
95: BVCheckSizes(X,1);
97: BVCheckSizes(Y,2);
99: PetscCheckSameTypeAndComm(X,1,Y,2);
100: SlepcMatCheckSeq(M);
102: PetscCall(MatGetSize(M,&m,&n));
103: PetscCheck(m>=Y->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,Y->k);
104: PetscCheck(n>=X->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,X->k);
105: PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
106: PetscCheck(X->matrix==Y->matrix,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"X and Y must have the same inner product matrix");
107: if (X->l==X->k || Y->l==Y->k) PetscFunctionReturn(PETSC_SUCCESS);
109: PetscCall(PetscLogEventBegin(BV_Dot,X,Y,0,0));
110: if (X->matrix) { /* non-standard inner product */
111: /* compute BX first */
112: PetscCall(BV_IPMatMultBV(X));
113: if (X->vmm==BV_MATMULT_VECS) {
114: /* perform computation column by column */
115: PetscCall(BVDot_Private(X,Y,M));
116: } else PetscUseTypeMethod(X->cached,dot,Y,M);
117: } else PetscUseTypeMethod(X,dot,Y,M);
118: PetscCall(PetscLogEventEnd(BV_Dot,X,Y,0,0));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /*@
123: BVDotVec - Computes multiple dot products of a vector against all the
124: column vectors of a BV.
126: Collective
128: Input Parameters:
129: + X - basis vectors
130: - y - a vector
132: Output Parameter:
133: . m - an array where the result must be placed
135: Notes:
136: This is analogue to VecMDot(), but using BV to represent a collection
137: of vectors. The result is m = X^H*y, so m_i is equal to x_i^H y. Note
138: that here X is transposed as opposed to BVDot().
140: If a non-standard inner product has been specified with BVSetMatrix(),
141: then the result is m = X^H*B*y.
143: The length of array m must be equal to the number of active columns of X
144: minus the number of leading columns, i.e. the first entry of m is the
145: product of the first non-leading column with y.
147: Level: intermediate
149: .seealso: BVDot(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
150: @*/
151: PetscErrorCode BVDotVec(BV X,Vec y,PetscScalar m[])
152: {
153: PetscInt n;
155: PetscFunctionBegin;
159: BVCheckSizes(X,1);
160: BVCheckOp(X,1,dotvec);
162: PetscCheckSameTypeAndComm(X,1,y,2);
164: PetscCall(VecGetLocalSize(y,&n));
165: PetscCheck(X->n==n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", y %" PetscInt_FMT,X->n,n);
167: PetscCall(PetscLogEventBegin(BV_DotVec,X,y,0,0));
168: PetscUseTypeMethod(X,dotvec,y,m);
169: PetscCall(PetscLogEventEnd(BV_DotVec,X,y,0,0));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@
174: BVDotVecBegin - Starts a split phase dot product computation.
176: Input Parameters:
177: + X - basis vectors
178: . y - a vector
179: - m - an array where the result will go (can be NULL)
181: Note:
182: Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
184: Level: advanced
186: .seealso: BVDotVecEnd(), BVDotVec()
187: @*/
188: PetscErrorCode BVDotVecBegin(BV X,Vec y,PetscScalar *m)
189: {
190: PetscInt i,n,nv;
191: PetscSplitReduction *sr;
192: MPI_Comm comm;
194: PetscFunctionBegin;
198: BVCheckSizes(X,1);
200: PetscCheckSameTypeAndComm(X,1,y,2);
202: PetscCall(VecGetLocalSize(y,&n));
203: PetscCheck(X->n==n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", y %" PetscInt_FMT,X->n,n);
205: if (X->ops->dotvec_begin) PetscUseTypeMethod(X,dotvec_begin,y,m);
206: else {
207: BVCheckOp(X,1,dotvec_local);
208: nv = X->k-X->l;
209: PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
210: PetscCall(PetscSplitReductionGet(comm,&sr));
211: PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
212: for (i=0;i<nv;i++) {
213: if (sr->numopsbegin+i >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
214: sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
215: sr->invecs[sr->numopsbegin+i] = (void*)X;
216: }
217: PetscCall(PetscLogEventBegin(BV_DotVec,X,y,0,0));
218: PetscUseTypeMethod(X,dotvec_local,y,sr->lvalues+sr->numopsbegin);
219: sr->numopsbegin += nv;
220: PetscCall(PetscLogEventEnd(BV_DotVec,X,y,0,0));
221: }
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /*@
226: BVDotVecEnd - Ends a split phase dot product computation.
228: Input Parameters:
229: + X - basis vectors
230: . y - a vector
231: - m - an array where the result will go
233: Note:
234: Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
236: Level: advanced
238: .seealso: BVDotVecBegin(), BVDotVec()
239: @*/
240: PetscErrorCode BVDotVecEnd(BV X,Vec y,PetscScalar *m)
241: {
242: PetscInt i,nv;
243: PetscSplitReduction *sr;
244: MPI_Comm comm;
246: PetscFunctionBegin;
249: BVCheckSizes(X,1);
251: if (X->ops->dotvec_end) PetscUseTypeMethod(X,dotvec_end,y,m);
252: else {
253: nv = X->k-X->l;
254: PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
255: PetscCall(PetscSplitReductionGet(comm,&sr));
256: PetscCall(PetscSplitReductionEnd(sr));
258: PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
259: PetscCheck((void*)X==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
260: PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_SUM,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
261: for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
263: /* Finished getting all the results so reset to no outstanding requests */
264: if (sr->numopsend == sr->numopsbegin) {
265: sr->state = STATE_BEGIN;
266: sr->numopsend = 0;
267: sr->numopsbegin = 0;
268: }
269: }
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: BVDotColumn - Computes multiple dot products of a column against all the
275: previous columns of a BV.
277: Collective
279: Input Parameters:
280: + X - basis vectors
281: - j - the column index
283: Output Parameter:
284: . q - an array where the result must be placed
286: Notes:
287: This operation is equivalent to BVDotVec() but it uses column j of X
288: rather than taking a Vec as an argument. The number of active columns of
289: X is set to j before the computation, and restored afterwards.
290: If X has leading columns specified, then these columns do not participate
291: in the computation. Therefore, the length of array q must be equal to j
292: minus the number of leading columns.
294: Developer Notes:
295: If q is NULL, then the result is written in position nc+l of the internal
296: buffer vector, see BVGetBufferVec().
298: Level: advanced
300: .seealso: BVDot(), BVDotVec(), BVSetActiveColumns(), BVSetMatrix()
301: @*/
302: PetscErrorCode BVDotColumn(BV X,PetscInt j,PetscScalar *q)
303: {
304: PetscInt ksave;
305: Vec y;
307: PetscFunctionBegin;
311: BVCheckSizes(X,1);
312: BVCheckOp(X,1,dotvec);
314: PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
315: PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
317: PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
318: ksave = X->k;
319: X->k = j;
320: if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
321: PetscCall(BVGetColumn(X,j,&y));
322: PetscUseTypeMethod(X,dotvec,y,q);
323: PetscCall(BVRestoreColumn(X,j,&y));
324: X->k = ksave;
325: PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: /*@
330: BVDotColumnBegin - Starts a split phase dot product computation.
332: Input Parameters:
333: + X - basis vectors
334: . j - the column index
335: - m - an array where the result will go (can be NULL)
337: Note:
338: Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
340: Level: advanced
342: .seealso: BVDotColumnEnd(), BVDotColumn()
343: @*/
344: PetscErrorCode BVDotColumnBegin(BV X,PetscInt j,PetscScalar *m)
345: {
346: PetscInt i,nv,ksave;
347: PetscSplitReduction *sr;
348: MPI_Comm comm;
349: Vec y;
351: PetscFunctionBegin;
355: BVCheckSizes(X,1);
357: PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
358: PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
359: ksave = X->k;
360: X->k = j;
361: PetscCall(BVGetColumn(X,j,&y));
363: if (X->ops->dotvec_begin) PetscUseTypeMethod(X,dotvec_begin,y,m);
364: else {
365: BVCheckOp(X,1,dotvec_local);
366: nv = X->k-X->l;
367: PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
368: PetscCall(PetscSplitReductionGet(comm,&sr));
369: PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
370: for (i=0;i<nv;i++) {
371: if (sr->numopsbegin+i >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
372: sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
373: sr->invecs[sr->numopsbegin+i] = (void*)X;
374: }
375: PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
376: PetscUseTypeMethod(X,dotvec_local,y,sr->lvalues+sr->numopsbegin);
377: sr->numopsbegin += nv;
378: PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
379: }
380: PetscCall(BVRestoreColumn(X,j,&y));
381: X->k = ksave;
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: BVDotColumnEnd - Ends a split phase dot product computation.
388: Input Parameters:
389: + X - basis vectors
390: . j - the column index
391: - m - an array where the result will go
393: Notes:
394: Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
396: Level: advanced
398: .seealso: BVDotColumnBegin(), BVDotColumn()
399: @*/
400: PetscErrorCode BVDotColumnEnd(BV X,PetscInt j,PetscScalar *m)
401: {
402: PetscInt i,nv,ksave;
403: PetscSplitReduction *sr;
404: MPI_Comm comm;
405: Vec y;
407: PetscFunctionBegin;
411: BVCheckSizes(X,1);
413: PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
414: PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
415: ksave = X->k;
416: X->k = j;
418: if (X->ops->dotvec_end) {
419: PetscCall(BVGetColumn(X,j,&y));
420: PetscUseTypeMethod(X,dotvec_end,y,m);
421: PetscCall(BVRestoreColumn(X,j,&y));
422: } else {
423: nv = X->k-X->l;
424: PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
425: PetscCall(PetscSplitReductionGet(comm,&sr));
426: PetscCall(PetscSplitReductionEnd(sr));
428: PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
429: PetscCheck((void*)X==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
430: PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_SUM,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
431: for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
433: /* Finished getting all the results so reset to no outstanding requests */
434: if (sr->numopsend == sr->numopsbegin) {
435: sr->state = STATE_BEGIN;
436: sr->numopsend = 0;
437: sr->numopsbegin = 0;
438: }
439: }
440: X->k = ksave;
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: static inline PetscErrorCode BVNorm_Private(BV bv,Vec z,NormType type,PetscReal *val)
445: {
446: PetscScalar p;
448: PetscFunctionBegin;
449: PetscCall(BV_IPMatMult(bv,z));
450: PetscCall(VecDot(bv->Bx,z,&p));
451: PetscCall(BV_SafeSqrt(bv,p,val));
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: static inline PetscErrorCode BVNorm_Begin_Private(BV bv,Vec z,NormType type,PetscReal *val)
456: {
457: PetscScalar p;
459: PetscFunctionBegin;
460: PetscCall(BV_IPMatMult(bv,z));
461: PetscCall(VecDotBegin(bv->Bx,z,&p));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: static inline PetscErrorCode BVNorm_End_Private(BV bv,Vec z,NormType type,PetscReal *val)
466: {
467: PetscScalar p;
469: PetscFunctionBegin;
470: PetscCall(VecDotEnd(bv->Bx,z,&p));
471: PetscCall(BV_SafeSqrt(bv,p,val));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: /*@
476: BVNorm - Computes the matrix norm of the BV.
478: Collective
480: Input Parameters:
481: + bv - basis vectors
482: - type - the norm type
484: Output Parameter:
485: . val - the norm
487: Notes:
488: All active columns (except the leading ones) are considered as a matrix.
489: The allowed norms are NORM_1, NORM_FROBENIUS, and NORM_INFINITY.
491: This operation fails if a non-standard inner product has been
492: specified with BVSetMatrix().
494: Level: intermediate
496: .seealso: BVNormVec(), BVNormColumn(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
497: @*/
498: PetscErrorCode BVNorm(BV bv,NormType type,PetscReal *val)
499: {
500: PetscFunctionBegin;
503: PetscAssertPointer(val,3);
505: BVCheckSizes(bv,1);
507: PetscCheck(type!=NORM_2 && type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
508: PetscCheck(!bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Matrix norm not available for non-standard inner product");
510: PetscCall(PetscLogEventBegin(BV_Norm,bv,0,0,0));
511: PetscUseTypeMethod(bv,norm,-1,type,val);
512: PetscCall(PetscLogEventEnd(BV_Norm,bv,0,0,0));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: /*@
517: BVNormVec - Computes the norm of a given vector.
519: Collective
521: Input Parameters:
522: + bv - basis vectors
523: . v - the vector
524: - type - the norm type
526: Output Parameter:
527: . val - the norm
529: Notes:
530: This is the analogue of BVNormColumn() but for a vector that is not in the BV.
531: If a non-standard inner product has been specified with BVSetMatrix(),
532: then the returned value is sqrt(v'*B*v), where B is the inner product
533: matrix (argument 'type' is ignored). Otherwise, VecNorm() is called.
535: Level: developer
537: .seealso: BVNorm(), BVNormColumn(), BVSetMatrix()
538: @*/
539: PetscErrorCode BVNormVec(BV bv,Vec v,NormType type,PetscReal *val)
540: {
541: PetscInt n;
543: PetscFunctionBegin;
547: PetscAssertPointer(val,4);
549: BVCheckSizes(bv,1);
551: PetscCheckSameComm(bv,1,v,2);
553: PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
555: PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
556: if (bv->matrix) { /* non-standard inner product */
557: PetscCall(VecGetLocalSize(v,&n));
558: PetscCheck(bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %" PetscInt_FMT ", v %" PetscInt_FMT,bv->n,n);
559: PetscCall(BVNorm_Private(bv,v,type,val));
560: } else PetscCall(VecNorm(v,type,val));
561: PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: /*@
566: BVNormVecBegin - Starts a split phase norm computation.
568: Input Parameters:
569: + bv - basis vectors
570: . v - the vector
571: . type - the norm type
572: - val - the norm
574: Note:
575: Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
577: Level: advanced
579: .seealso: BVNormVecEnd(), BVNormVec()
580: @*/
581: PetscErrorCode BVNormVecBegin(BV bv,Vec v,NormType type,PetscReal *val)
582: {
583: PetscInt n;
585: PetscFunctionBegin;
589: PetscAssertPointer(val,4);
591: BVCheckSizes(bv,1);
593: PetscCheckSameTypeAndComm(bv,1,v,2);
595: PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
597: PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
598: if (bv->matrix) { /* non-standard inner product */
599: PetscCall(VecGetLocalSize(v,&n));
600: PetscCheck(bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %" PetscInt_FMT ", v %" PetscInt_FMT,bv->n,n);
601: PetscCall(BVNorm_Begin_Private(bv,v,type,val));
602: } else PetscCall(VecNormBegin(v,type,val));
603: PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: /*@
608: BVNormVecEnd - Ends a split phase norm computation.
610: Input Parameters:
611: + bv - basis vectors
612: . v - the vector
613: . type - the norm type
614: - val - the norm
616: Note:
617: Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
619: Level: advanced
621: .seealso: BVNormVecBegin(), BVNormVec()
622: @*/
623: PetscErrorCode BVNormVecEnd(BV bv,Vec v,NormType type,PetscReal *val)
624: {
625: PetscFunctionBegin;
628: PetscAssertPointer(val,4);
630: BVCheckSizes(bv,1);
632: PetscCheck(type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
634: if (bv->matrix) PetscCall(BVNorm_End_Private(bv,v,type,val)); /* non-standard inner product */
635: else PetscCall(VecNormEnd(v,type,val));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: /*@
640: BVNormColumn - Computes the vector norm of a selected column.
642: Collective
644: Input Parameters:
645: + bv - basis vectors
646: . j - column number to be used
647: - type - the norm type
649: Output Parameter:
650: . val - the norm
652: Notes:
653: The norm of V[j] is computed (NORM_1, NORM_2, or NORM_INFINITY).
654: If a non-standard inner product has been specified with BVSetMatrix(),
655: then the returned value is sqrt(V[j]'*B*V[j]),
656: where B is the inner product matrix (argument 'type' is ignored).
658: Level: intermediate
660: .seealso: BVNorm(), BVNormVec(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
661: @*/
662: PetscErrorCode BVNormColumn(BV bv,PetscInt j,NormType type,PetscReal *val)
663: {
664: Vec z;
666: PetscFunctionBegin;
670: PetscAssertPointer(val,4);
672: BVCheckSizes(bv,1);
674: PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
675: PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
677: PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
678: if (bv->matrix) { /* non-standard inner product */
679: PetscCall(BVGetColumn(bv,j,&z));
680: PetscCall(BVNorm_Private(bv,z,type,val));
681: PetscCall(BVRestoreColumn(bv,j,&z));
682: } else PetscUseTypeMethod(bv,norm,j,type,val);
683: PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@
688: BVNormColumnBegin - Starts a split phase norm computation.
690: Input Parameters:
691: + bv - basis vectors
692: . j - column number to be used
693: . type - the norm type
694: - val - the norm
696: Note:
697: Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
699: Level: advanced
701: .seealso: BVNormColumnEnd(), BVNormColumn()
702: @*/
703: PetscErrorCode BVNormColumnBegin(BV bv,PetscInt j,NormType type,PetscReal *val)
704: {
705: PetscSplitReduction *sr;
706: PetscReal lresult;
707: MPI_Comm comm;
708: Vec z;
710: PetscFunctionBegin;
714: PetscAssertPointer(val,4);
716: BVCheckSizes(bv,1);
718: PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
719: PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
721: PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
722: PetscCall(BVGetColumn(bv,j,&z));
723: if (bv->matrix) PetscCall(BVNorm_Begin_Private(bv,z,type,val)); /* non-standard inner product */
724: else if (bv->ops->norm_begin) PetscUseTypeMethod(bv,norm_begin,j,type,val);
725: else {
726: BVCheckOp(bv,1,norm_local);
727: PetscCall(PetscObjectGetComm((PetscObject)z,&comm));
728: PetscCall(PetscSplitReductionGet(comm,&sr));
729: PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
730: if (sr->numopsbegin >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
731: sr->invecs[sr->numopsbegin] = (void*)bv;
732: PetscUseTypeMethod(bv,norm_local,j,type,&lresult);
733: if (type == NORM_2) lresult = lresult*lresult;
734: if (type == NORM_MAX) sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_MAX;
735: else sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_SUM;
736: sr->lvalues[sr->numopsbegin++] = lresult;
737: }
738: PetscCall(BVRestoreColumn(bv,j,&z));
739: PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: BVNormColumnEnd - Ends a split phase norm computation.
746: Input Parameters:
747: + bv - basis vectors
748: . j - column number to be used
749: . type - the norm type
750: - val - the norm
752: Note:
753: Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
755: Level: advanced
757: .seealso: BVNormColumnBegin(), BVNormColumn()
758: @*/
759: PetscErrorCode BVNormColumnEnd(BV bv,PetscInt j,NormType type,PetscReal *val)
760: {
761: PetscSplitReduction *sr;
762: MPI_Comm comm;
763: Vec z;
765: PetscFunctionBegin;
769: PetscAssertPointer(val,4);
771: BVCheckSizes(bv,1);
773: PetscCheck(type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
775: PetscCall(BVGetColumn(bv,j,&z));
776: if (bv->matrix) PetscCall(BVNorm_End_Private(bv,z,type,val)); /* non-standard inner product */
777: else if (bv->ops->norm_end) PetscUseTypeMethod(bv,norm_end,j,type,val);
778: else {
779: PetscCall(PetscObjectGetComm((PetscObject)z,&comm));
780: PetscCall(PetscSplitReductionGet(comm,&sr));
781: PetscCall(PetscSplitReductionEnd(sr));
783: PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times then VecxxxBegin()");
784: PetscCheck((void*)bv==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() in a different order or with a different vector than VecxxxBegin()");
785: PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_MAX || type!=NORM_MAX,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called BVNormEnd(,NORM_MAX,) on a reduction started with VecDotBegin() or NORM_1 or NORM_2");
786: *val = PetscRealPart(sr->gvalues[sr->numopsend++]);
787: if (type == NORM_2) *val = PetscSqrtReal(*val);
788: if (sr->numopsend == sr->numopsbegin) {
789: sr->state = STATE_BEGIN;
790: sr->numopsend = 0;
791: sr->numopsbegin = 0;
792: }
793: }
794: PetscCall(BVRestoreColumn(bv,j,&z));
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*@
799: BVNormalize - Normalize all columns (starting from the leading ones).
801: Collective
803: Input Parameters:
804: + bv - basis vectors
805: - eigi - (optional) imaginary parts of eigenvalues
807: Notes:
808: On output, all columns will have unit norm. The normalization is done with
809: respect to the 2-norm, or to the B-norm if a non-standard inner product has
810: been specified with BVSetMatrix(), see BVNormColumn().
812: If the optional argument eigi is passed (taken into account only in real
813: scalars) it is interpreted as the imaginary parts of the eigenvalues and
814: the BV is supposed to contain the corresponding eigenvectors. Suppose the
815: first three values are eigi = { 0, alpha, -alpha }, then the first column
816: is normalized as usual, but the second and third ones are normalized assuming
817: that they contain the real and imaginary parts of a complex conjugate pair of
818: eigenvectors.
820: If eigi is passed, the inner-product matrix is ignored.
822: If there are leading columns, they are not modified (are assumed to be already
823: normalized).
825: Level: intermediate
827: .seealso: BVNormColumn()
828: @*/
829: PetscErrorCode BVNormalize(BV bv,PetscScalar *eigi)
830: {
831: PetscReal norm;
832: PetscInt i;
834: PetscFunctionBegin;
837: BVCheckSizes(bv,1);
839: PetscCall(PetscLogEventBegin(BV_Normalize,bv,0,0,0));
840: if (bv->matrix && !eigi) {
841: for (i=bv->l;i<bv->k;i++) {
842: PetscCall(BVNormColumn(bv,i,NORM_2,&norm));
843: PetscCall(BVScaleColumn(bv,i,1.0/norm));
844: }
845: } else PetscTryTypeMethod(bv,normalize,eigi);
846: PetscCall(PetscLogEventEnd(BV_Normalize,bv,0,0,0));
847: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: /*
852: Compute Y^H*A*X: right part column by column (with MatMult) and bottom
853: part row by row (with MatMultHermitianTranspose); result placed in marray[*,ldm]
854: */
855: static inline PetscErrorCode BVMatProject_Vec(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
856: {
857: PetscInt i,j,lx,ly,kx,ky,ulim;
858: Vec z,f;
860: PetscFunctionBegin;
861: lx = X->l; kx = X->k;
862: ly = Y->l; ky = Y->k;
863: PetscCall(BVCreateVec(X,&f));
864: BVCheckOp(Y,3,dotvec);
865: for (j=lx;j<kx;j++) {
866: PetscCall(BVGetColumn(X,j,&z));
867: PetscCall(MatMult(A,z,f));
868: PetscCall(BVRestoreColumn(X,j,&z));
869: ulim = PetscMin(ly+(j-lx)+1,ky);
870: Y->l = 0; Y->k = ulim;
871: PetscUseTypeMethod(Y,dotvec,f,marray+j*ldm);
872: if (symm) {
873: for (i=0;i<j;i++) marray[j+i*ldm] = PetscConj(marray[i+j*ldm]);
874: }
875: }
876: if (!symm) {
877: BVCheckOp(X,1,dotvec);
878: PetscCall(BV_AllocateCoeffs(Y));
879: for (j=ly;j<ky;j++) {
880: PetscCall(BVGetColumn(Y,j,&z));
881: PetscCall(MatMultHermitianTranspose(A,z,f));
882: PetscCall(BVRestoreColumn(Y,j,&z));
883: ulim = PetscMin(lx+(j-ly),kx);
884: X->l = 0; X->k = ulim;
885: PetscUseTypeMethod(X,dotvec,f,Y->h);
886: for (i=0;i<ulim;i++) marray[j+i*ldm] = PetscConj(Y->h[i]);
887: }
888: }
889: PetscCall(VecDestroy(&f));
890: X->l = lx; X->k = kx;
891: Y->l = ly; Y->k = ky;
892: PetscFunctionReturn(PETSC_SUCCESS);
893: }
895: /*
896: Compute Y^H*A*X= [ -- | Y0'*W1 ]
897: [ Y1'*W0 | Y1'*W1 ]
898: Allocates auxiliary BV to store the result of A*X, then one BVDot
899: call for top-right part and another one for bottom part;
900: result placed in marray[*,ldm]
901: */
902: static inline PetscErrorCode BVMatProject_MatMult(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm)
903: {
904: PetscInt j,lx,ly,kx,ky;
905: const PetscScalar *harray;
906: Mat H;
907: BV W;
909: PetscFunctionBegin;
910: lx = X->l; kx = X->k;
911: ly = Y->l; ky = Y->k;
912: PetscCall(BVDuplicate(X,&W));
913: X->l = 0; X->k = kx;
914: W->l = 0; W->k = kx;
915: PetscCall(BVMatMult(X,A,W));
917: /* top-right part, Y0'*AX1 */
918: if (ly>0 && lx<kx) {
919: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H));
920: W->l = lx; W->k = kx;
921: Y->l = 0; Y->k = ly;
922: PetscCall(BVDot(W,Y,H));
923: PetscCall(MatDenseGetArrayRead(H,&harray));
924: for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+j*ly,ly));
925: PetscCall(MatDenseRestoreArrayRead(H,&harray));
926: PetscCall(MatDestroy(&H));
927: }
929: /* bottom part, Y1'*AX */
930: if (kx>0 && ly<ky) {
931: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H));
932: W->l = 0; W->k = kx;
933: Y->l = ly; Y->k = ky;
934: PetscCall(BVDot(W,Y,H));
935: PetscCall(MatDenseGetArrayRead(H,&harray));
936: for (j=0;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly));
937: PetscCall(MatDenseRestoreArrayRead(H,&harray));
938: PetscCall(MatDestroy(&H));
939: }
940: PetscCall(BVDestroy(&W));
941: X->l = lx; X->k = kx;
942: Y->l = ly; Y->k = ky;
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*
947: Compute Y^H*A*X= [ -- | Y0'*W1 ]
948: [ Y1'*W0 | Y1'*W1 ]
949: First stage: allocate auxiliary BV to store A*X1, one BVDot for right part;
950: Second stage: resize BV to accommodate A'*Y1, then call BVDot for transpose of
951: bottom-left part; result placed in marray[*,ldm]
952: */
953: static inline PetscErrorCode BVMatProject_MatMult_2(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
954: {
955: PetscInt i,j,lx,ly,kx,ky;
956: const PetscScalar *harray;
957: Mat H;
958: BV W;
960: PetscFunctionBegin;
961: lx = X->l; kx = X->k;
962: ly = Y->l; ky = Y->k;
964: /* right part, Y'*AX1 */
965: PetscCall(BVDuplicateResize(X,kx-lx,&W));
966: if (ky>0 && lx<kx) {
967: PetscCall(BVMatMult(X,A,W));
968: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx-lx,NULL,&H));
969: Y->l = 0; Y->k = ky;
970: PetscCall(BVDot(W,Y,H));
971: PetscCall(MatDenseGetArrayRead(H,&harray));
972: for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+(j-lx)*ky,ky));
973: PetscCall(MatDenseRestoreArrayRead(H,&harray));
974: PetscCall(MatDestroy(&H));
975: }
977: /* bottom-left part, Y1'*AX0 */
978: if (lx>0 && ly<ky) {
979: if (symm) {
980: /* do not compute, just copy symmetric elements */
981: for (i=ly;i<ky;i++) {
982: for (j=0;j<lx;j++) marray[i+j*ldm] = PetscConj(marray[j+i*ldm]);
983: }
984: } else {
985: PetscCall(BVResize(W,ky-ly,PETSC_FALSE));
986: Y->l = ly; Y->k = ky;
987: PetscCall(BVMatMultHermitianTranspose(Y,A,W));
988: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,lx,ky-ly,NULL,&H));
989: X->l = 0; X->k = lx;
990: PetscCall(BVDot(W,X,H));
991: PetscCall(MatDenseGetArrayRead(H,&harray));
992: for (i=0;i<ky-ly;i++) {
993: for (j=0;j<lx;j++) {
994: marray[i+j*ldm+ly] = PetscConj(harray[j+i*(ky-ly)]);
995: }
996: }
997: PetscCall(MatDenseRestoreArrayRead(H,&harray));
998: PetscCall(MatDestroy(&H));
999: }
1000: }
1001: PetscCall(BVDestroy(&W));
1002: X->l = lx; X->k = kx;
1003: Y->l = ly; Y->k = ky;
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: /*
1008: Compute Y^H*X = [ -- | Y0'*X1 ] (X contains A*X):
1009: [ Y1'*X0 | Y1'*X1 ]
1010: one BVDot call for top-right part and another one for bottom part;
1011: result placed in marray[*,ldm]
1012: */
1013: static inline PetscErrorCode BVMatProject_Dot(BV X,BV Y,PetscScalar *marray,PetscInt ldm)
1014: {
1015: PetscInt j,lx,ly,kx,ky;
1016: const PetscScalar *harray;
1017: Mat H;
1019: PetscFunctionBegin;
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: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H));
1026: X->l = lx; X->k = kx;
1027: Y->l = 0; Y->k = ly;
1028: PetscCall(BVDot(X,Y,H));
1029: PetscCall(MatDenseGetArrayRead(H,&harray));
1030: for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+j*ly,ly));
1031: PetscCall(MatDenseRestoreArrayRead(H,&harray));
1032: PetscCall(MatDestroy(&H));
1033: }
1035: /* bottom part, Y1'*X */
1036: if (kx>0 && ly<ky) {
1037: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H));
1038: X->l = 0; X->k = kx;
1039: Y->l = ly; Y->k = ky;
1040: PetscCall(BVDot(X,Y,H));
1041: PetscCall(MatDenseGetArrayRead(H,&harray));
1042: for (j=0;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly));
1043: PetscCall(MatDenseRestoreArrayRead(H,&harray));
1044: PetscCall(MatDestroy(&H));
1045: }
1046: X->l = lx; X->k = kx;
1047: Y->l = ly; Y->k = ky;
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@
1052: BVMatProject - Computes the projection of a matrix onto a subspace.
1054: Collective
1056: Input Parameters:
1057: + X - basis vectors
1058: . A - (optional) matrix to be projected
1059: - Y - left basis vectors, can be equal to X
1061: Output Parameter:
1062: . M - the resulting matrix
1064: Notes:
1065: If A=NULL, then it is assumed that X already contains A*X.
1067: This operation is similar to BVDot(), with important differences.
1068: The goal is to compute the matrix resulting from the orthogonal projection
1069: of A onto the subspace spanned by the columns of X, M = X^H*A*X, or the
1070: oblique projection onto X along Y, M = Y^H*A*X.
1072: A difference with respect to BVDot() is that the standard inner product
1073: is always used, regardless of a non-standard inner product being specified
1074: with BVSetMatrix().
1076: On entry, M must be a sequential dense Mat with dimensions ky,kx at least,
1077: where ky (resp. kx) is the number of active columns of Y (resp. X).
1078: Another difference with respect to BVDot() is that all entries of M are
1079: computed except the leading ly,lx part, where ly (resp. lx) is the
1080: number of leading columns of Y (resp. X). Hence, the leading columns of
1081: X and Y participate in the computation, as opposed to BVDot().
1082: The leading part of M is assumed to be already available from previous
1083: computations.
1085: In the orthogonal projection case, Y=X, some computation can be saved if
1086: A is real symmetric (or complex Hermitian). In order to exploit this
1087: property, the symmetry flag of A must be set with MatSetOption().
1089: Level: intermediate
1091: .seealso: BVDot(), BVSetActiveColumns(), BVSetMatrix()
1092: @*/
1093: PetscErrorCode BVMatProject(BV X,Mat A,BV Y,Mat M)
1094: {
1095: PetscBool set,flg,symm=PETSC_FALSE;
1096: PetscInt m,n,ldm;
1097: PetscScalar *marray;
1098: Mat Xmatrix,Ymatrix;
1099: PetscObjectId idx,idy;
1101: PetscFunctionBegin;
1107: BVCheckSizes(X,1);
1108: if (A) {
1110: PetscCheckSameComm(X,1,A,2);
1111: }
1113: BVCheckSizes(Y,3);
1115: PetscCheckSameTypeAndComm(X,1,Y,3);
1116: SlepcMatCheckSeq(M);
1118: PetscCall(MatGetSize(M,&m,&n));
1119: PetscCall(MatDenseGetLDA(M,&ldm));
1120: PetscCheck(m>=Y->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,Y->k);
1121: PetscCheck(n>=X->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,X->k);
1122: PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
1124: PetscCall(PetscLogEventBegin(BV_MatProject,X,A,Y,0));
1125: /* temporarily set standard inner product */
1126: Xmatrix = X->matrix;
1127: Ymatrix = Y->matrix;
1128: X->matrix = Y->matrix = NULL;
1130: PetscCall(PetscObjectGetId((PetscObject)X,&idx));
1131: PetscCall(PetscObjectGetId((PetscObject)Y,&idy));
1132: if (A && idx==idy) { /* check symmetry of M=X'AX */
1133: PetscCall(MatIsHermitianKnown(A,&set,&flg));
1134: symm = set? flg: PETSC_FALSE;
1135: }
1137: PetscCall(MatDenseGetArray(M,&marray));
1139: if (A) {
1140: if (X->vmm==BV_MATMULT_VECS) {
1141: /* perform computation column by column */
1142: PetscCall(BVMatProject_Vec(X,A,Y,marray,ldm,symm));
1143: } else {
1144: /* use BVMatMult, then BVDot */
1145: PetscCall(MatHasOperation(A,MATOP_MULT_TRANSPOSE,&flg));
1146: if (symm || (flg && X->l>=X->k/2 && Y->l>=Y->k/2)) PetscCall(BVMatProject_MatMult_2(X,A,Y,marray,ldm,symm));
1147: else PetscCall(BVMatProject_MatMult(X,A,Y,marray,ldm));
1148: }
1149: } else {
1150: /* use BVDot on subblocks */
1151: PetscCall(BVMatProject_Dot(X,Y,marray,ldm));
1152: }
1154: PetscCall(MatDenseRestoreArray(M,&marray));
1155: PetscCall(PetscLogEventEnd(BV_MatProject,X,A,Y,0));
1156: /* restore non-standard inner product */
1157: X->matrix = Xmatrix;
1158: Y->matrix = Ymatrix;
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }