Actual source code: bvimpl.h
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: */
11: #pragma once
13: #include <slepcbv.h>
14: #include <slepc/private/slepcimpl.h>
16: /* SUBMANSEC = BV */
18: SLEPC_EXTERN PetscBool BVRegisterAllCalled;
19: SLEPC_EXTERN PetscErrorCode BVRegisterAll(void);
21: SLEPC_EXTERN PetscLogEvent BV_Create,BV_Copy,BV_Mult,BV_MultVec,BV_MultInPlace,BV_Dot,BV_DotVec,BV_Orthogonalize,BV_OrthogonalizeVec,BV_Scale,BV_Norm,BV_NormVec,BV_Normalize,BV_SetRandom,BV_MatMult,BV_MatMultVec,BV_MatProject,BV_SVDAndRank;
23: typedef struct _BVOps *BVOps;
25: struct _BVOps {
26: PetscErrorCode (*mult)(BV,PetscScalar,PetscScalar,BV,Mat);
27: PetscErrorCode (*multvec)(BV,PetscScalar,PetscScalar,Vec,PetscScalar*);
28: PetscErrorCode (*multinplace)(BV,Mat,PetscInt,PetscInt);
29: PetscErrorCode (*multinplacetrans)(BV,Mat,PetscInt,PetscInt);
30: PetscErrorCode (*dot)(BV,BV,Mat);
31: PetscErrorCode (*dotvec)(BV,Vec,PetscScalar*);
32: PetscErrorCode (*dotvec_local)(BV,Vec,PetscScalar*);
33: PetscErrorCode (*dotvec_begin)(BV,Vec,PetscScalar*);
34: PetscErrorCode (*dotvec_end)(BV,Vec,PetscScalar*);
35: PetscErrorCode (*scale)(BV,PetscInt,PetscScalar);
36: PetscErrorCode (*norm)(BV,PetscInt,NormType,PetscReal*);
37: PetscErrorCode (*norm_local)(BV,PetscInt,NormType,PetscReal*);
38: PetscErrorCode (*norm_begin)(BV,PetscInt,NormType,PetscReal*);
39: PetscErrorCode (*norm_end)(BV,PetscInt,NormType,PetscReal*);
40: PetscErrorCode (*normalize)(BV,PetscScalar*);
41: PetscErrorCode (*matmult)(BV,Mat,BV);
42: PetscErrorCode (*copy)(BV,BV);
43: PetscErrorCode (*copycolumn)(BV,PetscInt,PetscInt);
44: PetscErrorCode (*resize)(BV,PetscInt,PetscBool);
45: PetscErrorCode (*getcolumn)(BV,PetscInt,Vec*);
46: PetscErrorCode (*restorecolumn)(BV,PetscInt,Vec*);
47: PetscErrorCode (*getarray)(BV,PetscScalar**);
48: PetscErrorCode (*restorearray)(BV,PetscScalar**);
49: PetscErrorCode (*getarrayread)(BV,const PetscScalar**);
50: PetscErrorCode (*restorearrayread)(BV,const PetscScalar**);
51: PetscErrorCode (*restoresplit)(BV,BV*,BV*);
52: PetscErrorCode (*restoresplitrows)(BV,IS,IS,BV*,BV*);
53: PetscErrorCode (*gramschmidt)(BV,PetscInt,Vec,PetscBool*,PetscScalar*,PetscScalar*,PetscReal*,PetscReal*);
54: PetscErrorCode (*getmat)(BV,Mat*);
55: PetscErrorCode (*restoremat)(BV,Mat*);
56: PetscErrorCode (*duplicate)(BV,BV);
57: PetscErrorCode (*create)(BV);
58: PetscErrorCode (*setfromoptions)(BV,PetscOptionItems);
59: PetscErrorCode (*view)(BV,PetscViewer);
60: PetscErrorCode (*destroy)(BV);
61: };
63: struct _p_BV {
64: PETSCHEADER(struct _BVOps);
65: /*------------------------- User parameters --------------------------*/
66: PetscLayout map; /* layout of columns */
67: VecType vtype; /* vector type */
68: PetscInt n,N; /* dimensions of vectors (local, global) */
69: PetscInt m; /* number of vectors */
70: PetscInt l; /* number of leading columns */
71: PetscInt k; /* number of active columns */
72: PetscInt nc; /* number of constraints */
73: PetscInt ld; /* leading dimension */
74: BVOrthogType orthog_type; /* the method of vector orthogonalization */
75: BVOrthogRefineType orthog_ref; /* refinement method */
76: PetscReal orthog_eta; /* refinement threshold */
77: BVOrthogBlockType orthog_block; /* the method of block orthogonalization */
78: Mat matrix; /* inner product matrix */
79: PetscBool indef; /* matrix is indefinite */
80: BVMatMultType vmm; /* version of matmult operation */
81: PetscBool rrandom; /* reproducible random vectors */
82: PetscReal deftol; /* tolerance for BV_SafeSqrt */
84: /*---------------------- Cached data and workspace -------------------*/
85: Vec buffer; /* buffer vector used in orthogonalization */
86: Mat Abuffer; /* auxiliary seqdense matrix that wraps the buffer */
87: Vec Bx; /* result of matrix times a vector x */
88: PetscObjectId xid; /* object id of vector x */
89: PetscObjectState xstate; /* state of vector x */
90: Vec cv[2]; /* column vectors obtained with BVGetColumn() */
91: PetscInt ci[2]; /* column indices of obtained vectors */
92: PetscObjectState st[2]; /* state of obtained vectors */
93: PetscObjectId id[2]; /* object id of obtained vectors */
94: PetscScalar *h,*c; /* orthogonalization coefficients */
95: Vec omega; /* signature matrix values for indefinite case */
96: PetscBool defersfo; /* deferred call to setfromoptions */
97: BV cached; /* cached BV to store result of matrix times BV */
98: PetscObjectState bvstate; /* state of BV when BVApplyMatrixBV() was called */
99: BV L,R; /* BV objects obtained with BVGetSplit/Rows() */
100: PetscObjectState lstate,rstate;/* state of L and R when BVGetSplit/Rows() was called */
101: PetscInt lsplit; /* value of l when BVGetSplit() was called (-1 if BVGetSplitRows()) */
102: PetscInt issplit; /* !=0 if BV is from split (1=left, 2=right, -1=top, -2=bottom) */
103: BV splitparent; /* my parent if I am a split BV */
104: PetscRandom rand; /* random number generator */
105: Mat Acreate; /* matrix given at BVCreateFromMat() */
106: Mat Aget; /* matrix returned for BVGetMat() */
107: PetscBool cuda; /* true if NVIDIA GPU must be used */
108: PetscBool hip; /* true if AMD GPU must be used */
109: PetscBool kokkos; /* true if Kokkos vectors are being used */
110: PetscBool sfocalled; /* setfromoptions has been called */
111: PetscScalar *work;
112: PetscInt lwork;
113: void *data;
114: };
116: /*
117: BV_SafeSqrt - Computes the square root of a scalar value alpha, which is
118: assumed to be z'*B*z. The result is
119: if definite inner product: res = sqrt(alpha)
120: if indefinite inner product: res = sgn(alpha)*sqrt(abs(alpha))
121: */
122: static inline PetscErrorCode BV_SafeSqrt(BV bv,PetscScalar alpha,PetscReal *res)
123: {
124: PetscReal absal,realp;
125: const char *msg;
127: PetscFunctionBegin;
128: absal = PetscAbsScalar(alpha);
129: realp = PetscRealPart(alpha);
130: if (PetscUnlikely(absal<PETSC_MACHINE_EPSILON)) PetscCall(PetscInfo(bv,"Zero norm %g, either the vector is zero or a semi-inner product is being used\n",(double)absal));
131: #if PetscDefined(USE_COMPLEX)
132: PetscCheck(PetscAbsReal(PetscImaginaryPart(alpha))<bv->deftol || PetscAbsReal(PetscImaginaryPart(alpha))/absal<10*bv->deftol,PetscObjectComm((PetscObject)bv),PETSC_ERR_USER_INPUT,"The inner product is not well defined: nonzero imaginary part %g",(double)PetscImaginaryPart(alpha));
133: #endif
134: if (PetscUnlikely(bv->indef)) {
135: *res = (realp<0.0)? -PetscSqrtReal(-realp): PetscSqrtReal(realp);
136: } else {
137: msg = bv->matrix? "The inner product is not well defined: indefinite matrix %g": "Invalid inner product: %g";
138: PetscCheck(realp>-bv->deftol,PetscObjectComm((PetscObject)bv),PETSC_ERR_USER_INPUT,msg,(double)realp);
139: *res = (realp<0.0)? 0.0: PetscSqrtReal(realp);
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /*
145: BV_IPMatMult - Multiply a vector x by the inner-product matrix, cache the
146: result in Bx.
147: */
148: static inline PetscErrorCode BV_IPMatMult(BV bv,Vec x)
149: {
150: PetscFunctionBegin;
151: if (((PetscObject)x)->id != bv->xid || ((PetscObject)x)->state != bv->xstate) {
152: if (PetscUnlikely(!bv->Bx)) PetscCall(MatCreateVecs(bv->matrix,&bv->Bx,NULL));
153: PetscCall(MatMult(bv->matrix,x,bv->Bx));
154: PetscCall(PetscObjectGetId((PetscObject)x,&bv->xid));
155: PetscCall(VecGetState(x,&bv->xstate));
156: }
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: /*
161: BV_IPMatMultBV - Multiply BV by the inner-product matrix, cache the
162: result internally in bv->cached.
163: */
164: static inline PetscErrorCode BV_IPMatMultBV(BV bv)
165: {
166: PetscFunctionBegin;
167: PetscCall(BVGetCachedBV(bv,&bv->cached));
168: if (((PetscObject)bv)->state != bv->bvstate || bv->l != bv->cached->l || bv->k != bv->cached->k) {
169: PetscCall(BVSetActiveColumns(bv->cached,bv->l,bv->k));
170: if (bv->matrix) PetscCall(BVMatMult(bv,bv->matrix,bv->cached));
171: else PetscCall(BVCopy(bv,bv->cached));
172: bv->bvstate = ((PetscObject)bv)->state;
173: }
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*
178: BV_AllocateCoeffs - Allocate orthogonalization coefficients if not done already.
179: */
180: static inline PetscErrorCode BV_AllocateCoeffs(BV bv)
181: {
182: PetscFunctionBegin;
183: if (!bv->h) PetscCall(PetscMalloc2(bv->nc+bv->m,&bv->h,bv->nc+bv->m,&bv->c));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /*
188: BV_AllocateSignature - Allocate signature coefficients if not done already.
189: */
190: static inline PetscErrorCode BV_AllocateSignature(BV bv)
191: {
192: PetscFunctionBegin;
193: if (bv->indef && !bv->omega) {
194: if (bv->cuda) {
195: #if PetscDefined(HAVE_CUDA)
196: PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,bv->nc+bv->m,&bv->omega));
197: #else
198: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
199: #endif
200: } else if (bv->hip) {
201: #if PetscDefined(HAVE_HIP)
202: PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,bv->nc+bv->m,&bv->omega));
203: #else
204: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
205: #endif
206: } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,bv->nc+bv->m,&bv->omega));
207: PetscCall(VecSet(bv->omega,1.0));
208: }
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*
213: BV_SetMatrixDiagonal - sets the inner product matrix for BV as a diagonal matrix
214: with the diagonal specified by vector vomega, using the same matrix type as matrix M
215: */
216: static inline PetscErrorCode BV_SetMatrixDiagonal(BV bv,Vec vomega,Mat M)
217: {
218: Mat Omega;
219: MatType Mtype;
221: PetscFunctionBegin;
222: PetscCall(MatGetType(M,&Mtype));
223: PetscCall(MatCreate(PetscObjectComm((PetscObject)bv),&Omega));
224: PetscCall(MatSetSizes(Omega,bv->n,bv->n,bv->N,bv->N));
225: PetscCall(MatSetType(Omega,Mtype));
226: PetscCall(MatDiagonalSet(Omega,vomega,INSERT_VALUES));
227: PetscCall(BVSetMatrix(bv,Omega,PETSC_TRUE));
228: PetscCall(MatDestroy(&Omega));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*
233: BVAvailableVec: First (0) or second (1) vector available for
234: getcolumn operation (or -1 if both vectors already fetched).
235: */
236: #define BVAvailableVec (((bv->ci[0]==-bv->nc-1)? 0: (bv->ci[1]==-bv->nc-1)? 1: -1))
238: /*
239: Macros to test valid BV arguments
240: */
241: #if !PetscDefined(USE_DEBUG)
243: #define BVCheckSizes(h,arg) do {(void)(h);} while (0)
244: #define BVCheckOp(h,arg,op) do {(void)(h);} while (0)
246: #else
248: #define BVCheckSizes(h,arg) \
249: do { \
250: PetscCheck((h)->m,PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"BV sizes have not been defined: Parameter #%d",arg); \
251: } while (0)
253: #define BVCheckOp(h,arg,op) \
254: do { \
255: PetscCheck((h)->ops->op,PetscObjectComm((PetscObject)(h)),PETSC_ERR_SUP,"Operation not implemented in this BV type: Parameter #%d",arg); \
256: } while (0)
258: #endif
260: SLEPC_INTERN PetscErrorCode BVView_Vecs(BV,PetscViewer);
262: SLEPC_INTERN PetscErrorCode BVAllocateWork_Private(BV,PetscInt);
264: SLEPC_INTERN PetscErrorCode BVMult_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscScalar,PetscScalar*,PetscInt);
265: SLEPC_INTERN PetscErrorCode BVMultVec_BLAS_Private(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,const PetscScalar*,PetscScalar,PetscScalar*);
266: SLEPC_INTERN PetscErrorCode BVMultInPlace_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscBool);
267: SLEPC_INTERN PetscErrorCode BVMultInPlace_Vecs_Private(BV,PetscInt,PetscInt,PetscInt,Vec*,const PetscScalar*,PetscBool);
268: SLEPC_INTERN PetscErrorCode BVAXPY_BLAS_Private(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,PetscScalar,PetscScalar*,PetscInt);
269: SLEPC_INTERN PetscErrorCode BVDot_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,const PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscScalar*,PetscInt,PetscBool);
270: SLEPC_INTERN PetscErrorCode BVDotVec_BLAS_Private(BV,PetscInt,PetscInt,const PetscScalar*,PetscInt,const PetscScalar*,PetscScalar*,PetscBool);
271: SLEPC_INTERN PetscErrorCode BVScale_BLAS_Private(BV,PetscInt,PetscScalar*,PetscScalar);
272: SLEPC_INTERN PetscErrorCode BVNorm_LAPACK_Private(BV,PetscInt,PetscInt,const PetscScalar*,PetscInt,NormType,PetscReal*,PetscBool);
273: SLEPC_INTERN PetscErrorCode BVNormalize_LAPACK_Private(BV,PetscInt,PetscInt,const PetscScalar*,PetscInt,PetscScalar*,PetscBool);
274: SLEPC_INTERN PetscErrorCode BVGetMat_Default(BV,Mat*);
275: SLEPC_INTERN PetscErrorCode BVRestoreMat_Default(BV,Mat*);
276: SLEPC_INTERN PetscErrorCode BVMatCholInv_LAPACK_Private(BV,Mat,Mat);
277: SLEPC_INTERN PetscErrorCode BVMatTriInv_LAPACK_Private(BV,Mat,Mat);
278: SLEPC_INTERN PetscErrorCode BVMatSVQB_LAPACK_Private(BV,Mat,Mat);
279: SLEPC_INTERN PetscErrorCode BVOrthogonalize_LAPACK_TSQR(BV,PetscInt,PetscInt,PetscScalar*,PetscInt,PetscScalar*,PetscInt);
280: SLEPC_INTERN PetscErrorCode BVOrthogonalize_LAPACK_TSQR_OnlyR(BV,PetscInt,PetscInt,PetscScalar*,PetscInt,PetscScalar*,PetscInt);
282: /* reduction operations used in BVOrthogonalize and BVNormalize */
283: SLEPC_EXTERN MPI_Op MPIU_TSQR, MPIU_LAPY2;
284: SLEPC_EXTERN void MPIAPI SlepcGivensPacked(void*,void*,PetscMPIInt*,MPI_Datatype*);
285: SLEPC_EXTERN void MPIAPI SlepcPythag(void*,void*,PetscMPIInt*,MPI_Datatype*);
287: /*
288: BV_CleanCoefficients_Default - Sets to zero all entries of column j of the bv buffer
289: */
290: static inline PetscErrorCode BV_CleanCoefficients_Default(BV bv,PetscInt j,PetscScalar *h)
291: {
292: PetscScalar *hh=h,*a;
293: PetscInt i;
295: PetscFunctionBegin;
296: if (!h) {
297: PetscCall(VecGetArray(bv->buffer,&a));
298: hh = a + j*(bv->nc+bv->m);
299: }
300: for (i=0;i<bv->nc+j;i++) hh[i] = 0.0;
301: if (!h) PetscCall(VecRestoreArray(bv->buffer,&a));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*
306: BV_AddCoefficients_Default - Add the contents of the scratch (0-th column) of the bv buffer
307: into column j of the bv buffer
308: */
309: static inline PetscErrorCode BV_AddCoefficients_Default(BV bv,PetscInt j,PetscScalar *h,PetscScalar *c)
310: {
311: PetscScalar *hh=h,*cc=c;
312: PetscInt i;
314: PetscFunctionBegin;
315: if (!h) {
316: PetscCall(VecGetArray(bv->buffer,&cc));
317: hh = cc + j*(bv->nc+bv->m);
318: }
319: for (i=0;i<bv->nc+j;i++) hh[i] += cc[i];
320: if (!h) PetscCall(VecRestoreArray(bv->buffer,&cc));
321: PetscCall(PetscLogFlops(1.0*(bv->nc+j)));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*
326: BV_SetValue_Default - Sets value in row j (counted after the constraints) of column k
327: of the coefficients array
328: */
329: static inline PetscErrorCode BV_SetValue_Default(BV bv,PetscInt j,PetscInt k,PetscScalar *h,PetscScalar value)
330: {
331: PetscScalar *hh=h,*a;
333: PetscFunctionBegin;
334: if (!h) {
335: PetscCall(VecGetArray(bv->buffer,&a));
336: hh = a + k*(bv->nc+bv->m);
337: }
338: hh[bv->nc+j] = value;
339: if (!h) PetscCall(VecRestoreArray(bv->buffer,&a));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*
344: BV_SquareSum_Default - Returns the value h'*h, where h represents the contents of the
345: coefficients array (up to position j)
346: */
347: static inline PetscErrorCode BV_SquareSum_Default(BV bv,PetscInt j,PetscScalar *h,PetscReal *sum)
348: {
349: PetscScalar *hh=h;
350: PetscInt i;
352: PetscFunctionBegin;
353: *sum = 0.0;
354: if (!h) PetscCall(VecGetArray(bv->buffer,&hh));
355: for (i=0;i<bv->nc+j;i++) *sum += PetscRealPart(hh[i]*PetscConj(hh[i]));
356: if (!h) PetscCall(VecRestoreArray(bv->buffer,&hh));
357: PetscCall(PetscLogFlops(2.0*(bv->nc+j)));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*
362: BV_ApplySignature_Default - Computes the pointwise product h*omega, where h represents
363: the contents of the coefficients array (up to position j) and omega is the signature;
364: if inverse=TRUE then the operation is h/omega
365: */
366: static inline PetscErrorCode BV_ApplySignature_Default(BV bv,PetscInt j,PetscScalar *h,PetscBool inverse)
367: {
368: PetscScalar *hh=h;
369: PetscInt i;
370: const PetscScalar *omega;
372: PetscFunctionBegin;
373: if (PetscUnlikely(!(bv->nc+j))) PetscFunctionReturn(PETSC_SUCCESS);
374: if (!h) PetscCall(VecGetArray(bv->buffer,&hh));
375: PetscCall(VecGetArrayRead(bv->omega,&omega));
376: if (inverse) for (i=0;i<bv->nc+j;i++) hh[i] /= PetscRealPart(omega[i]);
377: else for (i=0;i<bv->nc+j;i++) hh[i] *= PetscRealPart(omega[i]);
378: PetscCall(VecRestoreArrayRead(bv->omega,&omega));
379: if (!h) PetscCall(VecRestoreArray(bv->buffer,&hh));
380: PetscCall(PetscLogFlops(1.0*(bv->nc+j)));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*
385: BV_SquareRoot_Default - Returns the square root of position j (counted after the constraints)
386: of the coefficients array
387: */
388: static inline PetscErrorCode BV_SquareRoot_Default(BV bv,PetscInt j,PetscScalar *h,PetscReal *beta)
389: {
390: PetscScalar *hh=h;
392: PetscFunctionBegin;
393: if (!h) PetscCall(VecGetArray(bv->buffer,&hh));
394: PetscCall(BV_SafeSqrt(bv,hh[bv->nc+j],beta));
395: if (!h) PetscCall(VecRestoreArray(bv->buffer,&hh));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: /*
400: BV_StoreCoefficients_Default - Copy the contents of the coefficients array to an array dest
401: provided by the caller (only values from l to j are copied)
402: */
403: static inline PetscErrorCode BV_StoreCoefficients_Default(BV bv,PetscInt j,PetscScalar *h,PetscScalar *dest)
404: {
405: PetscScalar *hh=h,*a;
406: PetscInt i;
408: PetscFunctionBegin;
409: if (!h) {
410: PetscCall(VecGetArray(bv->buffer,&a));
411: hh = a + j*(bv->nc+bv->m);
412: }
413: for (i=bv->l;i<j;i++) dest[i-bv->l] = hh[bv->nc+i];
414: if (!h) PetscCall(VecRestoreArray(bv->buffer,&a));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /*
419: BV_GetEigenvector - retrieves k-th eigenvector from basis vectors V.
420: The argument eigi is the imaginary part of the corresponding eigenvalue.
421: */
422: static inline PetscErrorCode BV_GetEigenvector(BV V,PetscInt k,PetscScalar eigi,Vec Vr,Vec Vi)
423: {
424: PetscFunctionBegin;
425: #if PetscDefined(USE_COMPLEX)
426: (void)eigi;
427: if (Vr) PetscCall(BVCopyVec(V,k,Vr));
428: if (Vi) PetscCall(VecSet(Vi,0.0));
429: #else
430: if (eigi > 0.0) { /* first value of conjugate pair */
431: if (Vr) PetscCall(BVCopyVec(V,k,Vr));
432: if (Vi) PetscCall(BVCopyVec(V,k+1,Vi));
433: } else if (eigi < 0.0) { /* second value of conjugate pair */
434: if (Vr) PetscCall(BVCopyVec(V,k-1,Vr));
435: if (Vi) {
436: PetscCall(BVCopyVec(V,k,Vi));
437: PetscCall(VecScale(Vi,-1.0));
438: }
439: } else { /* real eigenvalue */
440: if (Vr) PetscCall(BVCopyVec(V,k,Vr));
441: if (Vi) PetscCall(VecSet(Vi,0.0));
442: }
443: #endif
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /*
448: BV_OrthogonalizeColumn_Safe - this is intended for cases where we know that
449: the resulting vector is going to be numerically zero, so normalization or
450: iterative refinement may cause problems in parallel (collective operation
451: not being called by all processes)
452: */
453: static inline PetscErrorCode BV_OrthogonalizeColumn_Safe(BV bv,PetscInt j,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
454: {
455: BVOrthogRefineType orthog_ref;
457: PetscFunctionBegin;
458: PetscCall(PetscInfo(bv,"Orthogonalizing column %" PetscInt_FMT " without refinement\n",j));
459: orthog_ref = bv->orthog_ref;
460: bv->orthog_ref = BV_ORTHOG_REFINE_NEVER; /* avoid refinement */
461: PetscCall(BVOrthogonalizeColumn(bv,j,H,NULL,NULL));
462: bv->orthog_ref = orthog_ref; /* restore refinement setting */
463: if (norm) *norm = 0.0;
464: if (lindep) *lindep = PETSC_TRUE;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*
469: BV_SetDefaultLD - set the default value of the leading dimension, based on
470: the local size.
471: */
472: static inline PetscErrorCode BV_SetDefaultLD(BV bv,PetscInt nloc)
473: {
474: size_t bytes,align;
476: PetscFunctionBegin;
477: if (bv->ld) { /* set by user */
478: PetscCheck(bv->ld>=nloc,PetscObjectComm((PetscObject)bv),PETSC_ERR_USER_INPUT,"The leading dimension %" PetscInt_FMT " should be larger or equal to the local number of rows %" PetscInt_FMT,bv->ld,nloc);
479: } else {
480: align = PetscMax(PETSC_MEMALIGN,16); /* assume that CUDA requires 16-byte alignment */
481: bytes = (nloc*sizeof(PetscScalar) + align - 1) & ~(align - 1);
482: PetscCall(PetscIntCast(bytes/sizeof(PetscScalar),&bv->ld));
483: }
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: #if PetscDefined(HAVE_CUDA)
488: /*
489: BV_MatDenseCUDAGetArrayRead - if Q is MATSEQDENSE it will allocate memory on the
490: GPU and copy the contents; otherwise, calls MatDenseCUDAGetArrayRead()
491: */
492: static inline PetscErrorCode BV_MatDenseCUDAGetArrayRead(PETSC_UNUSED BV bv,Mat Q,const PetscScalar **d_q)
493: {
494: const PetscScalar *q;
495: PetscInt ldq,mq;
496: PetscCuBLASInt ldq_=0;
497: PetscBool matiscuda;
499: PetscFunctionBegin;
500: PetscCall(MatGetSize(Q,NULL,&mq));
501: PetscCall(MatDenseGetLDA(Q,&ldq));
502: PetscCall(PetscCuBLASIntCast(ldq,&ldq_));
503: PetscCall(PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSECUDA,&matiscuda));
504: if (matiscuda) PetscCall(MatDenseCUDAGetArrayRead(Q,d_q));
505: else {
506: PetscCall(MatDenseGetArrayRead(Q,&q));
507: PetscCallCUDA(cudaMalloc((void**)d_q,ldq*mq*sizeof(PetscScalar)));
508: PetscCallCUDA(cudaMemcpy((void*)*d_q,q,ldq*mq*sizeof(PetscScalar),cudaMemcpyHostToDevice));
509: PetscCall(PetscLogCpuToGpu(ldq*mq*sizeof(PetscScalar)));
510: }
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /*
515: BV_MatDenseCUDARestoreArrayRead - restores the pointer obtained with BV_MatDenseCUDAGetArrayRead(),
516: freeing the GPU memory in case of MATSEQDENSE
517: */
518: static inline PetscErrorCode BV_MatDenseCUDARestoreArrayRead(PETSC_UNUSED BV bv,Mat Q,const PetscScalar **d_q)
519: {
520: PetscBool matiscuda;
522: PetscFunctionBegin;
523: PetscCall(PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSECUDA,&matiscuda));
524: if (matiscuda) PetscCall(MatDenseCUDARestoreArrayRead(Q,d_q));
525: else {
526: PetscCall(MatDenseRestoreArrayRead(Q,NULL));
527: PetscCallCUDA(cudaFree((void*)*d_q));
528: *d_q = NULL;
529: }
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*
534: BV_VecPlaceArray - allows the user to replace the device array in a vector with a
535: device array provided by the caller, calling the accessor that corresponds to the
536: vector type being used by the BV object
537: */
538: static inline PetscErrorCode BV_VecPlaceArray(PETSC_UNUSED BV bv,Vec v,PetscScalar *a)
539: {
540: PetscFunctionBegin;
541: #if PetscDefined(HAVE_KOKKOS_KERNELS)
542: if (bv->kokkos) {
543: PetscCall(VecKokkosPlaceArray(v,a));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
546: #endif
547: PetscCall(VecCUDAPlaceArray(v,a));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: /*
552: BV_VecResetArray - resets the device array in a vector to the one it had before the
553: call to BV_VecPlaceArray()
554: */
555: static inline PetscErrorCode BV_VecResetArray(PETSC_UNUSED BV bv,Vec v)
556: {
557: PetscFunctionBegin;
558: #if PetscDefined(HAVE_KOKKOS_KERNELS)
559: if (bv->kokkos) {
560: PetscCall(VecKokkosResetArray(v));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
563: #endif
564: PetscCall(VecCUDAResetArray(v));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: SLEPC_INTERN PetscErrorCode BVMult_BLAS_CUDA(BV,PetscInt,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscScalar,PetscScalar*,PetscInt);
569: SLEPC_INTERN PetscErrorCode BVMultVec_BLAS_CUDA(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,const PetscScalar*,PetscScalar,PetscScalar*);
570: SLEPC_INTERN PetscErrorCode BVMultInPlace_BLAS_CUDA(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscBool);
571: SLEPC_INTERN PetscErrorCode BVAXPY_BLAS_CUDA(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,PetscScalar,PetscScalar*,PetscInt);
572: SLEPC_INTERN PetscErrorCode BVDot_BLAS_CUDA(BV,PetscInt,PetscInt,PetscInt,const PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscScalar*,PetscInt,PetscBool);
573: SLEPC_INTERN PetscErrorCode BVDotVec_BLAS_CUDA(BV,PetscInt,PetscInt,const PetscScalar*,PetscInt,const PetscScalar*,PetscScalar*,PetscBool);
574: SLEPC_INTERN PetscErrorCode BVScale_BLAS_CUDA(BV,PetscInt,PetscScalar*,PetscScalar);
575: SLEPC_INTERN PetscErrorCode BVNorm_BLAS_CUDA(BV,PetscInt,const PetscScalar*,PetscReal*);
576: SLEPC_INTERN PetscErrorCode BVNormalize_BLAS_CUDA(BV,PetscInt,PetscInt,PetscScalar*,PetscInt,PetscScalar*);
578: SLEPC_INTERN PetscErrorCode BV_CleanCoefficients_CUDA(BV,PetscInt,PetscScalar*);
579: SLEPC_INTERN PetscErrorCode BV_AddCoefficients_CUDA(BV,PetscInt,PetscScalar*,PetscScalar*);
580: SLEPC_INTERN PetscErrorCode BV_SetValue_CUDA(BV,PetscInt,PetscInt,PetscScalar*,PetscScalar);
581: SLEPC_INTERN PetscErrorCode BV_SquareSum_CUDA(BV,PetscInt,PetscScalar*,PetscReal*);
582: SLEPC_INTERN PetscErrorCode BV_ApplySignature_CUDA(BV,PetscInt,PetscScalar*,PetscBool);
583: SLEPC_INTERN PetscErrorCode BV_SquareRoot_CUDA(BV,PetscInt,PetscScalar*,PetscReal*);
584: SLEPC_INTERN PetscErrorCode BV_StoreCoefficients_CUDA(BV,PetscInt,PetscScalar*,PetscScalar*);
585: #define BV_CleanCoefficients(a,b,c) ((a)->cuda?BV_CleanCoefficients_CUDA:BV_CleanCoefficients_Default)((a),(b),(c))
586: #define BV_AddCoefficients(a,b,c,d) ((a)->cuda?BV_AddCoefficients_CUDA:BV_AddCoefficients_Default)((a),(b),(c),(d))
587: #define BV_SetValue(a,b,c,d,e) ((a)->cuda?BV_SetValue_CUDA:BV_SetValue_Default)((a),(b),(c),(d),(e))
588: #define BV_SquareSum(a,b,c,d) ((a)->cuda?BV_SquareSum_CUDA:BV_SquareSum_Default)((a),(b),(c),(d))
589: #define BV_ApplySignature(a,b,c,d) ((a)->cuda?BV_ApplySignature_CUDA:BV_ApplySignature_Default)((a),(b),(c),(d))
590: #define BV_SquareRoot(a,b,c,d) ((a)->cuda?BV_SquareRoot_CUDA:BV_SquareRoot_Default)((a),(b),(c),(d))
591: #define BV_StoreCoefficients(a,b,c,d) ((a)->cuda?BV_StoreCoefficients_CUDA:BV_StoreCoefficients_Default)((a),(b),(c),(d))
593: #elif PetscDefined(HAVE_HIP)
594: #include <petscdevice_cupm.h>
595: /*
596: BV_MatDenseHIPGetArrayRead - if Q is MATSEQDENSE it will allocate memory on the
597: GPU and copy the contents; otherwise, calls MatDenseHIPGetArrayRead()
598: */
599: static inline PetscErrorCode BV_MatDenseHIPGetArrayRead(PETSC_UNUSED BV bv,Mat Q,const PetscScalar **d_q)
600: {
601: const PetscScalar *q;
602: PetscInt ldq,mq;
603: PetscCuBLASInt ldq_=0;
604: PetscBool matiship;
606: PetscFunctionBegin;
607: PetscCall(MatGetSize(Q,NULL,&mq));
608: PetscCall(MatDenseGetLDA(Q,&ldq));
609: PetscCall(PetscHipBLASIntCast(ldq,&ldq_));
610: PetscCall(PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSEHIP,&matiship));
611: if (matiship) PetscCall(MatDenseHIPGetArrayRead(Q,d_q));
612: else {
613: PetscCall(MatDenseGetArrayRead(Q,&q));
614: PetscCallHIP(hipMalloc((void**)d_q,ldq*mq*sizeof(PetscScalar)));
615: PetscCallHIP(hipMemcpy((void*)*d_q,q,ldq*mq*sizeof(PetscScalar),hipMemcpyHostToDevice));
616: PetscCall(PetscLogCpuToGpu(ldq*mq*sizeof(PetscScalar)));
617: }
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /*
622: BV_MatDenseHIPRestoreArrayRead - restores the pointer obtained with BV_MatDenseHIPGetArrayRead(),
623: freeing the GPU memory in case of MATSEQDENSE
624: */
625: static inline PetscErrorCode BV_MatDenseHIPRestoreArrayRead(PETSC_UNUSED BV bv,Mat Q,const PetscScalar **d_q)
626: {
627: PetscBool matiship;
629: PetscFunctionBegin;
630: PetscCall(PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSEHIP,&matiship));
631: if (matiship) PetscCall(MatDenseHIPRestoreArrayRead(Q,d_q));
632: else {
633: PetscCall(MatDenseRestoreArrayRead(Q,NULL));
634: PetscCallHIP(hipFree((void*)*d_q));
635: *d_q = NULL;
636: }
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: /*
641: BV_VecPlaceArray - allows the user to replace the device array in a vector with a
642: device array provided by the caller, calling the accessor that corresponds to the
643: vector type being used by the BV object
644: */
645: static inline PetscErrorCode BV_VecPlaceArray(PETSC_UNUSED BV bv,Vec v,PetscScalar *a)
646: {
647: PetscFunctionBegin;
648: #if PetscDefined(HAVE_KOKKOS_KERNELS)
649: if (bv->kokkos) {
650: PetscCall(VecKokkosPlaceArray(v,a));
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
653: #endif
654: PetscCall(VecHIPPlaceArray(v,a));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*
659: BV_VecResetArray - resets the device array in a vector to the one it had before the
660: call to BV_VecPlaceArray()
661: */
662: static inline PetscErrorCode BV_VecResetArray(PETSC_UNUSED BV bv,Vec v)
663: {
664: PetscFunctionBegin;
665: #if PetscDefined(HAVE_KOKKOS_KERNELS)
666: if (bv->kokkos) {
667: PetscCall(VecKokkosResetArray(v));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
670: #endif
671: PetscCall(VecHIPResetArray(v));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: SLEPC_INTERN PetscErrorCode BVMult_BLAS_HIP(BV,PetscInt,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscScalar,PetscScalar*,PetscInt);
676: SLEPC_INTERN PetscErrorCode BVMultVec_BLAS_HIP(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,const PetscScalar*,PetscScalar,PetscScalar*);
677: SLEPC_INTERN PetscErrorCode BVMultInPlace_BLAS_HIP(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscBool);
678: SLEPC_INTERN PetscErrorCode BVAXPY_BLAS_HIP(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscInt,PetscScalar,PetscScalar*,PetscInt);
679: SLEPC_INTERN PetscErrorCode BVDot_BLAS_HIP(BV,PetscInt,PetscInt,PetscInt,const PetscScalar*,PetscInt,const PetscScalar*,PetscInt,PetscScalar*,PetscInt,PetscBool);
680: SLEPC_INTERN PetscErrorCode BVDotVec_BLAS_HIP(BV,PetscInt,PetscInt,const PetscScalar*,PetscInt,const PetscScalar*,PetscScalar*,PetscBool);
681: SLEPC_INTERN PetscErrorCode BVScale_BLAS_HIP(BV,PetscInt,PetscScalar*,PetscScalar);
682: SLEPC_INTERN PetscErrorCode BVNorm_BLAS_HIP(BV,PetscInt,const PetscScalar*,PetscReal*);
683: SLEPC_INTERN PetscErrorCode BVNormalize_BLAS_HIP(BV,PetscInt,PetscInt,PetscScalar*,PetscInt,PetscScalar*);
685: SLEPC_INTERN PetscErrorCode BV_CleanCoefficients_HIP(BV,PetscInt,PetscScalar*);
686: SLEPC_INTERN PetscErrorCode BV_AddCoefficients_HIP(BV,PetscInt,PetscScalar*,PetscScalar*);
687: SLEPC_INTERN PetscErrorCode BV_SetValue_HIP(BV,PetscInt,PetscInt,PetscScalar*,PetscScalar);
688: SLEPC_INTERN PetscErrorCode BV_SquareSum_HIP(BV,PetscInt,PetscScalar*,PetscReal*);
689: SLEPC_INTERN PetscErrorCode BV_ApplySignature_HIP(BV,PetscInt,PetscScalar*,PetscBool);
690: SLEPC_INTERN PetscErrorCode BV_SquareRoot_HIP(BV,PetscInt,PetscScalar*,PetscReal*);
691: SLEPC_INTERN PetscErrorCode BV_StoreCoefficients_HIP(BV,PetscInt,PetscScalar*,PetscScalar*);
692: #define BV_CleanCoefficients(a,b,c) ((a)->hip?BV_CleanCoefficients_HIP:BV_CleanCoefficients_Default)((a),(b),(c))
693: #define BV_AddCoefficients(a,b,c,d) ((a)->hip?BV_AddCoefficients_HIP:BV_AddCoefficients_Default)((a),(b),(c),(d))
694: #define BV_SetValue(a,b,c,d,e) ((a)->hip?BV_SetValue_HIP:BV_SetValue_Default)((a),(b),(c),(d),(e))
695: #define BV_SquareSum(a,b,c,d) ((a)->hip?BV_SquareSum_HIP:BV_SquareSum_Default)((a),(b),(c),(d))
696: #define BV_ApplySignature(a,b,c,d) ((a)->hip?BV_ApplySignature_HIP:BV_ApplySignature_Default)((a),(b),(c),(d))
697: #define BV_SquareRoot(a,b,c,d) ((a)->hip?BV_SquareRoot_HIP:BV_SquareRoot_Default)((a),(b),(c),(d))
698: #define BV_StoreCoefficients(a,b,c,d) ((a)->hip?BV_StoreCoefficients_HIP:BV_StoreCoefficients_Default)((a),(b),(c),(d))
700: #else /* CPU */
701: #define BV_CleanCoefficients(a,b,c) BV_CleanCoefficients_Default((a),(b),(c))
702: #define BV_AddCoefficients(a,b,c,d) BV_AddCoefficients_Default((a),(b),(c),(d))
703: #define BV_SetValue(a,b,c,d,e) BV_SetValue_Default((a),(b),(c),(d),(e))
704: #define BV_SquareSum(a,b,c,d) BV_SquareSum_Default((a),(b),(c),(d))
705: #define BV_ApplySignature(a,b,c,d) BV_ApplySignature_Default((a),(b),(c),(d))
706: #define BV_SquareRoot(a,b,c,d) BV_SquareRoot_Default((a),(b),(c),(d))
707: #define BV_StoreCoefficients(a,b,c,d) BV_StoreCoefficients_Default((a),(b),(c),(d))
708: #endif /* PETSC_HAVE_CUDA */