Actual source code: bvimpl.h

slepc-3.6.3 2016-03-29
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) , Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.

  8:    SLEPc is free software: you can redistribute it and/or modify it under  the
  9:    terms of version 3 of the GNU Lesser General Public License as published by
 10:    the Free Software Foundation.

 12:    SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
 13:    WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
 14:    FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
 15:    more details.

 17:    You  should have received a copy of the GNU Lesser General  Public  License
 18:    along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
 19:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 20: */

 22: #if !defined(_BVIMPL)
 23: #define _BVIMPL

 25: #include <slepcbv.h>
 26: #include <slepc/private/slepcimpl.h>

 28: PETSC_EXTERN PetscBool BVRegisterAllCalled;
 29: PETSC_EXTERN PetscErrorCode BVRegisterAll(void);

 31: PETSC_EXTERN PetscLogEvent BV_Create,BV_Copy,BV_Mult,BV_Dot,BV_Orthogonalize,BV_Scale,BV_Norm,BV_SetRandom,BV_MatMult,BV_MatProject,BV_AXPY;

 33: typedef struct _BVOps *BVOps;

 35: struct _BVOps {
 36:   PetscErrorCode (*mult)(BV,PetscScalar,PetscScalar,BV,Mat);
 37:   PetscErrorCode (*multvec)(BV,PetscScalar,PetscScalar,Vec,PetscScalar*);
 38:   PetscErrorCode (*multinplace)(BV,Mat,PetscInt,PetscInt);
 39:   PetscErrorCode (*multinplacetrans)(BV,Mat,PetscInt,PetscInt);
 40:   PetscErrorCode (*axpy)(BV,PetscScalar,BV);
 41:   PetscErrorCode (*dot)(BV,BV,Mat);
 42:   PetscErrorCode (*dotvec)(BV,Vec,PetscScalar*);
 43:   PetscErrorCode (*dotvec_local)(BV,Vec,PetscScalar*);
 44:   PetscErrorCode (*dotvec_begin)(BV,Vec,PetscScalar*);
 45:   PetscErrorCode (*dotvec_end)(BV,Vec,PetscScalar*);
 46:   PetscErrorCode (*scale)(BV,PetscInt,PetscScalar);
 47:   PetscErrorCode (*norm)(BV,PetscInt,NormType,PetscReal*);
 48:   PetscErrorCode (*norm_local)(BV,PetscInt,NormType,PetscReal*);
 49:   PetscErrorCode (*norm_begin)(BV,PetscInt,NormType,PetscReal*);
 50:   PetscErrorCode (*norm_end)(BV,PetscInt,NormType,PetscReal*);
 51:   PetscErrorCode (*orthogonalize)(BV,Mat);
 52:   PetscErrorCode (*matmult)(BV,Mat,BV);
 53:   PetscErrorCode (*copy)(BV,BV);
 54:   PetscErrorCode (*resize)(BV,PetscInt,PetscBool);
 55:   PetscErrorCode (*getcolumn)(BV,PetscInt,Vec*);
 56:   PetscErrorCode (*restorecolumn)(BV,PetscInt,Vec*);
 57:   PetscErrorCode (*getarray)(BV,PetscScalar**);
 58:   PetscErrorCode (*restorearray)(BV,PetscScalar**);
 59:   PetscErrorCode (*create)(BV);
 60:   PetscErrorCode (*setfromoptions)(PetscOptions*,BV);
 61:   PetscErrorCode (*view)(BV,PetscViewer);
 62:   PetscErrorCode (*destroy)(BV);
 63: };

 65: struct _p_BV {
 66:   PETSCHEADER(struct _BVOps);
 67:   /*------------------------- User parameters --------------------------*/
 68:   Vec                t;            /* template vector */
 69:   PetscInt           n,N;          /* dimensions of vectors (local, global) */
 70:   PetscInt           m;            /* number of vectors */
 71:   PetscInt           l;            /* number of leading columns */
 72:   PetscInt           k;            /* number of active columns */
 73:   PetscInt           nc;           /* number of constraints */
 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 */

 82:   /*---------------------- Cached data and workspace -------------------*/
 83:   Vec                Bx;           /* result of matrix times a vector x */
 84:   PetscObjectId      xid;          /* object id of vector x */
 85:   PetscObjectState   xstate;       /* state of vector x */
 86:   Vec                cv[2];        /* column vectors obtained with BVGetColumn() */
 87:   PetscInt           ci[2];        /* column indices of obtained vectors */
 88:   PetscObjectState   st[2];        /* state of obtained vectors */
 89:   PetscObjectId      id[2];        /* object id of obtained vectors */
 90:   PetscScalar        *h,*c;        /* orthogonalization coefficients */
 91:   PetscReal          *omega;       /* signature matrix values for indefinite case */
 92:   Mat                B,C;          /* auxiliary dense matrices for matmult operation */
 93:   PetscObjectId      Aid;          /* object id of matrix A of matmult operation */
 94:   PetscBool          defersfo;     /* deferred call to setfromoptions */
 95:   BV                 cached;       /* cached BV to store result of matrix times BV */
 96:   PetscObjectState   bvstate;      /* state of BV when BVApplyMatrixBV() was called */
 97:   PetscScalar        *work;
 98:   PetscInt           lwork;
 99:   void               *data;
100: };

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

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

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

144:   if (((PetscObject)x)->id != bv->xid || ((PetscObject)x)->state != bv->xstate) {
145:     MatMult(bv->matrix,x,bv->Bx);
146:     bv->xid = ((PetscObject)x)->id;
147:     bv->xstate = ((PetscObject)x)->state;
148:   }
149:   return(0);
150: }

154: /*
155:   BV_AllocateCachedBV - Allocate auxiliary BV required for BVApplyMatrixBV if not available.
156: */
157: PETSC_STATIC_INLINE PetscErrorCode BV_AllocateCachedBV(BV V)
158: {

162:   if (!V->cached) {
163:     BVCreate(PetscObjectComm((PetscObject)V),&V->cached);
164:     BVSetSizesFromVec(V->cached,V->t,V->m);
165:     BVSetType(V->cached,((PetscObject)V)->type_name);
166:     BVSetOrthogonalization(V->cached,V->orthog_type,V->orthog_ref,V->orthog_eta,V->orthog_block);
167:   }
168:   return(0);
169: }

173: /*
174:   BV_IPMatMultBV - Multiply BV by the inner-product matrix, cache the
175:   result internally in bv->cached.
176: */
177: PETSC_STATIC_INLINE PetscErrorCode BV_IPMatMultBV(BV bv)
178: {

182:   BV_AllocateCachedBV(bv);
183:   BVSetActiveColumns(bv->cached,bv->l,bv->k);
184:   if (((PetscObject)bv)->state != bv->bvstate) {
185:     if (bv->matrix) {
186:       BVMatMult(bv,bv->matrix,bv->cached);
187:     } else {
188:       BVCopy(bv,bv->cached);
189:     }
190:     bv->bvstate = ((PetscObject)bv)->state;
191:   }
192:   return(0);
193: }

197: /*
198:   BV_AllocateCoeffs - Allocate orthogonalization coefficients if not done already.
199: */
200: PETSC_STATIC_INLINE PetscErrorCode BV_AllocateCoeffs(BV bv)
201: {

205:   if (!bv->h) {
206:     PetscMalloc2(bv->nc+bv->m,&bv->h,bv->nc+bv->m,&bv->c);
207:     PetscLogObjectMemory((PetscObject)bv,2*bv->m*sizeof(PetscScalar));
208:   }
209:   return(0);
210: }

214: /*
215:   BV_AllocateSignature - Allocate signature coefficients if not done already.
216: */
217: PETSC_STATIC_INLINE PetscErrorCode BV_AllocateSignature(BV bv)
218: {
220:   PetscInt       i;

223:   if (bv->indef && !bv->omega) {
224:     PetscMalloc1(bv->nc+bv->m,&bv->omega);
225:     PetscLogObjectMemory((PetscObject)bv,bv->m*sizeof(PetscReal));
226:     for (i=-bv->nc;i<bv->m;i++) bv->omega[i] = 1.0;
227:   }
228:   return(0);
229: }

233: /*
234:   BV_AllocateMatMult - Allocate auxiliary matrices required for BVMatMult if not available.
235: */
236: PETSC_STATIC_INLINE PetscErrorCode BV_AllocateMatMult(BV bv,Mat A,PetscInt m)
237: {
239:   PetscObjectId  Aid;
240:   PetscBool      create=PETSC_FALSE;
241:   PetscInt       cols;

244:   if (!bv->B) create=PETSC_TRUE;
245:   else {
246:     MatGetSize(bv->B,NULL,&cols);
247:     PetscObjectGetId((PetscObject)A,&Aid);
248:     if (cols!=m || bv->Aid!=Aid) {
249:       MatDestroy(&bv->B);
250:       MatDestroy(&bv->C);
251:       create=PETSC_TRUE;
252:     }
253:   }
254:   if (create) {
255:     MatCreateDense(PetscObjectComm((PetscObject)bv),bv->n,PETSC_DECIDE,bv->N,m,NULL,&bv->B);
256:     PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->B);
257:     MatAssemblyBegin(bv->B,MAT_FINAL_ASSEMBLY);
258:     MatAssemblyEnd(bv->B,MAT_FINAL_ASSEMBLY);
259:   }
260:   return(0);
261: }

263: /*
264:   BVAvailableVec: First (0) or second (1) vector available for
265:   getcolumn operation (or -1 if both vectors already fetched).
266: */
267: #define BVAvailableVec (((bv->ci[0]==-bv->nc-1)? 0: (bv->ci[1]==-bv->nc-1)? 1: -1))

269: /*
270:     Macros to test valid BV arguments
271: */
272: #if !defined(PETSC_USE_DEBUG)

274: #define BVCheckSizes(h,arg) do {} while (0)

276: #else

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

283: #endif

285: PETSC_INTERN PetscErrorCode BVView_Vecs(BV,PetscViewer);

287: PETSC_INTERN PetscErrorCode BVAllocateWork_Private(BV,PetscInt);

289: PETSC_INTERN PetscErrorCode BVMult_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar,const PetscScalar*,const PetscScalar*,PetscScalar,PetscScalar*);
290: PETSC_INTERN PetscErrorCode BVMultVec_BLAS_Private(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,const PetscScalar*,PetscScalar,PetscScalar*);
291: PETSC_INTERN PetscErrorCode BVMultInPlace_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscScalar*,const PetscScalar*,PetscBool);
292: PETSC_INTERN PetscErrorCode BVMultInPlace_Vecs_Private(BV,PetscInt,PetscInt,PetscInt,Vec*,const PetscScalar*,PetscBool);
293: PETSC_INTERN PetscErrorCode BVAXPY_BLAS_Private(BV,PetscInt,PetscInt,PetscScalar,const PetscScalar*,PetscScalar*);
294: PETSC_INTERN PetscErrorCode BVDot_BLAS_Private(BV,PetscInt,PetscInt,PetscInt,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscBool);
295: PETSC_INTERN PetscErrorCode BVDotVec_BLAS_Private(BV,PetscInt,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscBool);
296: PETSC_INTERN PetscErrorCode BVScale_BLAS_Private(BV,PetscInt,PetscScalar*,PetscScalar);
297: PETSC_INTERN PetscErrorCode BVNorm_LAPACK_Private(BV,PetscInt,PetscInt,const PetscScalar*,NormType,PetscReal*,PetscBool);
298: PETSC_INTERN PetscErrorCode BVOrthogonalize_LAPACK_Private(BV,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscBool);

300: #endif