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:   PetscCtxDestroyFn         *arbitrarydestroy;
    114:   void                      *convergedctx;
    115:   void                      *stoppingctx;
    116:   void                      *arbitraryctx;
    117:   EPSMonitorFn              *monitor[MAXEPSMONITORS];
    118:   PetscCtxDestroyFn         *monitordestroy[MAXEPSMONITORS];
    119:   void                      *monitorcontext[MAXEPSMONITORS];
    120:   PetscInt                  numbermonitors;

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

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

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

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

    162: #else

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

    169: #endif

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

    175: /* EPSCheckHermitianDefinite: the problem is HEP or GHEP */
    176: #define EPSCheckHermitianDefiniteCondition(eps,condition,msg) \
    177:   do { \
    178:     if (condition) { \
    179:       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); \
    180:       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); \
    181:     } \
    182:   } while (0)
    183: #define EPSCheckHermitianDefinite(eps) EPSCheckHermitianDefiniteCondition(eps,PETSC_TRUE,"")

    185: /* EPSCheckHermitian: the problem is HEP, GHEP, or GHIEP */
    186: #define EPSCheckHermitianCondition(eps,condition,msg) \
    187:   do { \
    188:     if (condition) { \
    189:       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); \
    190:     } \
    191:   } while (0)
    192: #define EPSCheckHermitian(eps) EPSCheckHermitianCondition(eps,PETSC_TRUE,"")

    194: /* EPSCheckDefinite: the problem is not GHIEP */
    195: #define EPSCheckDefiniteCondition(eps,condition,msg) \
    196:   do { \
    197:     if (condition) { \
    198:       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); \
    199:     } \
    200:   } while (0)
    201: #define EPSCheckDefinite(eps) EPSCheckDefiniteCondition(eps,PETSC_TRUE,"")

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

    212: /* EPSCheckNotStructured: the problem is not structured */
    213: #define EPSCheckNotStructuredCondition(eps,condition,msg) \
    214:   do { \
    215:     if (condition) { \
    216:       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)); \
    217:     } \
    218:   } while (0)
    219: #define EPSCheckNotStructured(eps) EPSCheckNotStructuredCondition(eps,PETSC_TRUE,"")

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

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

    243: /* Check for unsupported features */
    244: #define EPSCheckUnsupportedCondition(eps,mask,condition,msg) \
    245:   do { \
    246:     if (condition) { \
    247:       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)); \
    248:       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)); \
    249:       if ((mask) & EPS_FEATURE_REGION) { \
    250:         PetscBool      __istrivial; \
    251:         PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
    252:         PetscCheck(__istrivial,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support region filtering",((PetscObject)(eps))->type_name,(msg)); \
    253:       } \
    254:       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)); \
    255:       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)); \
    256:       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)); \
    257:       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)); \
    258:       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)); \
    259:     } \
    260:   } while (0)
    261: #define EPSCheckUnsupported(eps,mask) EPSCheckUnsupportedCondition(eps,mask,PETSC_TRUE,"")

    263: /* Check for ignored features */
    264: #define EPSCheckIgnoredCondition(eps,mask,condition,msg) \
    265:   do { \
    266:     if (condition) { \
    267:       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))); \
    268:       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))); \
    269:       if ((mask) & EPS_FEATURE_REGION) { \
    270:         PetscBool __istrivial; \
    271:         PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
    272:         if (!__istrivial) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the specified region\n",((PetscObject)(eps))->type_name,(msg))); \
    273:       } \
    274:       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))); \
    275:       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))); \
    276:       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))); \
    277:       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))); \
    278:     } \
    279:   } while (0)
    280: #define EPSCheckIgnored(eps,mask) EPSCheckIgnoredCondition(eps,mask,PETSC_TRUE,"")

    282: /*
    283:     EPSSetCtxThreshold - Fills EPSStoppingCtx with data needed for the threshold stopping test

    285:     k = number of converged approximations, n = total number of available approximations
    286: */
    287: #define EPSSetCtxThreshold(eps,eigr,eigi,err_est,k,n) \
    288:   do { \
    289:     if ((eps)->stop==EPS_STOP_THRESHOLD && k) { \
    290:       PetscScalar __kr=(eigr)[k-1],__ki=(eigi)[k-1],__krn=0.0,__kin=0.0,__kr0=(eigr)[0],__ki0=(eigi)[0]; \
    291:       PetscCall(STBackTransform((eps)->st,1,&__kr,&__ki)); \
    292:       PetscCall(STBackTransform((eps)->st,1,&__kr0,&__ki0)); \
    293:       if (n>k) { \
    294:         __krn=(eigr)[k];__kin=(eigi)[k]; \
    295:         PetscCall(STBackTransform((eps)->st,1,&__krn,&__kin)); \
    296:       } \
    297:       if ((eps)->which==EPS_LARGEST_MAGNITUDE || (eps)->which==EPS_SMALLEST_MAGNITUDE) { \
    298:         ((EPSStoppingCtx)(eps)->stoppingctx)->firstev = SlepcAbsEigenvalue(__kr0,__ki0); \
    299:         ((EPSStoppingCtx)(eps)->stoppingctx)->lastev  = SlepcAbsEigenvalue(__kr,__ki); \
    300:         ((EPSStoppingCtx)(eps)->stoppingctx)->firstnc = SlepcAbsEigenvalue(__krn,__kin); \
    301:       } else { \
    302:         ((EPSStoppingCtx)(eps)->stoppingctx)->firstev = PetscRealPart(__kr0); \
    303:         ((EPSStoppingCtx)(eps)->stoppingctx)->lastev  = PetscRealPart(__kr); \
    304:         ((EPSStoppingCtx)(eps)->stoppingctx)->firstnc = PetscRealPart(__krn); \
    305:       } \
    306:       ((EPSStoppingCtx)(eps)->stoppingctx)->errest  = (err_est)[k]; \
    307:       ((EPSStoppingCtx)(eps)->stoppingctx)->napprox = n; \
    308:     } \
    309:   } while (0)

    311: /*
    312:   EPS_SetInnerProduct - set B matrix for inner product if appropriate.
    313: */
    314: static inline PetscErrorCode EPS_SetInnerProduct(EPS eps)
    315: {
    316:   Mat            B;

    318:   PetscFunctionBegin;
    319:   if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
    320:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
    321:     PetscCall(STGetBilinearForm(eps->st,&B));
    322:     PetscCall(BVSetMatrix(eps->V,B,PetscNot(eps->ispositive)));
    323:     if (eps->twosided) PetscCall(BVSetMatrix(eps->W,B,PetscNot(eps->ispositive)));
    324:     PetscCall(MatDestroy(&B));
    325:   } else PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: /*
    330:   EPS_Purify - purify the first k vectors in the V basis
    331: */
    332: static inline PetscErrorCode EPS_Purify(EPS eps,PetscInt k)
    333: {
    334:   PetscInt       i;
    335:   Vec            v,z;

    337:   PetscFunctionBegin;
    338:   PetscCall(BVCreateVec(eps->V,&v));
    339:   for (i=0;i<k;i++) {
    340:     PetscCall(BVCopyVec(eps->V,i,v));
    341:     PetscCall(BVGetColumn(eps->V,i,&z));
    342:     PetscCall(STApply(eps->st,v,z));
    343:     PetscCall(BVRestoreColumn(eps->V,i,&z));
    344:   }
    345:   PetscCall(VecDestroy(&v));
    346:   PetscFunctionReturn(PETSC_SUCCESS);
    347: }

    349: /*
    350:   EPS_KSPSetOperators - Sets the KSP matrices, see also ST_KSPSetOperators()
    351: */
    352: static inline PetscErrorCode EPS_KSPSetOperators(KSP ksp,Mat A,Mat B)
    353: {
    354:   const char     *prefix;

    356:   PetscFunctionBegin;
    357:   PetscCall(KSPSetOperators(ksp,A,B));
    358:   PetscCall(MatGetOptionsPrefix(B,&prefix));
    359:   if (!prefix) {
    360:     /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS)
    361:        only applies if the Mat has no user-defined prefix */
    362:     PetscCall(KSPGetOptionsPrefix(ksp,&prefix));
    363:     PetscCall(MatSetOptionsPrefix(B,prefix));
    364:   }
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }

    368: /*
    369:   EPS_GetActualConverged - Gets the actual value of nconv; in special cases the
    370:   number of available eigenvalues is larger than the computed ones
    371: */
    372: static inline PetscErrorCode EPS_GetActualConverged(EPS eps,PetscInt *nconv)
    373: {
    374:   PetscFunctionBegin;
    375:   *nconv = eps->nconv;
    376:   if (eps->isstructured) {
    377:     if (eps->problem_type == EPS_BSE && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
    378:     else if (eps->problem_type == EPS_HAMILT && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
    379:     else if (eps->problem_type == EPS_LREP && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
    380:   }
    381:   PetscFunctionReturn(PETSC_SUCCESS);
    382: }

    384: static inline PetscErrorCode EPS_GetEigenvector_BSE(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
    385: {
    386:   PetscInt  k;
    387:   Vec       v0,v1,w,w0,w1;
    388:   Mat       H;
    389:   IS        is[2];

    391:   PetscFunctionBegin;
    392:   PetscCheck(eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_PLIB,"Inconsistent state");
    393:   /* BSE problem, even index is +lambda, odd index is -lambda */
    394:   k = eps->perm[i/2];
    395:   if (i%2) {
    396:     /* eigenvector of -lambda is J*conj(x) where J=[0 I; I 0] and x is eigenvector of lambda */
    397:     PetscCall(VecDuplicate(Vr?Vr:Vi,&w));
    398:     PetscCall(STGetMatrix(eps->st,0,&H));
    399:     PetscCall(MatNestGetISs(H,is,NULL));
    400:     if (Vr) {
    401:       PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],w,NULL));
    402:       PetscCall(VecConjugate(w));
    403:       PetscCall(VecGetSubVector(w,is[0],&w0));
    404:       PetscCall(VecGetSubVector(w,is[1],&w1));
    405:       PetscCall(VecGetSubVector(Vr,is[0],&v0));
    406:       PetscCall(VecGetSubVector(Vr,is[1],&v1));
    407:       PetscCall(VecCopy(w1,v0));
    408:       PetscCall(VecCopy(w0,v1));
    409:       PetscCall(VecRestoreSubVector(w,is[0],&w0));
    410:       PetscCall(VecRestoreSubVector(w,is[1],&w1));
    411:       PetscCall(VecRestoreSubVector(Vr,is[0],&v0));
    412:       PetscCall(VecRestoreSubVector(Vr,is[1],&v1));
    413:     }
    414: #if !defined(PETSC_USE_COMPLEX)
    415:     if (Vi) {
    416:       PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],NULL,w));
    417:       PetscCall(VecScale(w,-1.0));
    418:       PetscCall(VecGetSubVector(w,is[0],&w0));
    419:       PetscCall(VecGetSubVector(w,is[1],&w1));
    420:       PetscCall(VecGetSubVector(Vi,is[0],&v0));
    421:       PetscCall(VecGetSubVector(Vi,is[1],&v1));
    422:       PetscCall(VecCopy(w1,v0));
    423:       PetscCall(VecCopy(w0,v1));
    424:       PetscCall(VecRestoreSubVector(w,is[0],&w0));
    425:       PetscCall(VecRestoreSubVector(w,is[1],&w1));
    426:       PetscCall(VecRestoreSubVector(Vi,is[0],&v0));
    427:       PetscCall(VecRestoreSubVector(Vi,is[1],&v1));
    428:     }
    429: #endif
    430:     PetscCall(VecDestroy(&w));
    431:   } else {
    432:     PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
    433:   }
    434:   PetscFunctionReturn(PETSC_SUCCESS);
    435: }

    437: static inline PetscErrorCode EPS_GetEigenvector_HAMILT(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
    438: {
    439: #if !defined(PETSC_USE_COMPLEX)
    440:   PetscInt  k;
    441:   Vec       w;
    442:   PetscInt  k0,k1,k2,iquad;
    443:   PetscReal nrm,nrmr=0.0,nrmi=0.0,sgn;
    444: #endif

    446:   PetscFunctionBegin;
    447: #if !defined(PETSC_USE_COMPLEX)
    448:   k = eps->perm[i/2];
    449:   if (eps->eigi[k]==0.0) { /* real eigenvalue */
    450:     if (Vr) {
    451:       PetscCall(BVCopyVec(V,k+eps->ncv/2+1,Vr));
    452:       PetscCall(BVGetColumn(V,k,&w));
    453:       PetscCall(VecAXPY(Vr,(i%2)?-eps->eigr[k]:eps->eigr[k],w));
    454:       PetscCall(BVRestoreColumn(V,k,&w));
    455:       PetscCall(VecNorm(Vr,NORM_2,&nrmr));
    456:     }
    457:     if (Vi) PetscCall(VecZeroEntries(Vi));
    458:     nrm = nrmr;
    459:   } else if (eps->eigr[k]==0.0 ) { /* purely imaginary eigenvalue */
    460:     if (Vr) {
    461:       PetscCall(BVCopyVec(V,k+eps->ncv/2+1,Vr));
    462:       PetscCall(VecNorm(Vr,NORM_2,&nrmr));
    463:     }
    464:     if (Vi) {
    465:       PetscCall(BVCopyVec(V,k,Vi));
    466:       PetscCall(VecScale(Vi,(i%2)?-eps->eigi[k]:eps->eigi[k]));
    467:       PetscCall(VecNorm(Vi,NORM_2,&nrmi));
    468:     }
    469:     nrm = SlepcAbs(nrmr,nrmi);
    470:   } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
    471:     iquad = i%2;  /* index within the 4 values */
    472:     if (i>=2) {
    473:       k2 = eps->perm[(i-2)/2];
    474:       if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
    475:     }
    476:     k0 = (iquad<2)? k: k2;
    477:     k1 = k0+1;
    478:     /* Vr+Vi*i obtained as eig*u+v where u=ur+ui*i is stored in cols k0 (ur) and k1 (ui) and
    479:        v=vr+vi*i is in cols shifted by ncv/2+1.
    480:        For lambda=eigr+eigi*i:
    481:         Vr+Vi*i = (eigr+eigi*i)(ur+ui*i) + vr+vi*i
    482:         Vr+Vi*i = eigr*(ur+ui*i) - eigi*ui+eigi*ur*i + vr+vi*i
    483:         Vr+Vi*i = eigr*ur-eigi*ui+vr + (eigi*ur+eigr*ui+vi)*i
    484:        For -conj(lambda): eigr, ui and vi have the signs changed
    485:        For       -lambda: eigr and eigi have the signs changed
    486:        For  conj(lambda): eigi, ui and vi have the signs changed   */
    487:     if (Vr) {
    488:       sgn = (iquad<2)? -1.0: 1.0;
    489:       PetscCall(BVCopyVec(V,k0,Vr));                  /* ur */
    490:       PetscCall(VecScale(Vr,sgn*eps->eigr[k0]));
    491:       PetscCall(BVGetColumn(V,k1,&w));                /* ui */
    492:       PetscCall(VecAXPY(Vr,-sgn*eps->eigi[k0],w));
    493:       PetscCall(BVRestoreColumn(V,k1,&w));
    494:       PetscCall(BVGetColumn(V,k0+eps->ncv/2+1,&w));   /* vr */
    495:       PetscCall(VecAXPY(Vr,1.0,w));
    496:       PetscCall(BVRestoreColumn(V,k0+eps->ncv/2+1,&w));
    497:       PetscCall(VecNorm(Vr,NORM_2,&nrmr));
    498:     }
    499:     if (Vi) {
    500:       sgn = (iquad%2)? -1.0: 1.0;
    501:       PetscCall(BVCopyVec(V,k0,Vi));                  /* ur */
    502:       PetscCall(VecScale(Vi,sgn*eps->eigi[k0]));
    503:       PetscCall(BVGetColumn(V,k1,&w));                /* ui */
    504:       PetscCall(VecAXPY(Vi,sgn*eps->eigr[k0],w));
    505:       PetscCall(BVRestoreColumn(V,k1,&w));
    506:       PetscCall(BVGetColumn(V,k1+eps->ncv/2+1,&w));   /* vi */
    507:       sgn = (iquad%3)? 1.0: -1.0;
    508:       PetscCall(VecAXPY(Vi,sgn,w));
    509:       PetscCall(BVRestoreColumn(V,k1+eps->ncv/2+1,&w));
    510:       PetscCall(VecNorm(Vi,NORM_2,&nrmi));
    511:     }
    512:     nrm = SlepcAbs(nrmr,nrmi);
    513:   }
    514:   if (Vr) PetscCall(VecScale(Vr,1.0/nrm));
    515:   if (Vi) PetscCall(VecScale(Vi,1.0/nrm));
    516: #endif
    517:   PetscFunctionReturn(PETSC_SUCCESS);
    518: }

    520: static inline PetscErrorCode EPS_GetEigenvector_LREP(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
    521: {
    522:   PetscInt  k;
    523:   Vec       v;
    524:   Mat       H;
    525:   IS        is[2];

    527:   PetscFunctionBegin;
    528:   /* LREP problem, even index is +lambda, odd index is -lambda */
    529:   k = eps->perm[i/2];
    530:   PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
    531:   if (i%2) {
    532:     /* eigenvector of -lambda is S*x where S=[I 0; 0 -I] and x is eigenvector of lambda */
    533:     PetscCall(STGetMatrix(eps->st,0,&H));
    534:     PetscCall(MatNestGetISs(H,is,NULL));
    535:     PetscCall(VecGetSubVector(Vr,is[1],&v));
    536:     PetscCall(VecScale(v,-1.0));
    537:     PetscCall(VecRestoreSubVector(Vr,is[1],&v));
    538:   }
    539:   PetscFunctionReturn(PETSC_SUCCESS);
    540: }

    542: /*
    543:   EPS_GetEigenvector - Gets the i-th eigenvector taking into account the case
    544:   where i exceeds the number of computed vectors (structure-preserving solver).
    545:   The argument V should be eps->V for right eigenvectors, eps->W for left ones.
    546: */
    547: static inline PetscErrorCode EPS_GetEigenvector(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
    548: {
    549:   PetscInt  k;
    550:   PetscBool reduced;
    551:   Mat       H;

    553:   PetscFunctionBegin;
    554:   if (!eps->isstructured) {
    555:     k = eps->perm[i];
    556:     PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
    557:   } else {
    558:     switch (eps->problem_type) {
    559:       case EPS_BSE:
    560:         PetscCall(EPS_GetEigenvector_BSE(eps,V,i,Vr,Vi));
    561:         break;
    562:       case EPS_HAMILT:
    563:         PetscCall(EPS_GetEigenvector_HAMILT(eps,V,i,Vr,Vi));
    564:         break;
    565:       case EPS_LREP:
    566:         PetscCall(STGetMatrix(eps->st,0,&H));
    567:         PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
    568:         if (reduced) PetscCall(EPS_GetEigenvector_LREP(eps,V,i,Vr,Vi));
    569:         else PetscCall(EPS_GetEigenvector_BSE(eps,V,i,Vr,Vi));
    570:         break;
    571:       default:
    572:         SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Inconsistent state");
    573:     }
    574:   }
    575:   PetscFunctionReturn(PETSC_SUCCESS);
    576: }

    578: SLEPC_INTERN PetscErrorCode EPSSetWhichEigenpairs_Default(EPS);
    579: SLEPC_INTERN PetscErrorCode EPSSetDimensions_Default(EPS,PetscInt*,PetscInt*,PetscInt*);
    580: SLEPC_INTERN PetscErrorCode EPSBackTransform_Default(EPS);
    581: SLEPC_INTERN PetscErrorCode EPSComputeVectors(EPS);
    582: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Hermitian(EPS);
    583: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Schur(EPS);
    584: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Indefinite(EPS);
    585: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Twosided(EPS);
    586: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Slice(EPS);
    587: SLEPC_INTERN PetscErrorCode EPSComputeResidualNorm_Private(EPS,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
    588: SLEPC_INTERN PetscErrorCode EPSComputeRitzVector(EPS,PetscScalar*,PetscScalar*,BV,Vec,Vec);
    589: SLEPC_INTERN PetscErrorCode EPSGetStartVector(EPS,PetscInt,PetscBool*);
    590: SLEPC_INTERN PetscErrorCode EPSGetLeftStartVector(EPS,PetscInt,PetscBool*);
    591: SLEPC_INTERN PetscErrorCode MatEstimateSpectralRange_EPS(Mat,PetscReal*,PetscReal*);

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

    595: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
    596: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi1(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
    597: SLEPC_INTERN PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscInt,PetscInt,PetscReal,PetscReal,PetscReal,PetscInt*);
    598: SLEPC_INTERN PetscErrorCode EPSPseudoLanczos(EPS,PetscReal*,PetscReal*,PetscReal*,PetscInt,PetscInt*,PetscBool*,PetscBool*,PetscReal*,Vec);
    599: SLEPC_INTERN PetscErrorCode EPSBuildBalance_Krylov(EPS);
    600: SLEPC_INTERN PetscErrorCode EPSSetDefaultST(EPS);
    601: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_Precond(EPS);
    602: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_GMRES(EPS);
    603: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_NoFactor(EPS);
    604: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Basic(EPS);
    605: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Default(EPS);