Actual source code: epsimpl.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 <slepceps.h>
 14: #include <slepc/private/bvimpl.h>

 16: /* SUBMANSEC = EPS */

 18: SLEPC_EXTERN PetscBool EPSRegisterAllCalled;
 19: SLEPC_EXTERN PetscBool EPSMonitorRegisterAllCalled;
 20: SLEPC_EXTERN PetscErrorCode EPSRegisterAll(void);
 21: SLEPC_EXTERN PetscErrorCode EPSMonitorRegisterAll(void);
 22: SLEPC_EXTERN PetscLogEvent EPS_SetUp,EPS_Solve,EPS_CISS_SVD;

 24: typedef struct _EPSOps *EPSOps;

 26: struct _EPSOps {
 27:   PetscErrorCode (*solve)(EPS);
 28:   PetscErrorCode (*setup)(EPS);
 29:   PetscErrorCode (*setupsort)(EPS);
 30:   PetscErrorCode (*setfromoptions)(EPS,PetscOptionItems);
 31:   PetscErrorCode (*publishoptions)(EPS);
 32:   PetscErrorCode (*destroy)(EPS);
 33:   PetscErrorCode (*reset)(EPS);
 34:   PetscErrorCode (*view)(EPS,PetscViewer);
 35:   PetscErrorCode (*backtransform)(EPS);
 36:   PetscErrorCode (*computevectors)(EPS);
 37:   PetscErrorCode (*setdefaultst)(EPS);
 38:   PetscErrorCode (*setdstype)(EPS);
 39: };

 41: /*
 42:    Maximum number of monitors you can run with a single EPS
 43: */
 44: #define MAXEPSMONITORS 5

 46: /*
 47:    The solution process goes through several states
 48: */
 49: typedef enum { EPS_STATE_INITIAL,
 50:                EPS_STATE_SETUP,
 51:                EPS_STATE_SOLVED,
 52:                EPS_STATE_EIGENVECTORS } EPSStateType;

 54: /*
 55:    To classify the different solvers into categories
 56: */
 57: typedef enum { EPS_CATEGORY_KRYLOV,      /* Krylov solver: relies on STApply and STBackTransform (same as OTHER) */
 58:                EPS_CATEGORY_PRECOND,     /* Preconditioned solver: uses ST only to manage preconditioner */
 59:                EPS_CATEGORY_CONTOUR,     /* Contour integral: ST used to solve linear systems at integration points */
 60:                EPS_CATEGORY_OTHER } EPSSolverType;

 62: /*
 63:    To check for unsupported features at EPSSetUp_XXX()
 64: */
 65: typedef enum { EPS_FEATURE_BALANCE=1,       /* balancing */
 66:                EPS_FEATURE_ARBITRARY=2,     /* arbitrary selection of eigepairs */
 67:                EPS_FEATURE_REGION=4,        /* nontrivial region for filtering */
 68:                EPS_FEATURE_EXTRACTION=8,    /* extraction technique different from Ritz */
 69:                EPS_FEATURE_CONVERGENCE=16,  /* convergence test selected by user */
 70:                EPS_FEATURE_STOPPING=32,     /* stopping test */
 71:                EPS_FEATURE_THRESHOLD=64,    /* threshold stopping test */
 72:                EPS_FEATURE_TWOSIDED=128     /* two-sided variant */
 73:              } EPSFeatureType;

 75: /*
 76:    Defines the EPS data structure
 77: */
 78: struct _p_EPS {
 79:   PETSCHEADER(struct _EPSOps);
 80:   /*------------------------- User parameters ---------------------------*/
 81:   PetscInt       max_it;           /* maximum number of iterations */
 82:   PetscInt       nev;              /* number of eigenvalues to compute */
 83:   PetscInt       ncv;              /* number of basis vectors */
 84:   PetscInt       mpd;              /* maximum dimension of projected problem */
 85:   PetscInt       nini,ninil;       /* number of initial vectors (negative means not copied yet) */
 86:   PetscInt       nds;              /* number of basis vectors of deflation space */
 87:   PetscScalar    target;           /* target value */
 88:   PetscReal      tol;              /* tolerance */
 89:   PetscReal      thres;            /* threshold */
 90:   PetscBool      threlative;       /* threshold is relative */
 91:   EPSConv        conv;             /* convergence test */
 92:   EPSStop        stop;             /* stopping test */
 93:   EPSWhich       which;            /* which part of the spectrum to be sought */
 94:   PetscReal      inta,intb;        /* interval [a,b] for spectrum slicing */
 95:   EPSProblemType problem_type;     /* which kind of problem to be solved */
 96:   EPSExtraction  extraction;       /* which kind of extraction to be applied */
 97:   EPSBalance     balance;          /* the balancing method */
 98:   PetscInt       balance_its;      /* number of iterations of the balancing method */
 99:   PetscReal      balance_cutoff;   /* cutoff value for balancing */
100:   PetscBool      trueres;          /* whether the true residual norm must be computed */
101:   PetscBool      trackall;         /* whether all the residuals must be computed */
102:   PetscBool      purify;           /* whether eigenvectors need to be purified */
103:   PetscBool      twosided;         /* whether to compute left eigenvectors (two-sided solver) */

105:   /*-------------- User-provided functions and contexts -----------------*/
106:   EPSConvergenceTestFn      *converged;
107:   EPSConvergenceTestFn      *convergeduser;
108:   PetscCtxDestroyFn         *convergeddestroy;
109:   EPSStoppingTestFn         *stopping;
110:   EPSStoppingTestFn         *stoppinguser;
111:   PetscCtxDestroyFn         *stoppingdestroy;
112:   SlepcArbitrarySelectionFn *arbitrary;
113:   void                      *convergedctx;
114:   void                      *stoppingctx;
115:   void                      *arbitraryctx;
116:   EPSMonitorFn              *monitor[MAXEPSMONITORS];
117:   PetscCtxDestroyFn         *monitordestroy[MAXEPSMONITORS];
118:   void                      *monitorcontext[MAXEPSMONITORS];
119:   PetscInt                  numbermonitors;

121:   /*----------------- Child objects and working data -------------------*/
122:   ST             st;               /* spectral transformation object */
123:   DS             ds;               /* direct solver object */
124:   BV             V;                /* set of basis vectors and computed eigenvectors */
125:   BV             W;                /* left basis vectors (if left eigenvectors requested) */
126:   RG             rg;               /* optional region for filtering */
127:   SlepcSC        sc;               /* sorting criterion data */
128:   Vec            D;                /* diagonal matrix for balancing */
129:   Vec            *IS,*ISL;         /* references to user-provided initial spaces */
130:   Vec            *defl;            /* references to user-provided deflation space */
131:   PetscScalar    *eigr,*eigi;      /* real and imaginary parts of eigenvalues */
132:   PetscReal      *errest;          /* error estimates */
133:   PetscScalar    *rr,*ri;          /* values computed by user's arbitrary selection function */
134:   PetscInt       *perm;            /* permutation for eigenvalue ordering */
135:   PetscInt       nwork;            /* number of work vectors */
136:   Vec            *work;            /* work vectors */
137:   void           *data;            /* placeholder for solver-specific stuff */

139:   /* ----------------------- Status variables --------------------------*/
140:   EPSStateType   state;            /* initial -> setup -> solved -> eigenvectors */
141:   EPSSolverType  categ;            /* solver category */
142:   PetscInt       nconv;            /* number of converged eigenvalues */
143:   PetscInt       its;              /* number of iterations so far computed */
144:   PetscInt       n,nloc;           /* problem dimensions (global, local) */
145:   PetscReal      nrma,nrmb;        /* computed matrix norms */
146:   PetscBool      useds;            /* whether the solver uses the DS object or not */
147:   PetscBool      isgeneralized;
148:   PetscBool      ispositive;
149:   PetscBool      ishermitian;
150:   PetscBool      isstructured;
151:   EPSConvergedReason reason;
152: };

154: /*
155:     Macros to test valid EPS arguments
156: */
157: #if !defined(PETSC_USE_DEBUG)

159: #define EPSCheckSolved(h,arg) do {(void)(h);} while (0)

161: #else

163: #define EPSCheckSolved(h,arg) \
164:   do { \
165:     PetscCheck((h)->state>=EPS_STATE_SOLVED,PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSolve() first: Parameter #%d",arg); \
166:   } while (0)

168: #endif

170: /*
171:     Macros to check settings at EPSSetUp()
172: */

174: /* EPSCheckHermitianDefinite: the problem is HEP or GHEP */
175: #define EPSCheckHermitianDefiniteCondition(eps,condition,msg) \
176:   do { \
177:     if (condition) { \
178:       PetscCheck((eps)->ishermitian,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for non-%s problems",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
179:       PetscCheck(!(eps)->isgeneralized || (eps)->ispositive,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires that the problem is %s-definite",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
180:     } \
181:   } while (0)
182: #define EPSCheckHermitianDefinite(eps) EPSCheckHermitianDefiniteCondition(eps,PETSC_TRUE,"")

184: /* EPSCheckHermitian: the problem is HEP, GHEP, or GHIEP */
185: #define EPSCheckHermitianCondition(eps,condition,msg) \
186:   do { \
187:     if (condition) { \
188:       PetscCheck((eps)->ishermitian,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for non-%s problems",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
189:     } \
190:   } while (0)
191: #define EPSCheckHermitian(eps) EPSCheckHermitianCondition(eps,PETSC_TRUE,"")

193: /* EPSCheckDefinite: the problem is not GHIEP */
194: #define EPSCheckDefiniteCondition(eps,condition,msg) \
195:   do { \
196:     if (condition) { \
197:       PetscCheck(!(eps)->isgeneralized || !(eps)->ishermitian || (eps)->ispositive,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for %s-indefinite problems",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
198:     } \
199:   } while (0)
200: #define EPSCheckDefinite(eps) EPSCheckDefiniteCondition(eps,PETSC_TRUE,"")

202: /* EPSCheckStandard: the problem is HEP or NHEP */
203: #define EPSCheckStandardCondition(eps,condition,msg) \
204:   do { \
205:     if (condition) { \
206:       PetscCheck(!(eps)->isgeneralized,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for generalized problems",((PetscObject)(eps))->type_name,(msg)); \
207:     } \
208:   } while (0)
209: #define EPSCheckStandard(eps) EPSCheckStandardCondition(eps,PETSC_TRUE,"")

211: /* EPSCheckNotStructured: the problem is not structured */
212: #define EPSCheckNotStructuredCondition(eps,condition,msg) \
213:   do { \
214:     if (condition) { \
215:       PetscCheck(!(eps)->isstructured,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not provide support for structured eigenproblems",((PetscObject)(eps))->type_name,(msg)); \
216:     } \
217:   } while (0)
218: #define EPSCheckNotStructured(eps) EPSCheckNotStructuredCondition(eps,PETSC_TRUE,"")

220: /* EPSCheckSinvert: shift-and-invert ST */
221: #define EPSCheckSinvertCondition(eps,condition,msg) \
222:   do { \
223:     if (condition) { \
224:       PetscBool __flg; \
225:       PetscCall(PetscObjectTypeCompare((PetscObject)(eps)->st,STSINVERT,&__flg)); \
226:       PetscCheck(__flg,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires a shift-and-invert spectral transform",((PetscObject)(eps))->type_name,(msg)); \
227:     } \
228:   } while (0)
229: #define EPSCheckSinvert(eps) EPSCheckSinvertCondition(eps,PETSC_TRUE,"")

231: /* EPSCheckSinvertCayley: shift-and-invert or Cayley ST */
232: #define EPSCheckSinvertCayleyCondition(eps,condition,msg) \
233:   do { \
234:     if (condition) { \
235:       PetscBool __flg; \
236:       PetscCall(PetscObjectTypeCompareAny((PetscObject)(eps)->st,&__flg,STSINVERT,STCAYLEY,"")); \
237:       PetscCheck(__flg,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires shift-and-invert or Cayley transform",((PetscObject)(eps))->type_name,(msg)); \
238:     } \
239:   } while (0)
240: #define EPSCheckSinvertCayley(eps) EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE,"")

242: /* Check for unsupported features */
243: #define EPSCheckUnsupportedCondition(eps,mask,condition,msg) \
244:   do { \
245:     if (condition) { \
246:       PetscCheck(!((mask) & EPS_FEATURE_BALANCE) || (eps)->balance==EPS_BALANCE_NONE,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support balancing",((PetscObject)(eps))->type_name,(msg)); \
247:       PetscCheck(!((mask) & EPS_FEATURE_ARBITRARY) || !(eps)->arbitrary,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support arbitrary selection of eigenpairs",((PetscObject)(eps))->type_name,(msg)); \
248:       if ((mask) & EPS_FEATURE_REGION) { \
249:         PetscBool      __istrivial; \
250:         PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
251:         PetscCheck(__istrivial,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support region filtering",((PetscObject)(eps))->type_name,(msg)); \
252:       } \
253:       PetscCheck(!((mask) & EPS_FEATURE_EXTRACTION) || (eps)->extraction==EPS_RITZ,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports Ritz extraction",((PetscObject)(eps))->type_name,(msg)); \
254:       PetscCheck(!((mask) & EPS_FEATURE_CONVERGENCE) || (eps)->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default convergence test",((PetscObject)(eps))->type_name,(msg)); \
255:       PetscCheck(!((mask) & EPS_FEATURE_STOPPING) || (eps)->stopping==EPSStoppingBasic,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default stopping test",((PetscObject)(eps))->type_name,(msg)); \
256:       PetscCheck(!((mask) & EPS_FEATURE_THRESHOLD) || (eps)->stopping!=EPSStoppingThreshold,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support the threshold stopping test",((PetscObject)(eps))->type_name,(msg)); \
257:       PetscCheck(!((mask) & EPS_FEATURE_TWOSIDED) || !(eps)->twosided,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot compute left eigenvectors (no two-sided variant)",((PetscObject)(eps))->type_name,(msg)); \
258:     } \
259:   } while (0)
260: #define EPSCheckUnsupported(eps,mask) EPSCheckUnsupportedCondition(eps,mask,PETSC_TRUE,"")

262: /* Check for ignored features */
263: #define EPSCheckIgnoredCondition(eps,mask,condition,msg) \
264:   do { \
265:     if (condition) { \
266:       if (((mask) & EPS_FEATURE_BALANCE) && (eps)->balance!=EPS_BALANCE_NONE) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the balancing settings\n",((PetscObject)(eps))->type_name,(msg))); \
267:       if (((mask) & EPS_FEATURE_ARBITRARY) && (eps)->arbitrary) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the settings for arbitrary selection of eigenpairs\n",((PetscObject)(eps))->type_name,(msg))); \
268:       if ((mask) & EPS_FEATURE_REGION) { \
269:         PetscBool __istrivial; \
270:         PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
271:         if (!__istrivial) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the specified region\n",((PetscObject)(eps))->type_name,(msg))); \
272:       } \
273:       if (((mask) & EPS_FEATURE_EXTRACTION) && (eps)->extraction!=EPS_RITZ) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the extraction settings\n",((PetscObject)(eps))->type_name,(msg))); \
274:       if (((mask) & EPS_FEATURE_CONVERGENCE) && (eps)->converged!=EPSConvergedRelative) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the convergence test settings\n",((PetscObject)(eps))->type_name,(msg))); \
275:       if (((mask) & EPS_FEATURE_STOPPING) && (eps)->stopping!=EPSStoppingBasic) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the stopping test settings\n",((PetscObject)(eps))->type_name,(msg))); \
276:       if (((mask) & EPS_FEATURE_TWOSIDED) && (eps)->twosided) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the two-sided flag\n",((PetscObject)(eps))->type_name,(msg))); \
277:     } \
278:   } while (0)
279: #define EPSCheckIgnored(eps,mask) EPSCheckIgnoredCondition(eps,mask,PETSC_TRUE,"")

281: /*
282:     EPSSetCtxThreshold - Fills EPSStoppingCtx with data needed for the threshold stopping test
283: */
284: #define EPSSetCtxThreshold(eps,eigr,eigi,k) \
285:   do { \
286:     if (eps->stop==EPS_STOP_THRESHOLD && k) { \
287:       PetscScalar __kr=eigr[k-1],__ki=eigi[k-1],__kr0=eigr[0],__ki0=eigi[0]; \
288:       PetscCall(STBackTransform(eps->st,1,&__kr,&__ki)); \
289:       PetscCall(STBackTransform(eps->st,1,&__kr0,&__ki0)); \
290:       if (eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_SMALLEST_MAGNITUDE) { \
291:         ((EPSStoppingCtx)eps->stoppingctx)->firstev = SlepcAbsEigenvalue(__kr0,__ki0); \
292:         ((EPSStoppingCtx)eps->stoppingctx)->lastev  = SlepcAbsEigenvalue(__kr,__ki); \
293:       } else { \
294:         ((EPSStoppingCtx)eps->stoppingctx)->firstev = PetscRealPart(__kr0); \
295:         ((EPSStoppingCtx)eps->stoppingctx)->lastev  = PetscRealPart(__kr); \
296:       } \
297:     } \
298:   } while (0)

300: /*
301:   EPS_SetInnerProduct - set B matrix for inner product if appropriate.
302: */
303: static inline PetscErrorCode EPS_SetInnerProduct(EPS eps)
304: {
305:   Mat            B;

307:   PetscFunctionBegin;
308:   if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
309:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
310:     PetscCall(STGetBilinearForm(eps->st,&B));
311:     PetscCall(BVSetMatrix(eps->V,B,PetscNot(eps->ispositive)));
312:     if (eps->twosided) PetscCall(BVSetMatrix(eps->W,B,PetscNot(eps->ispositive)));
313:     PetscCall(MatDestroy(&B));
314:   } else PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
315:   PetscFunctionReturn(PETSC_SUCCESS);
316: }

318: /*
319:   EPS_Purify - purify the first k vectors in the V basis
320: */
321: static inline PetscErrorCode EPS_Purify(EPS eps,PetscInt k)
322: {
323:   PetscInt       i;
324:   Vec            v,z;

326:   PetscFunctionBegin;
327:   PetscCall(BVCreateVec(eps->V,&v));
328:   for (i=0;i<k;i++) {
329:     PetscCall(BVCopyVec(eps->V,i,v));
330:     PetscCall(BVGetColumn(eps->V,i,&z));
331:     PetscCall(STApply(eps->st,v,z));
332:     PetscCall(BVRestoreColumn(eps->V,i,&z));
333:   }
334:   PetscCall(VecDestroy(&v));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: /*
339:   EPS_KSPSetOperators - Sets the KSP matrices, see also ST_KSPSetOperators()
340: */
341: static inline PetscErrorCode EPS_KSPSetOperators(KSP ksp,Mat A,Mat B)
342: {
343:   const char     *prefix;

345:   PetscFunctionBegin;
346:   PetscCall(KSPSetOperators(ksp,A,B));
347:   PetscCall(MatGetOptionsPrefix(B,&prefix));
348:   if (!prefix) {
349:     /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS)
350:        only applies if the Mat has no user-defined prefix */
351:     PetscCall(KSPGetOptionsPrefix(ksp,&prefix));
352:     PetscCall(MatSetOptionsPrefix(B,prefix));
353:   }
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }

357: /*
358:   EPS_GetActualConverged - Gets the actual value of nconv; in special cases the
359:   number of available eigenvalues is larger than the computed ones
360: */
361: static inline PetscErrorCode EPS_GetActualConverged(EPS eps,PetscInt *nconv)
362: {
363:   PetscFunctionBegin;
364:   *nconv = eps->nconv;
365:   if (eps->isstructured) {
366:     if (eps->problem_type == EPS_BSE && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
367:     if (eps->problem_type == EPS_HAMILT && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
368:   }
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: /*
373:   EPS_GetEigenvector - Gets the i-th eigenvector taking into account the case
374:   where i exceeds the number of computed vectors (structure-preserving solver).
375:   The argument V should be eps->V for right eigenvectors, eps->W for left ones.
376: */
377: static inline PetscErrorCode EPS_GetEigenvector(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
378: {
379:   PetscInt  k;
380:   Vec       v0,v1,w,w0,w1;
381:   Mat       H;
382:   IS        is[2];
383: #if !defined(PETSC_USE_COMPLEX)
384:   PetscInt  k0,k1,k2,iquad;
385:   PetscReal nrm,nrmr=0.0,nrmi=0.0,sgn;
386: #endif

388:   PetscFunctionBegin;
389:   if (!eps->isstructured) {
390:     k = eps->perm[i];
391:     PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
392:   } else {
393:     if (eps->problem_type == EPS_BSE && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) {
394:       /* BSE problem, even index is +lambda, odd index is -lambda */
395:       k = eps->perm[i/2];
396:       if (i%2) {
397:         /* eigenvector of -lambda is J*conj(X) where J=[0 I; I 0] and x is eigenvector of lambda */
398:         PetscCall(VecDuplicate(Vr?Vr:Vi,&w));
399:         PetscCall(STGetMatrix(eps->st,0,&H));
400:         PetscCall(MatNestGetISs(H,is,NULL));
401:         if (Vr) {
402:           PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],w,NULL));
403:           PetscCall(VecConjugate(w));
404:           PetscCall(VecGetSubVector(w,is[0],&w0));
405:           PetscCall(VecGetSubVector(w,is[1],&w1));
406:           PetscCall(VecGetSubVector(Vr,is[0],&v0));
407:           PetscCall(VecGetSubVector(Vr,is[1],&v1));
408:           PetscCall(VecCopy(w1,v0));
409:           PetscCall(VecCopy(w0,v1));
410:           PetscCall(VecRestoreSubVector(w,is[0],&w0));
411:           PetscCall(VecRestoreSubVector(w,is[1],&w1));
412:           PetscCall(VecRestoreSubVector(Vr,is[0],&v0));
413:           PetscCall(VecRestoreSubVector(Vr,is[1],&v1));
414:         }
415: #if !defined(PETSC_USE_COMPLEX)
416:         if (Vi) {
417:           PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],NULL,w));
418:           PetscCall(VecScale(w,-1.0));
419:           PetscCall(VecGetSubVector(w,is[0],&w0));
420:           PetscCall(VecGetSubVector(w,is[1],&w1));
421:           PetscCall(VecGetSubVector(Vi,is[0],&v0));
422:           PetscCall(VecGetSubVector(Vi,is[1],&v1));
423:           PetscCall(VecCopy(w1,v0));
424:           PetscCall(VecCopy(w0,v1));
425:           PetscCall(VecRestoreSubVector(w,is[0],&w0));
426:           PetscCall(VecRestoreSubVector(w,is[1],&w1));
427:           PetscCall(VecRestoreSubVector(Vi,is[0],&v0));
428:           PetscCall(VecRestoreSubVector(Vi,is[1],&v1));
429:         }
430: #endif
431:         PetscCall(VecDestroy(&w));
432:       } else {
433:         PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
434:       }
435:     } else if (eps->problem_type == EPS_HAMILT) {
436:       k = eps->perm[i/2];
437: #if !defined(PETSC_USE_COMPLEX)
438:       if (eps->eigi[k]==0.0) { /* real eigenvalue */
439:         if (Vr) {
440:           PetscCall(BVCopyVec(V,k+eps->ncv/2+1,Vr));
441:           PetscCall(BVGetColumn(V,k,&w));
442:           PetscCall(VecAXPY(Vr,(i%2)?-eps->eigr[k]:eps->eigr[k],w));
443:           PetscCall(BVRestoreColumn(V,k,&w));
444:           PetscCall(VecNorm(Vr,NORM_2,&nrmr));
445:         }
446:         if (Vi) PetscCall(VecZeroEntries(Vi));
447:         nrm = nrmr;
448:       } else if (eps->eigr[k]==0.0 ) { /* purely imaginary eigenvalue */
449:         if (Vr) {
450:           PetscCall(BVCopyVec(V,k+eps->ncv/2+1,Vr));
451:           PetscCall(VecNorm(Vr,NORM_2,&nrmr));
452:         }
453:         if (Vi) {
454:           PetscCall(BVCopyVec(V,k,Vi));
455:           PetscCall(VecScale(Vi,(i%2)?-eps->eigi[k]:eps->eigi[k]));
456:           PetscCall(VecNorm(Vi,NORM_2,&nrmi));
457:         }
458:         nrm = SlepcAbs(nrmr,nrmi);
459:       } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
460:         iquad = i%2;  /* index within the 4 values */
461:         if (i>=2) {
462:           k2 = eps->perm[(i-2)/2];
463:           if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
464:         }
465:         k0 = (iquad<2)? k: k2;
466:         k1 = k0+1;
467:         /* Vr+Vi*i obtained as eig*u+v where u=ur+ui*i is stored in cols k0 (ur) and k1 (ui) and
468:            v=vr+vi*i is in cols shifted by ncv/2+1.
469:            For lambda=eigr+eigi*i:
470:             Vr+Vi*i = (eigr+eigi*i)(ur+ui*i) + vr+vi*i
471:             Vr+Vi*i = eigr*(ur+ui*i) - eigi*ui+eigi*ur*i + vr+vi*i
472:             Vr+Vi*i = eigr*ur-eigi*ui+vr + (eigi*ur+eigr*ui+vi)*i
473:            For -conj(lambda): eigr, ui and vi have the signs changed
474:            For       -lambda: eigr and eigi have the signs changed
475:            For  conj(lambda): eigi, ui and vi have the signs changed   */
476:         if (Vr) {
477:           sgn = (iquad<2)? -1.0: 1.0;
478:           PetscCall(BVCopyVec(V,k0,Vr));                  /* ur */
479:           PetscCall(VecScale(Vr,sgn*eps->eigr[k0]));
480:           PetscCall(BVGetColumn(V,k1,&w));                /* ui */
481:           PetscCall(VecAXPY(Vr,-sgn*eps->eigi[k0],w));
482:           PetscCall(BVRestoreColumn(V,k1,&w));
483:           PetscCall(BVGetColumn(V,k0+eps->ncv/2+1,&w));   /* vr */
484:           PetscCall(VecAXPY(Vr,1.0,w));
485:           PetscCall(BVRestoreColumn(V,k0+eps->ncv/2+1,&w));
486:           PetscCall(VecNorm(Vr,NORM_2,&nrmr));
487:         }
488:         if (Vi) {
489:           sgn = (iquad%2)? -1.0: 1.0;
490:           PetscCall(BVCopyVec(V,k0,Vi));                  /* ur */
491:           PetscCall(VecScale(Vi,sgn*eps->eigi[k0]));
492:           PetscCall(BVGetColumn(V,k1,&w));                /* ui */
493:           PetscCall(VecAXPY(Vi,sgn*eps->eigr[k0],w));
494:           PetscCall(BVRestoreColumn(V,k1,&w));
495:           PetscCall(BVGetColumn(V,k1+eps->ncv/2+1,&w));   /* vi */
496:           sgn = (iquad%3)? 1.0: -1.0;
497:           PetscCall(VecAXPY(Vi,sgn,w));
498:           PetscCall(BVRestoreColumn(V,k1+eps->ncv/2+1,&w));
499:           PetscCall(VecNorm(Vi,NORM_2,&nrmi));
500:         }
501:         nrm = SlepcAbs(nrmr,nrmi);
502:       }
503:       if (Vr) PetscCall(VecScale(Vr,1.0/nrm));
504:       if (Vi) PetscCall(VecScale(Vi,1.0/nrm));
505: #endif
506:     } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Inconsistent state");
507:   }
508:   PetscFunctionReturn(PETSC_SUCCESS);
509: }

511: SLEPC_INTERN PetscErrorCode EPSSetWhichEigenpairs_Default(EPS);
512: SLEPC_INTERN PetscErrorCode EPSSetDimensions_Default(EPS,PetscInt*,PetscInt*,PetscInt*);
513: SLEPC_INTERN PetscErrorCode EPSBackTransform_Default(EPS);
514: SLEPC_INTERN PetscErrorCode EPSComputeVectors(EPS);
515: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Hermitian(EPS);
516: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Schur(EPS);
517: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Indefinite(EPS);
518: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Twosided(EPS);
519: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Slice(EPS);
520: SLEPC_INTERN PetscErrorCode EPSComputeResidualNorm_Private(EPS,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
521: SLEPC_INTERN PetscErrorCode EPSComputeRitzVector(EPS,PetscScalar*,PetscScalar*,BV,Vec,Vec);
522: SLEPC_INTERN PetscErrorCode EPSGetStartVector(EPS,PetscInt,PetscBool*);
523: SLEPC_INTERN PetscErrorCode EPSGetLeftStartVector(EPS,PetscInt,PetscBool*);
524: SLEPC_INTERN PetscErrorCode MatEstimateSpectralRange_EPS(Mat,PetscReal*,PetscReal*);

526: /* Private functions of the solver implementations */

528: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
529: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi1(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
530: SLEPC_INTERN PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscInt,PetscInt,PetscReal,PetscReal,PetscReal,PetscInt*);
531: SLEPC_INTERN PetscErrorCode EPSPseudoLanczos(EPS,PetscReal*,PetscReal*,PetscReal*,PetscInt,PetscInt*,PetscBool*,PetscBool*,PetscReal*,Vec);
532: SLEPC_INTERN PetscErrorCode EPSBuildBalance_Krylov(EPS);
533: SLEPC_INTERN PetscErrorCode EPSSetDefaultST(EPS);
534: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_Precond(EPS);
535: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_GMRES(EPS);
536: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_NoFactor(EPS);
537: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Basic(EPS);
538: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Default(EPS);