Actual source code: bvimpl.h

slepc-3.12.2 2020-01-13
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: #if !defined(SLEPCBVIMPL_H)
 12: #define SLEPCBVIMPL_H

 14: #include <slepcbv.h>
 15: #include <slepc/private/slepcimpl.h>

 17: SLEPC_EXTERN PetscBool BVRegisterAllCalled;
 18: SLEPC_EXTERN PetscErrorCode BVRegisterAll(void);

 20: 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_SetRandom,BV_MatMult,BV_MatMultVec,BV_MatProject;

 22: typedef struct _BVOps *BVOps;

 24: struct _BVOps {
 25:   PetscErrorCode (*mult)(BV,PetscScalar,PetscScalar,BV,Mat);
 26:   PetscErrorCode (*multvec)(BV,PetscScalar,PetscScalar,Vec,PetscScalar*);
 27:   PetscErrorCode (*multinplace)(BV,Mat,PetscInt,PetscInt);
 28:   PetscErrorCode (*multinplacetrans)(BV,Mat,PetscInt,PetscInt);
 29:   PetscErrorCode (*dot)(BV,BV,Mat);
 30:   PetscErrorCode (*dotvec)(BV,Vec,PetscScalar*);
 31:   PetscErrorCode (*dotvec_local)(BV,Vec,PetscScalar*);
 32:   PetscErrorCode (*dotvec_begin)(BV,Vec,PetscScalar*);
 33:   PetscErrorCode (*dotvec_end)(BV,Vec,PetscScalar*);
 34:   PetscErrorCode (*scale)(BV,PetscInt,PetscScalar);
 35:   PetscErrorCode (*norm)(BV,PetscInt,NormType,PetscReal*);
 36:   PetscErrorCode (*norm_local)(BV,PetscInt,NormType,PetscReal*);
 37:   PetscErrorCode (*norm_begin)(BV,PetscInt,NormType,PetscReal*);
 38:   PetscErrorCode (*norm_end)(BV,PetscInt,NormType,PetscReal*);
 39:   PetscErrorCode (*matmult)(BV,Mat,BV);
 40:   PetscErrorCode (*copy)(BV,BV);
 41:   PetscErrorCode (*copycolumn)(BV,PetscInt,PetscInt);
 42:   PetscErrorCode (*resize)(BV,PetscInt,PetscBool);
 43:   PetscErrorCode (*getcolumn)(BV,PetscInt,Vec*);
 44:   PetscErrorCode (*restorecolumn)(BV,PetscInt,Vec*);
 45:   PetscErrorCode (*getarray)(BV,PetscScalar**);
 46:   PetscErrorCode (*restorearray)(BV,PetscScalar**);
 47:   PetscErrorCode (*getarrayread)(BV,const PetscScalar**);
 48:   PetscErrorCode (*restorearrayread)(BV,const PetscScalar**);
 49:   PetscErrorCode (*restoresplit)(BV,BV*,BV*);
 50:   PetscErrorCode (*gramschmidt)(BV,PetscInt,Vec,PetscBool*,PetscScalar*,PetscScalar*,PetscReal*,PetscReal*);
 51:   PetscErrorCode (*duplicate)(BV,BV);
 52:   PetscErrorCode (*create)(BV);
 53:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,BV);
 54:   PetscErrorCode (*view)(BV,PetscViewer);
 55:   PetscErrorCode (*destroy)(BV);
 56: };

 58: struct _p_BV {
 59:   PETSCHEADER(struct _BVOps);
 60:   /*------------------------- User parameters --------------------------*/
 61:   Vec                t;            /* template vector */
 62:   PetscInt           n,N;          /* dimensions of vectors (local, global) */
 63:   PetscInt           m;            /* number of vectors */
 64:   PetscInt           l;            /* number of leading columns */
 65:   PetscInt           k;            /* number of active columns */
 66:   PetscInt           nc;           /* number of constraints */
 67:   BVOrthogType       orthog_type;  /* the method of vector orthogonalization */
 68:   BVOrthogRefineType orthog_ref;   /* refinement method */
 69:   PetscReal          orthog_eta;   /* refinement threshold */
 70:   BVOrthogBlockType  orthog_block; /* the method of block orthogonalization */
 71:   Mat                matrix;       /* inner product matrix */
 72:   PetscBool          indef;        /* matrix is indefinite */
 73:   BVMatMultType      vmm;          /* version of matmult operation */
 74:   PetscBool          rrandom;      /* reproducible random vectors */

 76:   /*---------------------- Cached data and workspace -------------------*/
 77:   Vec                buffer;       /* buffer vector used in orthogonalization */
 78:   Mat                Abuffer;      /* auxiliary seqdense matrix that wraps the buffer */
 79:   Vec                Bx;           /* result of matrix times a vector x */
 80:   PetscObjectId      xid;          /* object id of vector x */
 81:   PetscObjectState   xstate;       /* state of vector x */
 82:   Vec                cv[2];        /* column vectors obtained with BVGetColumn() */
 83:   PetscInt           ci[2];        /* column indices of obtained vectors */
 84:   PetscObjectState   st[2];        /* state of obtained vectors */
 85:   PetscObjectId      id[2];        /* object id of obtained vectors */
 86:   PetscScalar        *h,*c;        /* orthogonalization coefficients */
 87:   Vec                omega;        /* signature matrix values for indefinite case */
 88:   PetscBool          defersfo;     /* deferred call to setfromoptions */
 89:   BV                 cached;       /* cached BV to store result of matrix times BV */
 90:   PetscObjectState   bvstate;      /* state of BV when BVApplyMatrixBV() was called */
 91:   BV                 L,R;          /* BV objects obtained with BVGetSplit() */
 92:   PetscObjectState   lstate,rstate;/* state of L and R when BVGetSplit() was called */
 93:   PetscInt           lsplit;       /* the value of l when BVGetSplit() was called */
 94:   PetscInt           issplit;      /* >0 if this BV has been created by splitting (1=left, 2=right) */
 95:   BV                 splitparent;  /* my parent if I am a split BV */
 96:   PetscRandom        rand;         /* random number generator */
 97:   Mat                Acreate;      /* matrix given at BVCreateFromMat() */
 98:   Mat                Aget;         /* matrix returned for BVGetMat() */
 99:   Mat                Amult;        /* matrix argument of last call to BVMatMult() */
100:   PetscObjectState   Amultstate;   /* state of Amult */
101:   PetscBool          cuda;         /* true if GPU must be used in SVEC */
102:   PetscScalar        *work;
103:   PetscInt           lwork;
104:   void               *data;
105: };

107: /*
108:   BV_SafeSqrt - Computes the square root of a scalar value alpha, which is
109:   assumed to be z'*B*z. The result is
110:     if definite inner product:     res = sqrt(alpha)
111:     if indefinite inner product:   res = sgn(alpha)*sqrt(abs(alpha))
112: */
113: PETSC_STATIC_INLINE PetscErrorCode BV_SafeSqrt(BV bv,PetscScalar alpha,PetscReal *res)
114: {
116:   PetscReal      absal,realp;

119:   absal = PetscAbsScalar(alpha);
120:   realp = PetscRealPart(alpha);
121:   if (absal<PETSC_MACHINE_EPSILON) {
122:     PetscInfo(bv,"Zero norm, either the vector is zero or a semi-inner product is being used\n");
123:   }
124: #if defined(PETSC_USE_COMPLEX)
125:   if (PetscAbsReal(PetscImaginaryPart(alpha))>10*PETSC_MACHINE_EPSILON && PetscAbsReal(PetscImaginaryPart(alpha))/absal>100*PETSC_MACHINE_EPSILON) SETERRQ1(PetscObjectComm((PetscObject)bv),1,"The inner product is not well defined: nonzero imaginary part %g",PetscImaginaryPart(alpha));
126: #endif
127:   if (bv->indef) {
128:     *res = (realp<0.0)? -PetscSqrtReal(-realp): PetscSqrtReal(realp);
129:   } else {
130:     if (realp<-10*PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)bv),1,"The inner product is not well defined: indefinite matrix");
131:     *res = (realp<0.0)? 0.0: PetscSqrtReal(realp);
132:   }
133:   return(0);
134: }

136: /*
137:   BV_IPMatMult - Multiply a vector x by the inner-product matrix, cache the
138:   result in Bx.
139: */
140: PETSC_STATIC_INLINE PetscErrorCode BV_IPMatMult(BV bv,Vec x)
141: {

145:   if (((PetscObject)x)->id != bv->xid || ((PetscObject)x)->state != bv->xstate) {
146:     if (!bv->Bx) {
147:       MatCreateVecs(bv->matrix,&bv->Bx,NULL);
148:       PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Bx);
149:     }
150:     MatMult(bv->matrix,x,bv->Bx);
151:     PetscObjectGetId((PetscObject)x,&bv->xid);
152:     PetscObjectStateGet((PetscObject)x,&bv->xstate);
153:   }
154:   return(0);
155: }

157: /*
158:   BV_IPMatMultBV - Multiply BV by the inner-product matrix, cache the
159:   result internally in bv->cached.
160: */
161: PETSC_STATIC_INLINE PetscErrorCode BV_IPMatMultBV(BV bv)
162: {

166:   BVGetCachedBV(bv,&bv->cached);
167:   if (((PetscObject)bv)->state != bv->bvstate || bv->l != bv->cached->l || bv->k != bv->cached->k) {
168:     BVSetActiveColumns(bv->cached,bv->l,bv->k);
169:     if (bv->matrix) {
170:       BVMatMult(bv,bv->matrix,bv->cached);
171:     } else {
172:       BVCopy(bv,bv->cached);
173:     }
174:     bv->bvstate = ((PetscObject)bv)->state;
175:   }
176:   return(0);
177: }

179: /*
180:   BV_AllocateCoeffs - Allocate orthogonalization coefficients if not done already.
181: */
182: PETSC_STATIC_INLINE PetscErrorCode BV_AllocateCoeffs(BV bv)
183: {

187:   if (!bv->h) {
188:     PetscMalloc2(bv->nc+bv->m,&bv->h,bv->nc+bv->m,&bv->c);
189:     PetscLogObjectMemory((PetscObject)bv,2*bv->m*sizeof(PetscScalar));
190:   }
191:   return(0);
192: }

194: /*
195:   BV_AllocateSignature - Allocate signature coefficients if not done already.
196: */
197: PETSC_STATIC_INLINE PetscErrorCode BV_AllocateSignature(BV bv)
198: {

202:   if (bv->indef && !bv->omega) {
203:     if (bv->cuda) {
204: #if defined(PETSC_HAVE_CUDA)
205:       VecCreateSeqCUDA(PETSC_COMM_SELF,bv->nc+bv->m,&bv->omega);
206: #else
207:       SETERRQ(PetscObjectComm((PetscObject)bv),1,"Something wrong happened");
208: #endif
209:     } else {
210:       VecCreateSeq(PETSC_COMM_SELF,bv->nc+bv->m,&bv->omega);
211:     }
212:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->omega);
213:     VecSet(bv->omega,1.0);
214:   }
215:   return(0);
216: }

218: /*
219:   BVAvailableVec: First (0) or second (1) vector available for
220:   getcolumn operation (or -1 if both vectors already fetched).
221: */
222: #define BVAvailableVec (((bv->ci[0]==-bv->nc-1)? 0: (bv->ci[1]==-bv->nc-1)? 1: -1))

224: /*
225:     Macros to test valid BV arguments
226: */
227: #if !defined(PETSC_USE_DEBUG)

229: #define BVCheckSizes(h,arg) do {} while (0)
230: #define BVCheckOp(h,arg,op) do {} while (0)

232: #else

234: #define BVCheckSizes(h,arg) \
235:   do { \
236:     if (!h->m) SETERRQ1(PetscObjectComm((PetscObject)h),PETSC_ERR_ARG_WRONGSTATE,"BV sizes have not been defined: Parameter #%d",arg); \
237:   } while (0)

239: #define BVCheckOp(h,arg,op) \
240:   do { \
241:     if (!h->ops->op) SETERRQ1(PetscObjectComm((PetscObject)h),PETSC_ERR_SUP,"Operation not implemented in this BV type: Parameter #%d",arg); \
242:   } while (0)

244: #endif

246: SLEPC_INTERN PetscErrorCode BVView_Vecs(BV,PetscViewer);

248: SLEPC_INTERN PetscErrorCode BVAllocateWork_Private(BV,PetscInt);

250: SLEPC_INTERN PetscErrorCode BVMult_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar,const PetscScalar*,const PetscScalar*,PetscScalar,PetscScalar*);
251: SLEPC_INTERN PetscErrorCode BVMultVec_BLAS_Private(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,const PetscScalar*,PetscScalar,PetscScalar*);
252: SLEPC_INTERN PetscErrorCode BVMultInPlace_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar*,const PetscScalar*,PetscBool);
253: SLEPC_INTERN PetscErrorCode BVMultInPlace_Vecs_Private(BV,PetscInt,PetscInt,PetscInt,Vec*,const PetscScalar*,PetscBool);
254: SLEPC_INTERN PetscErrorCode BVAXPY_BLAS_Private(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscScalar,PetscScalar*);
255: SLEPC_INTERN PetscErrorCode BVDot_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscBool);
256: SLEPC_INTERN PetscErrorCode BVDotVec_BLAS_Private(BV,PetscInt,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscBool);
257: SLEPC_INTERN PetscErrorCode BVScale_BLAS_Private(BV,PetscInt,PetscScalar*,PetscScalar);
258: SLEPC_INTERN PetscErrorCode BVNorm_LAPACK_Private(BV,PetscInt,PetscInt,const PetscScalar*,NormType,PetscReal*,PetscBool);
259: SLEPC_INTERN PetscErrorCode BVMatCholInv_LAPACK_Private(BV,Mat,Mat);
260: SLEPC_INTERN PetscErrorCode BVMatTriInv_LAPACK_Private(BV,Mat,Mat);
261: SLEPC_INTERN PetscErrorCode BVMatSVQB_LAPACK_Private(BV,Mat,Mat);
262: SLEPC_INTERN PetscErrorCode BVOrthogonalize_LAPACK_TSQR(BV,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscInt);
263: SLEPC_INTERN PetscErrorCode BVOrthogonalize_LAPACK_TSQR_OnlyR(BV,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscInt);

265: /* reduction operation used in BVOrthogonalize */
266: SLEPC_EXTERN MPI_Op MPIU_TSQR;
267: SLEPC_EXTERN void MPIAPI SlepcGivensPacked(void*,void*,PetscMPIInt*,MPI_Datatype*);

269: #if defined(PETSC_HAVE_CUDA)

271: #define WaitForGPU() PetscCUDASynchronize ? cudaDeviceSynchronize() : 0

273: /* complex single */
274: #if defined(PETSC_USE_COMPLEX)
275: #if defined(PETSC_USE_REAL_SINGLE)
276: #define cublasXgemm(a,b,c,d,e,f,g,h,i,j,k,l,m,n) cublasCgemm((a),(b),(c),(d),(e),(f),(cuComplex*)(g),(cuComplex*)(h),(i),(cuComplex*)(j),(k),(cuComplex*)(l),(cuComplex*)(m),(n))
277: #define cublasXgemv(a,b,c,d,e,f,g,h,i,j,k,l) cublasCgemv((a),(b),(c),(d),(cuComplex*)(e),(cuComplex*)(f),(g),(cuComplex*)(h),(i),(cuComplex*)(j),(cuComplex*)(k),(l))
278: #define cublasXscal(a,b,c,d,e) cublasCscal(a,b,(const cuComplex*)(c),(cuComplex*)(d),e)
279: #define cublasXnrm2(a,b,c,d,e) cublasScnrm2(a,b,(const cuComplex*)(c),d,e)
280: #define cublasXaxpy(a,b,c,d,e,f,g) cublasCaxpy((a),(b),(cuComplex*)(c),(cuComplex*)(d),(e),(cuComplex*)(f),(g))
281: #define cublasXdotc(a,b,c,d,e,f,g) cublasCdotc((a),(b),(const cuComplex *)(c),(d),(const cuComplex *)(e),(f),(cuComplex *)(g))
282: #else /* complex double */
283: #define cublasXgemm(a,b,c,d,e,f,g,h,i,j,k,l,m,n) cublasZgemm((a),(b),(c),(d),(e),(f),(cuDoubleComplex*)(g),(cuDoubleComplex*)(h),(i),(cuDoubleComplex*)(j),(k),(cuDoubleComplex*)(l),(cuDoubleComplex*)(m),(n))
284: #define cublasXgemv(a,b,c,d,e,f,g,h,i,j,k,l) cublasZgemv((a),(b),(c),(d),(cuDoubleComplex*)(e),(cuDoubleComplex*)(f),(g),(cuDoubleComplex*)(h),(i),(cuDoubleComplex*)(j),(cuDoubleComplex*)(k),(l))
285: #define cublasXscal(a,b,c,d,e) cublasZscal(a,b,(const cuDoubleComplex*)(c),(cuDoubleComplex*)(d),e)
286: #define cublasXnrm2(a,b,c,d,e) cublasDznrm2(a,b,(const cuDoubleComplex*)(c),d,e)
287: #define cublasXaxpy(a,b,c,d,e,f,g) cublasZaxpy((a),(b),(cuDoubleComplex*)(c),(cuDoubleComplex*)(d),(e),(cuDoubleComplex*)(f),(g))
288: #define cublasXdotc(a,b,c,d,e,f,g) cublasZdotc((a),(b),(const cuDoubleComplex *)(c),(d),(const cuDoubleComplex *)(e),(f),(cuDoubleComplex *)(g))
289: #endif
290: #else /* real single */
291: #if defined(PETSC_USE_REAL_SINGLE)
292: #define cublasXgemm cublasSgemm
293: #define cublasXgemv cublasSgemv
294: #define cublasXscal cublasSscal
295: #define cublasXnrm2 cublasSnrm2
296: #define cublasXaxpy cublasSaxpy
297: #define cublasXdotc cublasSdot
298: #else /* real double */
299: #define cublasXgemm cublasDgemm
300: #define cublasXgemv cublasDgemv
301: #define cublasXscal cublasDscal
302: #define cublasXnrm2 cublasDnrm2
303: #define cublasXaxpy cublasDaxpy
304: #define cublasXdotc cublasDdot
305: #endif
306: #endif

308: SLEPC_INTERN PetscErrorCode BV_CleanCoefficients_CUDA(BV,PetscInt,PetscScalar*);
309: SLEPC_INTERN PetscErrorCode BV_AddCoefficients_CUDA(BV,PetscInt,PetscScalar*,PetscScalar*);
310: SLEPC_INTERN PetscErrorCode BV_SetValue_CUDA(BV,PetscInt,PetscInt,PetscScalar*,PetscScalar);
311: SLEPC_INTERN PetscErrorCode BV_SquareSum_CUDA(BV,PetscInt,PetscScalar*,PetscReal*);
312: SLEPC_INTERN PetscErrorCode BV_ApplySignature_CUDA(BV,PetscInt,PetscScalar*,PetscBool);
313: SLEPC_INTERN PetscErrorCode BV_SquareRoot_CUDA(BV,PetscInt,PetscScalar*,PetscReal*);
314: SLEPC_INTERN PetscErrorCode BV_StoreCoefficients_CUDA(BV,PetscInt,PetscScalar*,PetscScalar*);

316: #endif /* PETSC_HAVE_CUDA */

318: /*
319:    BV_CleanCoefficients_Default - Sets to zero all entries of column j of the bv buffer
320: */
321: PETSC_STATIC_INLINE PetscErrorCode BV_CleanCoefficients_Default(BV bv,PetscInt j,PetscScalar *h)
322: {
324:   PetscScalar    *hh=h,*a;
325:   PetscInt       i;

328:   if (!h) {
329:     VecGetArray(bv->buffer,&a);
330:     hh = a + j*(bv->nc+bv->m);
331:   }
332:   for (i=0;i<bv->nc+j;i++) hh[i] = 0.0;
333:   if (!h) { VecRestoreArray(bv->buffer,&a); }
334:   return(0);
335: }

337: /*
338:    BV_AddCoefficients_Default - Add the contents of the scratch (0-th column) of the bv buffer
339:    into column j of the bv buffer
340: */
341: PETSC_STATIC_INLINE PetscErrorCode BV_AddCoefficients_Default(BV bv,PetscInt j,PetscScalar *h,PetscScalar *c)
342: {
344:   PetscScalar    *hh=h,*cc=c;
345:   PetscInt       i;

348:   if (!h) {
349:     VecGetArray(bv->buffer,&cc);
350:     hh = cc + j*(bv->nc+bv->m);
351:   }
352:   for (i=0;i<bv->nc+j;i++) hh[i] += cc[i];
353:   if (!h) { VecRestoreArray(bv->buffer,&cc); }
354:   return(0);
355: }

357: /*
358:    BV_SetValue_Default - Sets value in row j (counted after the constraints) of column k
359:    of the coefficients array
360: */
361: PETSC_STATIC_INLINE PetscErrorCode BV_SetValue_Default(BV bv,PetscInt j,PetscInt k,PetscScalar *h,PetscScalar value)
362: {
364:   PetscScalar    *hh=h,*a;

367:   if (!h) {
368:     VecGetArray(bv->buffer,&a);
369:     hh = a + k*(bv->nc+bv->m);
370:   }
371:   hh[bv->nc+j] = value;
372:   if (!h) { VecRestoreArray(bv->buffer,&a); }
373:   return(0);
374: }

376: /*
377:    BV_SquareSum_Default - Returns the value h'*h, where h represents the contents of the
378:    coefficients array (up to position j)
379: */
380: PETSC_STATIC_INLINE PetscErrorCode BV_SquareSum_Default(BV bv,PetscInt j,PetscScalar *h,PetscReal *sum)
381: {
383:   PetscScalar    *hh=h;
384:   PetscInt       i;

387:   *sum = 0.0;
388:   if (!h) { VecGetArray(bv->buffer,&hh); }
389:   for (i=0;i<bv->nc+j;i++) *sum += PetscRealPart(hh[i]*PetscConj(hh[i]));
390:   if (!h) { VecRestoreArray(bv->buffer,&hh); }
391:   return(0);
392: }

394: /*
395:    BV_ApplySignature_Default - Computes the pointwise product h*omega, where h represents
396:    the contents of the coefficients array (up to position j) and omega is the signature;
397:    if inverse=TRUE then the operation is h/omega
398: */
399: PETSC_STATIC_INLINE PetscErrorCode BV_ApplySignature_Default(BV bv,PetscInt j,PetscScalar *h,PetscBool inverse)
400: {
401:   PetscErrorCode    ierr;
402:   PetscScalar       *hh=h;
403:   PetscInt          i;
404:   const PetscScalar *omega;

407:   if (!(bv->nc+j)) return(0);
408:   if (!h) { VecGetArray(bv->buffer,&hh); }
409:   VecGetArrayRead(bv->omega,&omega);
410:   if (inverse) for (i=0;i<bv->nc+j;i++) hh[i] /= PetscRealPart(omega[i]);
411:   else for (i=0;i<bv->nc+j;i++) hh[i] *= PetscRealPart(omega[i]);
412:   VecRestoreArrayRead(bv->omega,&omega);
413:   if (!h) { VecRestoreArray(bv->buffer,&hh); }
414:   return(0);
415: }

417: /*
418:    BV_SquareRoot_Default - Returns the square root of position j (counted after the constraints)
419:    of the coefficients array
420: */
421: PETSC_STATIC_INLINE PetscErrorCode BV_SquareRoot_Default(BV bv,PetscInt j,PetscScalar *h,PetscReal *beta)
422: {
424:   PetscScalar    *hh=h;

427:   if (!h) { VecGetArray(bv->buffer,&hh); }
428:   BV_SafeSqrt(bv,hh[bv->nc+j],beta);
429:   if (!h) { VecRestoreArray(bv->buffer,&hh); }
430:   return(0);
431: }

433: /*
434:    BV_StoreCoefficients_Default - Copy the contents of the coefficients array to an array dest
435:    provided by the caller (only values from l to j are copied)
436: */
437: PETSC_STATIC_INLINE PetscErrorCode BV_StoreCoefficients_Default(BV bv,PetscInt j,PetscScalar *h,PetscScalar *dest)
438: {
440:   PetscScalar    *hh=h,*a;
441:   PetscInt       i;

444:   if (!h) {
445:     VecGetArray(bv->buffer,&a);
446:     hh = a + j*(bv->nc+bv->m);
447:   }
448:   for (i=bv->l;i<j;i++) dest[i-bv->l] = hh[bv->nc+i];
449:   if (!h) { VecRestoreArray(bv->buffer,&a); }
450:   return(0);
451: }

453: /*
454:   BV_GetEigenvector - retrieves k-th eigenvector from basis vectors V.
455:   The argument eigi is the imaginary part of the corresponding eigenvalue.
456: */
457: PETSC_STATIC_INLINE PetscErrorCode BV_GetEigenvector(BV V,PetscInt k,PetscScalar eigi,Vec Vr,Vec Vi)
458: {

462: #if defined(PETSC_USE_COMPLEX)
463:   if (Vr) { BVCopyVec(V,k,Vr); }
464:   if (Vi) { VecSet(Vi,0.0); }
465: #else
466:   if (eigi > 0.0) { /* first value of conjugate pair */
467:     if (Vr) { BVCopyVec(V,k,Vr); }
468:     if (Vi) { BVCopyVec(V,k+1,Vi); }
469:   } else if (eigi < 0.0) { /* second value of conjugate pair */
470:     if (Vr) { BVCopyVec(V,k-1,Vr); }
471:     if (Vi) {
472:       BVCopyVec(V,k,Vi);
473:       VecScale(Vi,-1.0);
474:     }
475:   } else { /* real eigenvalue */
476:     if (Vr) { BVCopyVec(V,k,Vr); }
477:     if (Vi) { VecSet(Vi,0.0); }
478:   }
479: #endif
480:   return(0);
481: }

483: #if defined(PETSC_HAVE_CUDA)
484: #define BV_CleanCoefficients(a,b,c)   ((a)->cuda?BV_CleanCoefficients_CUDA:BV_CleanCoefficients_Default)((a),(b),(c))
485: #define BV_AddCoefficients(a,b,c,d)   ((a)->cuda?BV_AddCoefficients_CUDA:BV_AddCoefficients_Default)((a),(b),(c),(d))
486: #define BV_SetValue(a,b,c,d,e)        ((a)->cuda?BV_SetValue_CUDA:BV_SetValue_Default)((a),(b),(c),(d),(e))
487: #define BV_SquareSum(a,b,c,d)         ((a)->cuda?BV_SquareSum_CUDA:BV_SquareSum_Default)((a),(b),(c),(d))
488: #define BV_ApplySignature(a,b,c,d)    ((a)->cuda?BV_ApplySignature_CUDA:BV_ApplySignature_Default)((a),(b),(c),(d))
489: #define BV_SquareRoot(a,b,c,d)        ((a)->cuda?BV_SquareRoot_CUDA:BV_SquareRoot_Default)((a),(b),(c),(d))
490: #define BV_StoreCoefficients(a,b,c,d) ((a)->cuda?BV_StoreCoefficients_CUDA:BV_StoreCoefficients_Default)((a),(b),(c),(d))
491: #else
492: #define BV_CleanCoefficients(a,b,c)   BV_CleanCoefficients_Default((a),(b),(c))
493: #define BV_AddCoefficients(a,b,c,d)   BV_AddCoefficients_Default((a),(b),(c),(d))
494: #define BV_SetValue(a,b,c,d,e)        BV_SetValue_Default((a),(b),(c),(d),(e))
495: #define BV_SquareSum(a,b,c,d)         BV_SquareSum_Default((a),(b),(c),(d))
496: #define BV_ApplySignature(a,b,c,d)    BV_ApplySignature_Default((a),(b),(c),(d))
497: #define BV_SquareRoot(a,b,c,d)        BV_SquareRoot_Default((a),(b),(c),(d))
498: #define BV_StoreCoefficients(a,b,c,d) BV_StoreCoefficients_Default((a),(b),(c),(d))
499: #endif /* PETSC_HAVE_CUDA */

501: #endif