Actual source code: bvimpl.h
slepc-3.6.3 2016-03-29
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