Actual source code: bvops.c
slepc-3.18.1 2022-11-02
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, except those involving global communication
12: */
14: #include <slepc/private/bvimpl.h>
15: #include <slepcds.h>
17: /*@
18: BVMult - Computes Y = beta*Y + alpha*X*Q.
20: Logically Collective on Y
22: Input Parameters:
23: + Y - first basis vectors context (modified on output)
24: . alpha - first scalar
25: . beta - second scalar
26: . X - second basis vectors context
27: - Q - (optional) sequential dense matrix
29: Notes:
30: X and Y must be different objects. The case X=Y can be addressed with
31: BVMultInPlace().
33: If matrix Q is NULL, then an AXPY operation Y = beta*Y + alpha*X is done
34: (i.e. results as if Q = identity). If provided,
35: the matrix Q must be a sequential dense Mat, with all entries equal on
36: all processes (otherwise each process will compute a different update).
37: The dimensions of Q must be at least m,n where m is the number of active
38: columns of X and n is the number of active columns of Y.
40: The leading columns of Y are not modified. Also, if X has leading
41: columns specified, then these columns do not participate in the computation.
42: Hence, only rows (resp. columns) of Q starting from lx (resp. ly) are used,
43: where lx (resp. ly) is the number of leading columns of X (resp. Y).
45: Level: intermediate
47: .seealso: BVMultVec(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
48: @*/
49: PetscErrorCode BVMult(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
50: {
51: PetscInt m,n;
59: BVCheckSizes(Y,1);
60: BVCheckOp(Y,1,mult);
62: BVCheckSizes(X,4);
66: if (Q) {
68: MatGetSize(Q,&m,&n);
71: }
74: PetscLogEventBegin(BV_Mult,X,Y,0,0);
75: PetscUseTypeMethod(Y,mult,alpha,beta,X,Q);
76: PetscLogEventEnd(BV_Mult,X,Y,0,0);
77: PetscObjectStateIncrease((PetscObject)Y);
78: return 0;
79: }
81: /*@
82: BVMultVec - Computes y = beta*y + alpha*X*q.
84: Logically Collective on X
86: Input Parameters:
87: + X - a basis vectors object
88: . alpha - first scalar
89: . beta - second scalar
90: . y - a vector (modified on output)
91: - q - an array of scalars
93: Notes:
94: This operation is the analogue of BVMult() but with a BV and a Vec,
95: instead of two BV. Note that arguments are listed in different order
96: with respect to BVMult().
98: If X has leading columns specified, then these columns do not participate
99: in the computation.
101: The length of array q must be equal to the number of active columns of X
102: minus the number of leading columns, i.e. the first entry of q multiplies
103: the first non-leading column.
105: Level: intermediate
107: .seealso: BVMult(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
108: @*/
109: PetscErrorCode BVMultVec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar q[])
110: {
111: PetscInt n,N;
119: BVCheckSizes(X,1);
120: BVCheckOp(X,1,multvec);
124: VecGetSize(y,&N);
125: VecGetLocalSize(y,&n);
128: PetscLogEventBegin(BV_MultVec,X,y,0,0);
129: PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
130: PetscLogEventEnd(BV_MultVec,X,y,0,0);
131: return 0;
132: }
134: /*@
135: BVMultColumn - Computes y = beta*y + alpha*X*q, where y is the j-th column
136: of X.
138: Logically Collective on X
140: Input Parameters:
141: + X - a basis vectors object
142: . alpha - first scalar
143: . beta - second scalar
144: . j - the column index
145: - q - an array of scalars
147: Notes:
148: This operation is equivalent to BVMultVec() but it uses column j of X
149: rather than taking a Vec as an argument. The number of active columns of
150: X is set to j before the computation, and restored afterwards.
151: If X has leading columns specified, then these columns do not participate
152: in the computation. Therefore, the length of array q must be equal to j
153: minus the number of leading columns.
155: Developer Notes:
156: If q is NULL, then the coefficients are taken from position nc+l of the
157: internal buffer vector, see BVGetBufferVec().
159: Level: advanced
161: .seealso: BVMult(), BVMultVec(), BVMultInPlace(), BVSetActiveColumns()
162: @*/
163: PetscErrorCode BVMultColumn(BV X,PetscScalar alpha,PetscScalar beta,PetscInt j,PetscScalar *q)
164: {
165: PetscInt ksave;
166: Vec y;
173: BVCheckSizes(X,1);
178: PetscLogEventBegin(BV_MultVec,X,0,0,0);
179: ksave = X->k;
180: X->k = j;
181: if (!q && !X->buffer) BVGetBufferVec(X,&X->buffer);
182: BVGetColumn(X,j,&y);
183: PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
184: BVRestoreColumn(X,j,&y);
185: X->k = ksave;
186: PetscLogEventEnd(BV_MultVec,X,0,0,0);
187: PetscObjectStateIncrease((PetscObject)X);
188: return 0;
189: }
191: /*@
192: BVMultInPlace - Update a set of vectors as V(:,s:e-1) = V*Q(:,s:e-1).
194: Logically Collective on V
196: Input Parameters:
197: + Q - a sequential dense matrix
198: . s - first column of V to be overwritten
199: - e - first column of V not to be overwritten
201: Input/Output Parameter:
202: . V - basis vectors
204: Notes:
205: The matrix Q must be a sequential dense Mat, with all entries equal on
206: all processes (otherwise each process will compute a different update).
208: This function computes V(:,s:e-1) = V*Q(:,s:e-1), that is, given a set of
209: vectors V, columns from s to e-1 are overwritten with columns from s to
210: e-1 of the matrix-matrix product V*Q. Only columns s to e-1 of Q are
211: referenced.
213: Level: intermediate
215: .seealso: BVMult(), BVMultVec(), BVMultInPlaceHermitianTranspose(), BVSetActiveColumns()
216: @*/
217: PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
218: {
219: PetscInt m,n;
226: BVCheckSizes(V,1);
232: MatGetSize(Q,&m,&n);
235: if (s>=e) return 0;
237: PetscLogEventBegin(BV_MultInPlace,V,Q,0,0);
238: PetscUseTypeMethod(V,multinplace,Q,s,e);
239: PetscLogEventEnd(BV_MultInPlace,V,Q,0,0);
240: PetscObjectStateIncrease((PetscObject)V);
241: return 0;
242: }
244: /*@
245: BVMultInPlaceHermitianTranspose - Update a set of vectors as V(:,s:e-1) = V*Q'(:,s:e-1).
247: Logically Collective on V
249: Input Parameters:
250: + Q - a sequential dense matrix
251: . s - first column of V to be overwritten
252: - e - first column of V not to be overwritten
254: Input/Output Parameter:
255: . V - basis vectors
257: Notes:
258: This is a variant of BVMultInPlace() where the conjugate transpose
259: of Q is used.
261: Level: intermediate
263: .seealso: BVMultInPlace()
264: @*/
265: PetscErrorCode BVMultInPlaceHermitianTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
266: {
267: PetscInt m,n;
274: BVCheckSizes(V,1);
280: MatGetSize(Q,&m,&n);
283: if (s>=e || !V->n) return 0;
285: PetscLogEventBegin(BV_MultInPlace,V,Q,0,0);
286: PetscUseTypeMethod(V,multinplacetrans,Q,s,e);
287: PetscLogEventEnd(BV_MultInPlace,V,Q,0,0);
288: PetscObjectStateIncrease((PetscObject)V);
289: return 0;
290: }
292: /*@
293: BVScale - Multiply the BV entries by a scalar value.
295: Logically Collective on bv
297: Input Parameters:
298: + bv - basis vectors
299: - alpha - scaling factor
301: Note:
302: All active columns (except the leading ones) are scaled.
304: Level: intermediate
306: .seealso: BVScaleColumn(), BVSetActiveColumns()
307: @*/
308: PetscErrorCode BVScale(BV bv,PetscScalar alpha)
309: {
313: BVCheckSizes(bv,1);
314: if (alpha == (PetscScalar)1.0) return 0;
316: PetscLogEventBegin(BV_Scale,bv,0,0,0);
317: if (bv->n) PetscUseTypeMethod(bv,scale,-1,alpha);
318: PetscLogEventEnd(BV_Scale,bv,0,0,0);
319: PetscObjectStateIncrease((PetscObject)bv);
320: return 0;
321: }
323: /*@
324: BVScaleColumn - Scale one column of a BV.
326: Logically Collective on bv
328: Input Parameters:
329: + bv - basis vectors
330: . j - column number to be scaled
331: - alpha - scaling factor
333: Level: intermediate
335: .seealso: BVScale(), BVSetActiveColumns()
336: @*/
337: PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
338: {
343: BVCheckSizes(bv,1);
346: if (alpha == (PetscScalar)1.0) return 0;
348: PetscLogEventBegin(BV_Scale,bv,0,0,0);
349: if (bv->n) PetscUseTypeMethod(bv,scale,j,alpha);
350: PetscLogEventEnd(BV_Scale,bv,0,0,0);
351: PetscObjectStateIncrease((PetscObject)bv);
352: return 0;
353: }
355: static inline PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
356: {
357: PetscInt i,low,high;
358: PetscScalar *px,t;
359: Vec x;
361: BVGetColumn(bv,k,&x);
362: if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
363: VecGetOwnershipRange(x,&low,&high);
364: VecGetArray(x,&px);
365: for (i=0;i<bv->N;i++) {
366: PetscRandomGetValue(bv->rand,&t);
367: if (i>=low && i<high) px[i-low] = t;
368: }
369: VecRestoreArray(x,&px);
370: } else VecSetRandom(x,bv->rand);
371: BVRestoreColumn(bv,k,&x);
372: return 0;
373: }
375: static inline PetscErrorCode BVSetRandomNormalColumn_Private(BV bv,PetscInt k,Vec w1,Vec w2)
376: {
377: PetscInt i,low,high;
378: PetscScalar *px,s,t;
379: Vec x;
381: BVGetColumn(bv,k,&x);
382: if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
383: VecGetOwnershipRange(x,&low,&high);
384: VecGetArray(x,&px);
385: for (i=0;i<bv->N;i++) {
386: PetscRandomGetValue(bv->rand,&s);
387: PetscRandomGetValue(bv->rand,&t);
388: if (i>=low && i<high) {
389: #if defined(PETSC_USE_COMPLEX)
390: px[i-low] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(t)),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(t)));
391: #else
392: px[i-low] = PetscSqrtReal(-2.0*PetscLogReal(s))*PetscCosReal(2.0*PETSC_PI*t);
393: #endif
394: }
395: }
396: VecRestoreArray(x,&px);
397: } else VecSetRandomNormal(x,bv->rand,w1,w2);
398: BVRestoreColumn(bv,k,&x);
399: return 0;
400: }
402: static inline PetscErrorCode BVSetRandomSignColumn_Private(BV bv,PetscInt k)
403: {
404: PetscInt i,low,high;
405: PetscScalar *px,t;
406: Vec x;
408: BVGetColumn(bv,k,&x);
409: VecGetOwnershipRange(x,&low,&high);
410: if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
411: VecGetArray(x,&px);
412: for (i=0;i<bv->N;i++) {
413: PetscRandomGetValue(bv->rand,&t);
414: if (i>=low && i<high) px[i-low] = (PetscRealPart(t)<0.5)? -1.0: 1.0;
415: }
416: VecRestoreArray(x,&px);
417: } else {
418: VecSetRandom(x,bv->rand);
419: VecGetArray(x,&px);
420: for (i=low;i<high;i++) {
421: px[i-low] = (PetscRealPart(px[i-low])<0.5)? -1.0: 1.0;
422: }
423: VecRestoreArray(x,&px);
424: }
425: BVRestoreColumn(bv,k,&x);
426: return 0;
427: }
429: /*@
430: BVSetRandom - Set the columns of a BV to random numbers.
432: Logically Collective on bv
434: Input Parameters:
435: . bv - basis vectors
437: Note:
438: All active columns (except the leading ones) are modified.
440: Level: advanced
442: .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetRandomSign(), BVSetRandomCond(), BVSetActiveColumns()
443: @*/
444: PetscErrorCode BVSetRandom(BV bv)
445: {
446: PetscInt k;
450: BVCheckSizes(bv,1);
452: BVGetRandomContext(bv,&bv->rand);
453: PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
454: for (k=bv->l;k<bv->k;k++) BVSetRandomColumn_Private(bv,k);
455: PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
456: PetscObjectStateIncrease((PetscObject)bv);
457: return 0;
458: }
460: /*@
461: BVSetRandomColumn - Set one column of a BV to random numbers.
463: Logically Collective on bv
465: Input Parameters:
466: + bv - basis vectors
467: - j - column number to be set
469: Level: advanced
471: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomCond()
472: @*/
473: PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
474: {
478: BVCheckSizes(bv,1);
481: BVGetRandomContext(bv,&bv->rand);
482: PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
483: BVSetRandomColumn_Private(bv,j);
484: PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
485: PetscObjectStateIncrease((PetscObject)bv);
486: return 0;
487: }
489: /*@
490: BVSetRandomNormal - Set the columns of a BV to random numbers with a normal
491: distribution.
493: Logically Collective on bv
495: Input Parameter:
496: . bv - basis vectors
498: Notes:
499: All active columns (except the leading ones) are modified.
501: Other functions such as BVSetRandom(), BVSetRandomColumn(), and BVSetRandomCond()
502: produce random numbers with a uniform distribution. This function returns values
503: that fit a normal distribution (Gaussian).
505: Developer Notes:
506: The current implementation obtains each of the columns by applying the Box-Muller
507: transform on two random vectors with uniformly distributed entries.
509: Level: advanced
511: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
512: @*/
513: PetscErrorCode BVSetRandomNormal(BV bv)
514: {
515: PetscInt k;
516: Vec w1=NULL,w2=NULL;
520: BVCheckSizes(bv,1);
522: BVGetRandomContext(bv,&bv->rand);
523: if (!bv->rrandom) {
524: BVCreateVec(bv,&w1);
525: BVCreateVec(bv,&w2);
526: }
527: PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
528: for (k=bv->l;k<bv->k;k++) BVSetRandomNormalColumn_Private(bv,k,w1,w2);
529: PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
530: if (!bv->rrandom) {
531: VecDestroy(&w1);
532: VecDestroy(&w2);
533: }
534: PetscObjectStateIncrease((PetscObject)bv);
535: return 0;
536: }
538: /*@
539: BVSetRandomSign - Set the entries of a BV to values 1 or -1 with equal
540: probability.
542: Logically Collective on bv
544: Input Parameter:
545: . bv - basis vectors
547: Notes:
548: All active columns (except the leading ones) are modified.
550: This function is used, e.g., in contour integral methods when estimating
551: the number of eigenvalues enclosed by the contour via an unbiased
552: estimator of tr(f(A)) [Bai et al., JCAM 1996].
554: Developer Notes:
555: The current implementation obtains random numbers and then replaces them
556: with -1 or 1 depending on the value being less than 0.5 or not.
558: Level: advanced
560: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetActiveColumns()
561: @*/
562: PetscErrorCode BVSetRandomSign(BV bv)
563: {
564: PetscScalar low,high;
565: PetscInt k;
569: BVCheckSizes(bv,1);
571: BVGetRandomContext(bv,&bv->rand);
572: PetscRandomGetInterval(bv->rand,&low,&high);
574: PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
575: for (k=bv->l;k<bv->k;k++) BVSetRandomSignColumn_Private(bv,k);
576: PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
577: PetscObjectStateIncrease((PetscObject)bv);
578: return 0;
579: }
581: /*@
582: BVSetRandomCond - Set the columns of a BV to random numbers, in a way that
583: the generated matrix has a given condition number.
585: Logically Collective on bv
587: Input Parameters:
588: + bv - basis vectors
589: - condn - condition number
591: Note:
592: All active columns (except the leading ones) are modified.
594: Level: advanced
596: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetActiveColumns()
597: @*/
598: PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
599: {
600: PetscInt k,i;
601: PetscScalar *eig,*d;
602: DS ds;
603: Mat A,X,Xt,M,G;
607: BVCheckSizes(bv,1);
609: BVGetRandomContext(bv,&bv->rand);
610: PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
611: /* B = rand(n,k) */
612: for (k=bv->l;k<bv->k;k++) BVSetRandomColumn_Private(bv,k);
613: DSCreate(PetscObjectComm((PetscObject)bv),&ds);
614: DSSetType(ds,DSHEP);
615: DSAllocate(ds,bv->m);
616: DSSetDimensions(ds,bv->k,bv->l,bv->k);
617: /* [V,S] = eig(B'*B) */
618: DSGetMat(ds,DS_MAT_A,&A);
619: BVDot(bv,bv,A);
620: DSRestoreMat(ds,DS_MAT_A,&A);
621: PetscMalloc1(bv->m,&eig);
622: DSSolve(ds,eig,NULL);
623: DSSynchronize(ds,eig,NULL);
624: DSVectors(ds,DS_MAT_X,NULL,NULL);
625: /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
626: MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M);
627: MatZeroEntries(M);
628: MatDenseGetArray(M,&d);
629: for (i=0;i<bv->k;i++) d[i+i*bv->m] = (1.0/condn+(1.0-1.0/condn)/(bv->k-1)*i)/PetscSqrtScalar(eig[i]);
630: MatDenseRestoreArray(M,&d);
631: /* G = X*M*X' */
632: DSGetMat(ds,DS_MAT_X,&X);
633: MatTranspose(X,MAT_INITIAL_MATRIX,&Xt);
634: MatProductCreate(Xt,M,NULL,&G);
635: MatProductSetType(G,MATPRODUCT_PtAP);
636: MatProductSetFromOptions(G);
637: MatProductSymbolic(G);
638: MatProductNumeric(G);
639: MatProductClear(G);
640: DSRestoreMat(ds,DS_MAT_X,&X);
641: MatDestroy(&Xt);
642: MatDestroy(&M);
643: /* B = B*G */
644: BVMultInPlace(bv,G,bv->l,bv->k);
645: MatDestroy(&G);
646: PetscFree(eig);
647: DSDestroy(&ds);
648: PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
649: PetscObjectStateIncrease((PetscObject)bv);
650: return 0;
651: }
653: /*@
654: BVMatMult - Computes the matrix-vector product for each column, Y=A*V.
656: Neighbor-wise Collective on A
658: Input Parameters:
659: + V - basis vectors context
660: - A - the matrix
662: Output Parameter:
663: . Y - the result
665: Notes:
666: Both V and Y must be distributed in the same manner. Only active columns
667: (excluding the leading ones) are processed.
668: In the result Y, columns are overwritten starting from the leading ones.
669: The number of active columns in V and Y should match, although they need
670: not be the same columns.
672: It is possible to choose whether the computation is done column by column
673: or as a Mat-Mat product, see BVSetMatMultMethod().
675: Level: beginner
677: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
678: @*/
679: PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
680: {
681: PetscInt M,N,m,n;
685: BVCheckSizes(V,1);
686: BVCheckOp(V,1,matmult);
691: BVCheckSizes(Y,3);
695: MatGetSize(A,&M,&N);
696: MatGetLocalSize(A,&m,&n);
703: PetscLogEventBegin(BV_MatMult,V,A,Y,0);
704: PetscUseTypeMethod(V,matmult,A,Y);
705: PetscLogEventEnd(BV_MatMult,V,A,Y,0);
706: PetscObjectStateIncrease((PetscObject)Y);
707: return 0;
708: }
710: /*@
711: BVMatMultTranspose - Computes the matrix-vector product with the transpose
712: of a matrix for each column, Y=A^T*V.
714: Neighbor-wise Collective on A
716: Input Parameters:
717: + V - basis vectors context
718: - A - the matrix
720: Output Parameter:
721: . Y - the result
723: Notes:
724: Both V and Y must be distributed in the same manner. Only active columns
725: (excluding the leading ones) are processed.
726: In the result Y, columns are overwritten starting from the leading ones.
727: The number of active columns in V and Y should match, although they need
728: not be the same columns.
730: Currently implemented via MatCreateTranspose().
732: Level: beginner
734: .seealso: BVMatMult(), BVMatMultHermitianTranspose()
735: @*/
736: PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
737: {
738: PetscInt M,N,m,n;
739: Mat AT;
743: BVCheckSizes(V,1);
748: BVCheckSizes(Y,3);
752: MatGetSize(A,&M,&N);
753: MatGetLocalSize(A,&m,&n);
760: MatCreateTranspose(A,&AT);
761: BVMatMult(V,AT,Y);
762: MatDestroy(&AT);
763: return 0;
764: }
766: /*@
767: BVMatMultHermitianTranspose - Computes the matrix-vector product with the
768: conjugate transpose of a matrix for each column, Y=A^H*V.
770: Neighbor-wise Collective on A
772: Input Parameters:
773: + V - basis vectors context
774: - A - the matrix
776: Output Parameter:
777: . Y - the result
779: Note:
780: Both V and Y must be distributed in the same manner. Only active columns
781: (excluding the leading ones) are processed.
782: In the result Y, columns are overwritten starting from the leading ones.
783: The number of active columns in V and Y should match, although they need
784: not be the same columns.
786: Currently implemented via MatCreateHermitianTranspose().
788: Level: beginner
790: .seealso: BVMatMult(), BVMatMultTranspose()
791: @*/
792: PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
793: {
794: PetscInt j,M,N,m,n;
795: Vec v,y;
799: BVCheckSizes(V,1);
804: BVCheckSizes(Y,3);
808: MatGetSize(A,&M,&N);
809: MatGetLocalSize(A,&m,&n);
816: /* TODO: recover this code if PETSc ever gets MATPRODUCT_AhB done
817: MatCreateHermitianTranspose(A,&AH);
818: BVMatMult(V,AH,Y);
819: MatDestroy(&AH);
820: */
821: for (j=0;j<V->k-V->l;j++) {
822: BVGetColumn(V,V->l+j,&v);
823: BVGetColumn(Y,Y->l+j,&y);
824: MatMultHermitianTranspose(A,v,y);
825: BVRestoreColumn(V,V->l+j,&v);
826: BVRestoreColumn(Y,Y->l+j,&y);
827: }
828: return 0;
829: }
831: /*@
832: BVMatMultColumn - Computes the matrix-vector product for a specified
833: column, storing the result in the next column v_{j+1}=A*v_j.
835: Neighbor-wise Collective on A
837: Input Parameters:
838: + V - basis vectors context
839: . A - the matrix
840: - j - the column
842: Level: beginner
844: .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
845: @*/
846: PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
847: {
848: Vec vj,vj1;
852: BVCheckSizes(V,1);
859: PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
860: BVGetColumn(V,j,&vj);
861: BVGetColumn(V,j+1,&vj1);
862: MatMult(A,vj,vj1);
863: BVRestoreColumn(V,j,&vj);
864: BVRestoreColumn(V,j+1,&vj1);
865: PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
866: PetscObjectStateIncrease((PetscObject)V);
867: return 0;
868: }
870: /*@
871: BVMatMultTransposeColumn - Computes the transpose matrix-vector product for a
872: specified column, storing the result in the next column v_{j+1}=A^T*v_j.
874: Neighbor-wise Collective on A
876: Input Parameters:
877: + V - basis vectors context
878: . A - the matrix
879: - j - the column
881: Level: beginner
883: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
884: @*/
885: PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
886: {
887: Vec vj,vj1;
891: BVCheckSizes(V,1);
898: PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
899: BVGetColumn(V,j,&vj);
900: BVGetColumn(V,j+1,&vj1);
901: MatMultTranspose(A,vj,vj1);
902: BVRestoreColumn(V,j,&vj);
903: BVRestoreColumn(V,j+1,&vj1);
904: PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
905: PetscObjectStateIncrease((PetscObject)V);
906: return 0;
907: }
909: /*@
910: BVMatMultHermitianTransposeColumn - Computes the conjugate-transpose matrix-vector
911: product for a specified column, storing the result in the next column v_{j+1}=A^H*v_j.
913: Neighbor-wise Collective on A
915: Input Parameters:
916: + V - basis vectors context
917: . A - the matrix
918: - j - the column
920: Level: beginner
922: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
923: @*/
924: PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
925: {
926: Vec vj,vj1;
930: BVCheckSizes(V,1);
937: PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
938: BVGetColumn(V,j,&vj);
939: BVGetColumn(V,j+1,&vj1);
940: MatMultHermitianTranspose(A,vj,vj1);
941: BVRestoreColumn(V,j,&vj);
942: BVRestoreColumn(V,j+1,&vj1);
943: PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
944: PetscObjectStateIncrease((PetscObject)V);
945: return 0;
946: }