Actual source code: bvops.c
slepc-3.21.0 2024-03-30
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
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;
53: PetscFunctionBegin;
60: BVCheckSizes(Y,1);
61: BVCheckOp(Y,1,mult);
63: BVCheckSizes(X,4);
65: PetscCheckSameTypeAndComm(Y,1,X,4);
66: PetscCheck(X!=Y,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONG,"X and Y arguments must be different");
67: if (Q) {
68: PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
69: PetscCall(MatGetSize(Q,&m,&n));
70: PetscCheck(m>=X->k,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,X->k);
71: PetscCheck(n>=Y->k,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,Y->k);
72: }
73: PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
75: PetscCall(PetscLogEventBegin(BV_Mult,X,Y,0,0));
76: PetscUseTypeMethod(Y,mult,alpha,beta,X,Q);
77: PetscCall(PetscLogEventEnd(BV_Mult,X,Y,0,0));
78: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*@
83: BVMultVec - Computes y = beta*y + alpha*X*q.
85: Logically Collective
87: Input Parameters:
88: + X - a basis vectors object
89: . alpha - first scalar
90: . beta - second scalar
91: . y - a vector (modified on output)
92: - q - an array of scalars
94: Notes:
95: This operation is the analogue of BVMult() but with a BV and a Vec,
96: instead of two BV. Note that arguments are listed in different order
97: with respect to BVMult().
99: If X has leading columns specified, then these columns do not participate
100: in the computation.
102: The length of array q must be equal to the number of active columns of X
103: minus the number of leading columns, i.e. the first entry of q multiplies
104: the first non-leading column.
106: Level: intermediate
108: .seealso: BVMult(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
109: @*/
110: PetscErrorCode BVMultVec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar q[])
111: {
112: PetscInt n,N;
114: PetscFunctionBegin;
119: PetscAssertPointer(q,5);
121: BVCheckSizes(X,1);
122: BVCheckOp(X,1,multvec);
124: PetscCheckSameComm(X,1,y,4);
126: PetscCall(VecGetSize(y,&N));
127: PetscCall(VecGetLocalSize(y,&n));
128: PetscCheck(N==X->N && n==X->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ") do not match BV sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ")",N,n,X->N,X->n);
130: PetscCall(PetscLogEventBegin(BV_MultVec,X,y,0,0));
131: PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
132: PetscCall(PetscLogEventEnd(BV_MultVec,X,y,0,0));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*@
137: BVMultColumn - Computes y = beta*y + alpha*X*q, where y is the j-th column
138: of X.
140: Logically Collective
142: Input Parameters:
143: + X - a basis vectors object
144: . alpha - first scalar
145: . beta - second scalar
146: . j - the column index
147: - q - an array of scalars
149: Notes:
150: This operation is equivalent to BVMultVec() but it uses column j of X
151: rather than taking a Vec as an argument. The number of active columns of
152: X is set to j before the computation, and restored afterwards.
153: If X has leading columns specified, then these columns do not participate
154: in the computation. Therefore, the length of array q must be equal to j
155: minus the number of leading columns.
157: Developer Notes:
158: If q is NULL, then the coefficients are taken from position nc+l of the
159: internal buffer vector, see BVGetBufferVec().
161: Level: advanced
163: .seealso: BVMult(), BVMultVec(), BVMultInPlace(), BVSetActiveColumns()
164: @*/
165: PetscErrorCode BVMultColumn(BV X,PetscScalar alpha,PetscScalar beta,PetscInt j,PetscScalar *q)
166: {
167: PetscInt ksave;
168: Vec y;
170: PetscFunctionBegin;
176: BVCheckSizes(X,1);
178: PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
179: 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);
181: PetscCall(PetscLogEventBegin(BV_MultVec,X,0,0,0));
182: ksave = X->k;
183: X->k = j;
184: if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
185: PetscCall(BVGetColumn(X,j,&y));
186: PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
187: PetscCall(BVRestoreColumn(X,j,&y));
188: X->k = ksave;
189: PetscCall(PetscLogEventEnd(BV_MultVec,X,0,0,0));
190: PetscCall(PetscObjectStateIncrease((PetscObject)X));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*@
195: BVMultInPlace - Update a set of vectors as V(:,s:e-1) = V*Q(:,s:e-1).
197: Logically Collective
199: Input Parameters:
200: + Q - a sequential dense matrix
201: . s - first column of V to be overwritten
202: - e - first column of V not to be overwritten
204: Input/Output Parameter:
205: . V - basis vectors
207: Notes:
208: The matrix Q must be a sequential dense Mat, with all entries equal on
209: all processes (otherwise each process will compute a different update).
211: This function computes V(:,s:e-1) = V*Q(:,s:e-1), that is, given a set of
212: vectors V, columns from s to e-1 are overwritten with columns from s to
213: e-1 of the matrix-matrix product V*Q. Only columns s to e-1 of Q are
214: referenced.
216: Level: intermediate
218: .seealso: BVMult(), BVMultVec(), BVMultInPlaceHermitianTranspose(), BVSetActiveColumns()
219: @*/
220: PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
221: {
222: PetscInt m,n;
224: PetscFunctionBegin;
230: BVCheckSizes(V,1);
232: PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
234: PetscCheck(s>=V->l && s<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,s,V->l,V->m);
235: PetscCheck(e>=V->l && e<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,e,V->l,V->m);
236: PetscCall(MatGetSize(Q,&m,&n));
237: PetscCheck(m>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,V->k);
238: PetscCheck(e<=n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %" PetscInt_FMT " columns, the requested value of e is larger: %" PetscInt_FMT,n,e);
240: PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
241: PetscUseTypeMethod(V,multinplace,Q,s,e);
242: PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
243: PetscCall(PetscObjectStateIncrease((PetscObject)V));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@
248: BVMultInPlaceHermitianTranspose - Update a set of vectors as V(:,s:e-1) = V*Q'(:,s:e-1).
250: Logically Collective
252: Input Parameters:
253: + Q - a sequential dense matrix
254: . s - first column of V to be overwritten
255: - e - first column of V not to be overwritten
257: Input/Output Parameter:
258: . V - basis vectors
260: Notes:
261: This is a variant of BVMultInPlace() where the conjugate transpose
262: of Q is used.
264: Level: intermediate
266: .seealso: BVMultInPlace()
267: @*/
268: PetscErrorCode BVMultInPlaceHermitianTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
269: {
270: PetscInt m,n;
272: PetscFunctionBegin;
278: BVCheckSizes(V,1);
280: PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
282: PetscCheck(s>=V->l && s<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,s,V->l,V->m);
283: PetscCheck(e>=V->l && e<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,e,V->l,V->m);
284: PetscCall(MatGetSize(Q,&m,&n));
285: PetscCheck(n>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,V->k);
286: PetscCheck(e<=m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %" PetscInt_FMT " rows, the requested value of e is larger: %" PetscInt_FMT,m,e);
288: PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
289: PetscUseTypeMethod(V,multinplacetrans,Q,s,e);
290: PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
291: PetscCall(PetscObjectStateIncrease((PetscObject)V));
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@
296: BVScale - Multiply the BV entries by a scalar value.
298: Logically Collective
300: Input Parameters:
301: + bv - basis vectors
302: - alpha - scaling factor
304: Note:
305: All active columns (except the leading ones) are scaled.
307: Level: intermediate
309: .seealso: BVScaleColumn(), BVSetActiveColumns()
310: @*/
311: PetscErrorCode BVScale(BV bv,PetscScalar alpha)
312: {
313: PetscFunctionBegin;
317: BVCheckSizes(bv,1);
318: if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
320: PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
321: PetscUseTypeMethod(bv,scale,-1,alpha);
322: PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
323: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /*@
328: BVScaleColumn - Scale one column of a BV.
330: Logically Collective
332: Input Parameters:
333: + bv - basis vectors
334: . j - column number to be scaled
335: - alpha - scaling factor
337: Level: intermediate
339: .seealso: BVScale(), BVSetActiveColumns()
340: @*/
341: PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
342: {
343: PetscFunctionBegin;
348: BVCheckSizes(bv,1);
350: 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);
351: if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
353: PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
354: PetscUseTypeMethod(bv,scale,j,alpha);
355: PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
356: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static inline PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
361: {
362: PetscInt i,low,high;
363: PetscScalar *px,t;
364: Vec x;
366: PetscFunctionBegin;
367: PetscCall(BVGetColumn(bv,k,&x));
368: if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
369: PetscCall(VecGetOwnershipRange(x,&low,&high));
370: PetscCall(VecGetArray(x,&px));
371: for (i=0;i<bv->N;i++) {
372: PetscCall(PetscRandomGetValue(bv->rand,&t));
373: if (i>=low && i<high) px[i-low] = t;
374: }
375: PetscCall(VecRestoreArray(x,&px));
376: } else PetscCall(VecSetRandom(x,bv->rand));
377: PetscCall(BVRestoreColumn(bv,k,&x));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static inline PetscErrorCode BVSetRandomNormalColumn_Private(BV bv,PetscInt k,Vec w1,Vec w2)
382: {
383: PetscInt i,low,high;
384: PetscScalar *px,s,t;
385: Vec x;
387: PetscFunctionBegin;
388: PetscCall(BVGetColumn(bv,k,&x));
389: if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
390: PetscCall(VecGetOwnershipRange(x,&low,&high));
391: PetscCall(VecGetArray(x,&px));
392: for (i=0;i<bv->N;i++) {
393: PetscCall(PetscRandomGetValue(bv->rand,&s));
394: PetscCall(PetscRandomGetValue(bv->rand,&t));
395: if (i>=low && i<high) {
396: #if defined(PETSC_USE_COMPLEX)
397: 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)));
398: #else
399: px[i-low] = PetscSqrtReal(-2.0*PetscLogReal(s))*PetscCosReal(2.0*PETSC_PI*t);
400: #endif
401: }
402: }
403: PetscCall(VecRestoreArray(x,&px));
404: } else PetscCall(VecSetRandomNormal(x,bv->rand,w1,w2));
405: PetscCall(BVRestoreColumn(bv,k,&x));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static inline PetscErrorCode BVSetRandomSignColumn_Private(BV bv,PetscInt k)
410: {
411: PetscInt i,low,high;
412: PetscScalar *px,t;
413: Vec x;
415: PetscFunctionBegin;
416: PetscCall(BVGetColumn(bv,k,&x));
417: PetscCall(VecGetOwnershipRange(x,&low,&high));
418: if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
419: PetscCall(VecGetArray(x,&px));
420: for (i=0;i<bv->N;i++) {
421: PetscCall(PetscRandomGetValue(bv->rand,&t));
422: if (i>=low && i<high) px[i-low] = (PetscRealPart(t)<0.5)? -1.0: 1.0;
423: }
424: PetscCall(VecRestoreArray(x,&px));
425: } else {
426: PetscCall(VecSetRandom(x,bv->rand));
427: PetscCall(VecGetArray(x,&px));
428: for (i=low;i<high;i++) {
429: px[i-low] = (PetscRealPart(px[i-low])<0.5)? -1.0: 1.0;
430: }
431: PetscCall(VecRestoreArray(x,&px));
432: }
433: PetscCall(BVRestoreColumn(bv,k,&x));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: /*@
438: BVSetRandom - Set the columns of a BV to random numbers.
440: Logically Collective
442: Input Parameters:
443: . bv - basis vectors
445: Note:
446: All active columns (except the leading ones) are modified.
448: Level: advanced
450: .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetRandomSign(), BVSetRandomCond(), BVSetActiveColumns()
451: @*/
452: PetscErrorCode BVSetRandom(BV bv)
453: {
454: PetscInt k;
456: PetscFunctionBegin;
459: BVCheckSizes(bv,1);
461: PetscCall(BVGetRandomContext(bv,&bv->rand));
462: PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
463: for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
464: PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
465: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: BVSetRandomColumn - Set one column of a BV to random numbers.
472: Logically Collective
474: Input Parameters:
475: + bv - basis vectors
476: - j - column number to be set
478: Level: advanced
480: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomCond()
481: @*/
482: PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
483: {
484: PetscFunctionBegin;
488: BVCheckSizes(bv,1);
489: 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);
491: PetscCall(BVGetRandomContext(bv,&bv->rand));
492: PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
493: PetscCall(BVSetRandomColumn_Private(bv,j));
494: PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
495: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*@
500: BVSetRandomNormal - Set the columns of a BV to random numbers with a normal
501: distribution.
503: Logically Collective
505: Input Parameter:
506: . bv - basis vectors
508: Notes:
509: All active columns (except the leading ones) are modified.
511: Other functions such as BVSetRandom(), BVSetRandomColumn(), and BVSetRandomCond()
512: produce random numbers with a uniform distribution. This function returns values
513: that fit a normal distribution (Gaussian).
515: Developer Notes:
516: The current implementation obtains each of the columns by applying the Box-Muller
517: transform on two random vectors with uniformly distributed entries.
519: Level: advanced
521: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
522: @*/
523: PetscErrorCode BVSetRandomNormal(BV bv)
524: {
525: PetscInt k;
526: Vec w1=NULL,w2=NULL;
528: PetscFunctionBegin;
531: BVCheckSizes(bv,1);
533: PetscCall(BVGetRandomContext(bv,&bv->rand));
534: if (!bv->rrandom) {
535: PetscCall(BVCreateVec(bv,&w1));
536: PetscCall(BVCreateVec(bv,&w2));
537: }
538: PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
539: for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomNormalColumn_Private(bv,k,w1,w2));
540: PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
541: if (!bv->rrandom) {
542: PetscCall(VecDestroy(&w1));
543: PetscCall(VecDestroy(&w2));
544: }
545: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: /*@
550: BVSetRandomSign - Set the entries of a BV to values 1 or -1 with equal
551: probability.
553: Logically Collective
555: Input Parameter:
556: . bv - basis vectors
558: Notes:
559: All active columns (except the leading ones) are modified.
561: This function is used, e.g., in contour integral methods when estimating
562: the number of eigenvalues enclosed by the contour via an unbiased
563: estimator of tr(f(A)) [Bai et al., JCAM 1996].
565: Developer Notes:
566: The current implementation obtains random numbers and then replaces them
567: with -1 or 1 depending on the value being less than 0.5 or not.
569: Level: advanced
571: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetActiveColumns()
572: @*/
573: PetscErrorCode BVSetRandomSign(BV bv)
574: {
575: PetscScalar low,high;
576: PetscInt k;
578: PetscFunctionBegin;
581: BVCheckSizes(bv,1);
583: PetscCall(BVGetRandomContext(bv,&bv->rand));
584: PetscCall(PetscRandomGetInterval(bv->rand,&low,&high));
585: PetscCheck(PetscRealPart(low)==0.0 && PetscRealPart(high)==1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The PetscRandom object in the BV must have interval [0,1]");
586: PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
587: for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomSignColumn_Private(bv,k));
588: PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
589: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@
594: BVSetRandomCond - Set the columns of a BV to random numbers, in a way that
595: the generated matrix has a given condition number.
597: Logically Collective
599: Input Parameters:
600: + bv - basis vectors
601: - condn - condition number
603: Note:
604: All active columns (except the leading ones) are modified.
606: Level: advanced
608: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetActiveColumns()
609: @*/
610: PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
611: {
612: PetscInt k,i;
613: PetscScalar *eig,*d;
614: DS ds;
615: Mat A,X,Xt,M,G;
617: PetscFunctionBegin;
620: BVCheckSizes(bv,1);
622: PetscCall(BVGetRandomContext(bv,&bv->rand));
623: PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
624: /* B = rand(n,k) */
625: for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
626: PetscCall(DSCreate(PetscObjectComm((PetscObject)bv),&ds));
627: PetscCall(DSSetType(ds,DSHEP));
628: PetscCall(DSAllocate(ds,bv->m));
629: PetscCall(DSSetDimensions(ds,bv->k,bv->l,bv->k));
630: /* [V,S] = eig(B'*B) */
631: PetscCall(DSGetMat(ds,DS_MAT_A,&A));
632: PetscCall(BVDot(bv,bv,A));
633: PetscCall(DSRestoreMat(ds,DS_MAT_A,&A));
634: PetscCall(PetscMalloc1(bv->m,&eig));
635: PetscCall(DSSolve(ds,eig,NULL));
636: PetscCall(DSSynchronize(ds,eig,NULL));
637: PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
638: /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
639: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M));
640: PetscCall(MatZeroEntries(M));
641: PetscCall(MatDenseGetArray(M,&d));
642: 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]);
643: PetscCall(MatDenseRestoreArray(M,&d));
644: /* G = X*M*X' */
645: PetscCall(DSGetMat(ds,DS_MAT_X,&X));
646: PetscCall(MatTranspose(X,MAT_INITIAL_MATRIX,&Xt));
647: PetscCall(MatProductCreate(Xt,M,NULL,&G));
648: PetscCall(MatProductSetType(G,MATPRODUCT_PtAP));
649: PetscCall(MatProductSetFromOptions(G));
650: PetscCall(MatProductSymbolic(G));
651: PetscCall(MatProductNumeric(G));
652: PetscCall(MatProductClear(G));
653: PetscCall(DSRestoreMat(ds,DS_MAT_X,&X));
654: PetscCall(MatDestroy(&Xt));
655: PetscCall(MatDestroy(&M));
656: /* B = B*G */
657: PetscCall(BVMultInPlace(bv,G,bv->l,bv->k));
658: PetscCall(MatDestroy(&G));
659: PetscCall(PetscFree(eig));
660: PetscCall(DSDestroy(&ds));
661: PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
662: PetscCall(PetscObjectStateIncrease((PetscObject)bv));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: /*@
667: BVMatMult - Computes the matrix-vector product for each column, Y=A*V.
669: Neighbor-wise Collective
671: Input Parameters:
672: + V - basis vectors context
673: - A - the matrix
675: Output Parameter:
676: . Y - the result
678: Notes:
679: Both V and Y must be distributed in the same manner. Only active columns
680: (excluding the leading ones) are processed.
681: In the result Y, columns are overwritten starting from the leading ones.
682: The number of active columns in V and Y should match, although they need
683: not be the same columns.
685: It is possible to choose whether the computation is done column by column
686: or as a Mat-Mat product, see BVSetMatMultMethod().
688: Level: beginner
690: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
691: @*/
692: PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
693: {
694: PetscInt M,N,m,n;
696: PetscFunctionBegin;
699: BVCheckSizes(V,1);
700: BVCheckOp(V,1,matmult);
705: BVCheckSizes(Y,3);
706: PetscCheckSameComm(V,1,A,2);
707: PetscCheckSameTypeAndComm(V,1,Y,3);
709: PetscCall(MatGetSize(A,&M,&N));
710: PetscCall(MatGetLocalSize(A,&m,&n));
711: PetscCheck(M==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,M,Y->N);
712: PetscCheck(m==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,m,Y->n);
713: PetscCheck(N==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,N,V->N);
714: PetscCheck(n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,n,V->n);
715: PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);
717: PetscCall(PetscLogEventBegin(BV_MatMult,V,A,Y,0));
718: PetscUseTypeMethod(V,matmult,A,Y);
719: PetscCall(PetscLogEventEnd(BV_MatMult,V,A,Y,0));
720: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: BVMatMultTranspose - Computes the matrix-vector product with the transpose
726: of a matrix for each column, Y=A^T*V.
728: Neighbor-wise Collective
730: Input Parameters:
731: + V - basis vectors context
732: - A - the matrix
734: Output Parameter:
735: . Y - the result
737: Notes:
738: Both V and Y must be distributed in the same manner. Only active columns
739: (excluding the leading ones) are processed.
740: In the result Y, columns are overwritten starting from the leading ones.
741: The number of active columns in V and Y should match, although they need
742: not be the same columns.
744: Currently implemented via MatCreateTranspose().
746: Level: beginner
748: .seealso: BVMatMult(), BVMatMultHermitianTranspose()
749: @*/
750: PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
751: {
752: PetscInt M,N,m,n;
753: Mat AT;
755: PetscFunctionBegin;
758: BVCheckSizes(V,1);
763: BVCheckSizes(Y,3);
764: PetscCheckSameComm(V,1,A,2);
765: PetscCheckSameTypeAndComm(V,1,Y,3);
767: PetscCall(MatGetSize(A,&M,&N));
768: PetscCall(MatGetLocalSize(A,&m,&n));
769: PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
770: PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
771: PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
772: PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
773: PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);
775: PetscCall(MatCreateTranspose(A,&AT));
776: PetscCall(BVMatMult(V,AT,Y));
777: PetscCall(MatDestroy(&AT));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@
782: BVMatMultHermitianTranspose - Computes the matrix-vector product with the
783: conjugate transpose of a matrix for each column, Y=A^H*V.
785: Neighbor-wise Collective
787: Input Parameters:
788: + V - basis vectors context
789: - A - the matrix
791: Output Parameter:
792: . Y - the result
794: Note:
795: Both V and Y must be distributed in the same manner. Only active columns
796: (excluding the leading ones) are processed.
797: In the result Y, columns are overwritten starting from the leading ones.
798: The number of active columns in V and Y should match, although they need
799: not be the same columns.
801: Currently implemented via MatCreateHermitianTranspose().
803: Level: beginner
805: .seealso: BVMatMult(), BVMatMultTranspose()
806: @*/
807: PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
808: {
809: PetscInt j,M,N,m,n;
810: Vec v,y;
812: PetscFunctionBegin;
815: BVCheckSizes(V,1);
820: BVCheckSizes(Y,3);
821: PetscCheckSameComm(V,1,A,2);
822: PetscCheckSameTypeAndComm(V,1,Y,3);
824: PetscCall(MatGetSize(A,&M,&N));
825: PetscCall(MatGetLocalSize(A,&m,&n));
826: PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
827: PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
828: PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
829: PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
830: PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);
832: /* TODO: recover this code if PETSc ever gets MATPRODUCT_AhB done
833: PetscCall(MatCreateHermitianTranspose(A,&AH));
834: PetscCall(BVMatMult(V,AH,Y));
835: PetscCall(MatDestroy(&AH));
836: */
837: for (j=0;j<V->k-V->l;j++) {
838: PetscCall(BVGetColumn(V,V->l+j,&v));
839: PetscCall(BVGetColumn(Y,Y->l+j,&y));
840: PetscCall(MatMultHermitianTranspose(A,v,y));
841: PetscCall(BVRestoreColumn(V,V->l+j,&v));
842: PetscCall(BVRestoreColumn(Y,Y->l+j,&y));
843: }
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*@
848: BVMatMultColumn - Computes the matrix-vector product for a specified
849: column, storing the result in the next column v_{j+1}=A*v_j.
851: Neighbor-wise Collective
853: Input Parameters:
854: + V - basis vectors context
855: . A - the matrix
856: - j - the column
858: Level: beginner
860: .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
861: @*/
862: PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
863: {
864: Vec vj,vj1;
866: PetscFunctionBegin;
869: BVCheckSizes(V,1);
871: PetscCheckSameComm(V,1,A,2);
873: PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
874: PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);
876: PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
877: PetscCall(BVGetColumn(V,j,&vj));
878: PetscCall(BVGetColumn(V,j+1,&vj1));
879: PetscCall(MatMult(A,vj,vj1));
880: PetscCall(BVRestoreColumn(V,j,&vj));
881: PetscCall(BVRestoreColumn(V,j+1,&vj1));
882: PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
883: PetscCall(PetscObjectStateIncrease((PetscObject)V));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: /*@
888: BVMatMultTransposeColumn - Computes the transpose matrix-vector product for a
889: specified column, storing the result in the next column v_{j+1}=A^T*v_j.
891: Neighbor-wise Collective
893: Input Parameters:
894: + V - basis vectors context
895: . A - the matrix
896: - j - the column
898: Level: beginner
900: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
901: @*/
902: PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
903: {
904: Vec vj,vj1;
906: PetscFunctionBegin;
909: BVCheckSizes(V,1);
911: PetscCheckSameComm(V,1,A,2);
913: PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
914: PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);
916: PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
917: PetscCall(BVGetColumn(V,j,&vj));
918: PetscCall(BVGetColumn(V,j+1,&vj1));
919: PetscCall(MatMultTranspose(A,vj,vj1));
920: PetscCall(BVRestoreColumn(V,j,&vj));
921: PetscCall(BVRestoreColumn(V,j+1,&vj1));
922: PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
923: PetscCall(PetscObjectStateIncrease((PetscObject)V));
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: /*@
928: BVMatMultHermitianTransposeColumn - Computes the conjugate-transpose matrix-vector
929: product for a specified column, storing the result in the next column v_{j+1}=A^H*v_j.
931: Neighbor-wise Collective
933: Input Parameters:
934: + V - basis vectors context
935: . A - the matrix
936: - j - the column
938: Level: beginner
940: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
941: @*/
942: PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
943: {
944: Vec vj,vj1;
946: PetscFunctionBegin;
949: BVCheckSizes(V,1);
951: PetscCheckSameComm(V,1,A,2);
953: PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
954: PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);
956: PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
957: PetscCall(BVGetColumn(V,j,&vj));
958: PetscCall(BVGetColumn(V,j+1,&vj1));
959: PetscCall(MatMultHermitianTranspose(A,vj,vj1));
960: PetscCall(BVRestoreColumn(V,j,&vj));
961: PetscCall(BVRestoreColumn(V,j+1,&vj1));
962: PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
963: PetscCall(PetscObjectStateIncrease((PetscObject)V));
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }