Actual source code: epssolve.c

      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: */
     10: /*
     11:    EPS routines related to the solution process
     12: */

     14: #include <slepc/private/epsimpl.h>
     15: #include <slepc/private/bvimpl.h>
     16: #include <petscdraw.h>

     18: PetscErrorCode EPSComputeVectors(EPS eps)
     19: {
     20:   PetscFunctionBegin;
     21:   EPSCheckSolved(eps,1);
     22:   if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
     23:   eps->state = EPS_STATE_EIGENVECTORS;
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

     27: static PetscErrorCode EPSComputeValues(EPS eps)
     28: {
     29:   PetscBool      injective,iscomp,isfilter;
     30:   PetscInt       i,n,aux,nconv0;
     31:   Mat            A,B=NULL,G,Z;

     33:   PetscFunctionBegin;
     34:   switch (eps->categ) {
     35:     case EPS_CATEGORY_KRYLOV:
     36:     case EPS_CATEGORY_OTHER:
     37:       PetscCall(STIsInjective(eps->st,&injective));
     38:       if (injective) {
     39:         /* one-to-one mapping: backtransform eigenvalues */
     40:         PetscUseTypeMethod(eps,backtransform);
     41:       } else {
     42:         /* compute eigenvalues from Rayleigh quotient */
     43:         PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
     44:         if (!n) break;
     45:         PetscCall(EPSGetOperators(eps,&A,&B));
     46:         PetscCall(BVSetActiveColumns(eps->V,0,n));
     47:         PetscCall(DSGetCompact(eps->ds,&iscomp));
     48:         PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
     49:         PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
     50:         PetscCall(BVMatProject(eps->V,A,eps->V,G));
     51:         PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
     52:         if (B) {
     53:           PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
     54:           PetscCall(BVMatProject(eps->V,B,eps->V,G));
     55:           PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
     56:         }
     57:         PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     58:         PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     59:         PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     60:         PetscCall(DSSetCompact(eps->ds,iscomp));
     61:         if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
     62:           PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
     63:           PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
     64:           PetscCall(BVMultInPlace(eps->V,Z,0,n));
     65:           PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
     66:         }
     67:         /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
     68:         PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
     69:         if (isfilter) {
     70:           nconv0 = eps->nconv;
     71:           for (i=0;i<eps->nconv;i++) {
     72:             if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
     73:               eps->nconv--;
     74:               if (i<eps->nconv) { SlepcSwap(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
     75:             }
     76:           }
     77:           if (nconv0>eps->nconv) PetscCall(PetscInfo(eps,"Discarded %" PetscInt_FMT " computed eigenvalues lying outside the interval\n",nconv0-eps->nconv));
     78:         }
     79:       }
     80:       break;
     81:     case EPS_CATEGORY_PRECOND:
     82:     case EPS_CATEGORY_CONTOUR:
     83:       /* eigenvalues already available as an output of the solver */
     84:       break;
     85:   }
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: /*@
     90:    EPSSolve - Solves the eigensystem.

     92:    Collective

     94:    Input Parameter:
     95: .  eps - the linear eigensolver context

     97:    Options Database Keys:
     98: +  -eps_view             - print information about the solver once the solve is complete
     99: .  -eps_view_pre         - print information about the solver before the solve starts
    100: .  -eps_view_mat0        - view the first matrix ($A$)
    101: .  -eps_view_mat1        - view the second matrix ($B$)
    102: .  -eps_view_vectors     - view the computed eigenvectors
    103: .  -eps_view_values      - view the computed eigenvalues
    104: .  -eps_converged_reason - print reason for convergence/divergence, and number of iterations
    105: .  -eps_error_absolute   - print absolute errors of each eigenpair
    106: .  -eps_error_relative   - print relative errors of each eigenpair
    107: -  -eps_error_backward   - print backward errors of each eigenpair

    109:    Notes:
    110:    The problem matrices are specified with `EPSSetOperators()`.

    112:    `EPSSolve()` will return without generating an error regardless of whether
    113:    all requested solutions were computed or not. Call `EPSGetConverged()` to get the
    114:    actual number of computed solutions, and `EPSGetConvergedReason()` to determine if
    115:    the solver converged or failed and why.

    117:    All the command-line options listed above admit an optional argument specifying
    118:    the viewer type and options. For instance, use `-eps_view_mat0 binary:amatrix.bin`
    119:    to save the $A$ matrix to a binary file, `-eps_view_values draw` to draw the computed
    120:    eigenvalues graphically, or `-eps_error_relative :myerr.m:ascii_matlab` to save
    121:    the errors in a file that can be executed in Matlab.
    122:    See `PetscObjectViewFromOptions()` for more details.

    124:    Level: beginner

    126: .seealso: [](ch:eps), `EPSCreate()`, `EPSSetUp()`, `EPSDestroy()`, `EPSSetOperators()`, `EPSGetConverged()`, `EPSGetConvergedReason()`
    127: @*/
    128: PetscErrorCode EPSSolve(EPS eps)
    129: {
    130:   PetscInt       i;
    131:   PetscBool      hasname;
    132:   STMatMode      matmode;
    133:   Mat            A,B;

    135:   PetscFunctionBegin;
    137:   if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
    138:   PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));

    140:   /* Call setup */
    141:   PetscCall(EPSSetUp(eps));

    143:   /* Safeguard for matrices of size 0 */
    144:   if (eps->n == 0) {
    145:     eps->nconv  = 0;
    146:     eps->reason = EPS_CONVERGED_TOL;
    147:     eps->state  = EPS_STATE_SOLVED;
    148:     PetscFunctionReturn(PETSC_SUCCESS);
    149:   }

    151:   eps->nconv = 0;
    152:   eps->its   = 0;
    153:   for (i=0;i<eps->ncv;i++) {
    154:     eps->eigr[i]   = 0.0;
    155:     eps->eigi[i]   = 0.0;
    156:     eps->errest[i] = 0.0;
    157:     eps->perm[i]   = i;
    158:   }
    159:   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
    160:   PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));

    162:   /* Call solver */
    163:   PetscUseTypeMethod(eps,solve);
    164:   PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
    165:   eps->state = EPS_STATE_SOLVED;

    167:   /* Only the first nconv columns contain useful information (except in CISS) */
    168:   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
    169:   if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));

    171:   /* If inplace, purify eigenvectors before reverting operator */
    172:   PetscCall(STGetMatMode(eps->st,&matmode));
    173:   if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
    174:   PetscCall(STPostSolve(eps->st));

    176:   /* Map eigenvalues back to the original problem if appropriate */
    177:   PetscCall(EPSComputeValues(eps));

    179: #if !defined(PETSC_USE_COMPLEX)
    180:   /* Reorder conjugate eigenvalues (positive imaginary first) */
    181:   for (i=0;i<eps->nconv-1;i++) {
    182:     if (eps->eigi[i] != 0 && (eps->problem_type!=EPS_HAMILT || eps->eigr[i]!=0)) {
    183:       /* conjugate eigenvalues */
    184:       if (eps->eigi[i] < 0) {
    185:         eps->eigi[i] = -eps->eigi[i];
    186:         eps->eigi[i+1] = -eps->eigi[i+1];
    187:         /* the next correction only works with eigenvectors */
    188:         PetscCall(EPSComputeVectors(eps));
    189:         PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
    190:         if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
    191:       }
    192:       i++;
    193:     }
    194:   }
    195: #endif

    197:   /* Sort eigenvalues according to eps->which parameter */
    198:   if (eps->problem_type==EPS_HAMILT) PetscCall(SlepcSortEigenvaluesSpecial(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
    199:   else PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
    200:   PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));

    202:   /* Various viewers */
    203:   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
    204:   PetscCall(EPSConvergedReasonViewFromOptions(eps));
    205:   PetscCall(EPSErrorViewFromOptions(eps));
    206:   PetscCall(EPSValuesViewFromOptions(eps));
    207:   PetscCall(EPSVectorsViewFromOptions(eps));

    209:   PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
    210:   if (hasname) {
    211:     PetscCall(EPSGetOperators(eps,&A,NULL));
    212:     PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
    213:   }
    214:   if (eps->isgeneralized) {
    215:     PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
    216:     if (hasname) {
    217:       PetscCall(EPSGetOperators(eps,NULL,&B));
    218:       PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
    219:     }
    220:   }

    222:   /* Remove deflation and initial subspaces */
    223:   if (eps->nds) {
    224:     PetscCall(BVSetNumConstraints(eps->V,0));
    225:     eps->nds = 0;
    226:   }
    227:   eps->nini = 0;
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: /*@
    232:    EPSGetIterationNumber - Gets the current iteration number. If the
    233:    call to `EPSSolve()` is complete, then it returns the number of iterations
    234:    carried out by the solution method.

    236:    Not Collective

    238:    Input Parameter:
    239: .  eps - the linear eigensolver context

    241:    Output Parameter:
    242: .  its - number of iterations

    244:    Note:
    245:    During the $i$-th iteration this call returns $i-1$. If `EPSSolve()` is
    246:    complete, then parameter `its` contains either the iteration number at
    247:    which convergence was successfully reached, or failure was detected.
    248:    Call `EPSGetConvergedReason()` to determine if the solver converged or
    249:    failed and why.

    251:    Level: intermediate

    253: .seealso: [](ch:eps), `EPSGetConvergedReason()`, `EPSSetTolerances()`
    254: @*/
    255: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
    256: {
    257:   PetscFunctionBegin;
    259:   PetscAssertPointer(its,2);
    260:   *its = eps->its;
    261:   PetscFunctionReturn(PETSC_SUCCESS);
    262: }

    264: /*@
    265:    EPSGetConverged - Gets the number of converged eigenpairs.

    267:    Not Collective

    269:    Input Parameter:
    270: .  eps - the linear eigensolver context

    272:    Output Parameter:
    273: .  nconv - number of converged eigenpairs

    275:    Notes:
    276:    This function should be called after `EPSSolve()` has finished.

    278:    The value `nconv` may be different from the number of requested solutions
    279:    `nev`, but not larger than `ncv`, see `EPSSetDimensions()`.

    281:    Level: beginner

    283: .seealso: [](ch:eps), `EPSSetDimensions()`, `EPSSolve()`, `EPSGetEigenpair()`
    284: @*/
    285: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
    286: {
    287:   PetscFunctionBegin;
    289:   PetscAssertPointer(nconv,2);
    290:   EPSCheckSolved(eps,1);
    291:   PetscCall(EPS_GetActualConverged(eps,nconv));
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

    295: /*@
    296:    EPSGetConvergedReason - Gets the reason why the `EPSSolve()` iteration was
    297:    stopped.

    299:    Not Collective

    301:    Input Parameter:
    302: .  eps - the linear eigensolver context

    304:    Output Parameter:
    305: .  reason - negative value indicates diverged, positive value converged, see
    306:             `EPSConvergedReason` for the possible values

    308:    Options Database Key:
    309: .  -eps_converged_reason - print reason for convergence/divergence, and number of iterations

    311:    Note:
    312:    If this routine is called before or doing the `EPSSolve()` the value of
    313:    `EPS_CONVERGED_ITERATING` is returned.

    315:    Level: intermediate

    317: .seealso: [](ch:eps), `EPSSetTolerances()`, `EPSSolve()`, `EPSConvergedReason`
    318: @*/
    319: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
    320: {
    321:   PetscFunctionBegin;
    323:   PetscAssertPointer(reason,2);
    324:   EPSCheckSolved(eps,1);
    325:   *reason = eps->reason;
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: /*@
    330:    EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
    331:    subspace.

    333:    Collective

    335:    Input Parameter:
    336: .  eps - the linear eigensolver context

    338:    Output Parameter:
    339: .  v - an array of vectors

    341:    Notes:
    342:    This function should be called after `EPSSolve()` has finished.

    344:    The user should provide in `v` an array of `nconv` vectors, where `nconv`
    345:    is the value returned by `EPSGetConverged()`.

    347:    The first $k$ vectors returned in `v` span an invariant subspace associated
    348:    with the first $k$ computed eigenvalues (note that this is not true if the
    349:    $k$-th eigenvalue is complex and matrix $A$ is real; in this case the first
    350:    $k+1$ vectors should be used). An invariant subspace $X$ of $A$ satisfies
    351:    $Ax\in X, \forall x\in X$ (a similar definition applies for generalized
    352:    eigenproblems).

    354:    Level: intermediate

    356: .seealso: [](ch:eps), `EPSGetEigenpair()`, `EPSGetConverged()`, `EPSSolve()`
    357: @*/
    358: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
    359: {
    360:   PetscInt       i;
    361:   BV             V=eps->V;
    362:   Vec            w;

    364:   PetscFunctionBegin;
    366:   PetscAssertPointer(v,2);
    368:   EPSCheckSolved(eps,1);
    369:   PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
    370:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
    371:     PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
    372:     PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
    373:     PetscCall(BVCopy(eps->V,V));
    374:     for (i=0;i<eps->nconv;i++) {
    375:       PetscCall(BVGetColumn(V,i,&w));
    376:       PetscCall(VecPointwiseDivide(w,w,eps->D));
    377:       PetscCall(BVRestoreColumn(V,i,&w));
    378:     }
    379:     PetscCall(BVOrthogonalize(V,NULL));
    380:   }
    381:   for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
    382:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
    383:   PetscFunctionReturn(PETSC_SUCCESS);
    384: }

    386: /*@
    387:    EPSGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
    388:    `EPSSolve()`. The solution consists in both the eigenvalue and the eigenvector.

    390:    Collective

    392:    Input Parameters:
    393: +  eps - the linear eigensolver context
    394: -  i   - index of the solution

    396:    Output Parameters:
    397: +  eigr - real part of eigenvalue
    398: .  eigi - imaginary part of eigenvalue
    399: .  Vr   - real part of eigenvector
    400: -  Vi   - imaginary part of eigenvector

    402:    Notes:
    403:    It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
    404:    required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
    405:    they must be created by the calling program with e.g. `MatCreateVecs()`.

    407:    If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
    408:    configured with complex scalars the eigenvalue is stored
    409:    directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi` is
    410:    set to zero). In both cases, the user can pass `NULL` in `eigi` and `Vi`.

    412:    The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
    413:    Eigenpairs are indexed according to the ordering criterion established
    414:    with `EPSSetWhichEigenpairs()`.

    416:    The 2-norm of the eigenvector is one unless the problem is generalized
    417:    Hermitian. In this case the eigenvector is normalized with respect to the
    418:    norm defined by the $B$ matrix.

    420:    In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
    421:    [](#sec:structured-vectors).

    423:    Level: beginner

    425: .seealso: [](ch:eps), `EPSGetEigenvalue()`, `EPSGetEigenvector()`, `EPSGetLeftEigenvector()`, `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetInvariantSubspace()`
    426: @*/
    427: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
    428: {
    429:   PetscInt nconv;

    431:   PetscFunctionBegin;
    434:   EPSCheckSolved(eps,1);
    435:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    436:   PetscCall(EPS_GetActualConverged(eps,&nconv));
    437:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
    438:   PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
    439:   if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
    440:   PetscFunctionReturn(PETSC_SUCCESS);
    441: }

    443: /*@
    444:    EPSGetEigenvalue - Gets the `i`-th eigenvalue as computed by `EPSSolve()`.

    446:    Not Collective

    448:    Input Parameters:
    449: +  eps - the linear eigensolver context
    450: -  i   - index of the solution

    452:    Output Parameters:
    453: +  eigr - real part of eigenvalue
    454: -  eigi - imaginary part of eigenvalue

    456:    Notes:
    457:    If the eigenvalue is real, then `eigi` is set to zero. If PETSc is
    458:    configured with complex scalars the eigenvalue is stored
    459:    directly in `eigr` (`eigi` is set to zero).

    461:    The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
    462:    Eigenpairs are indexed according to the ordering criterion established
    463:    with `EPSSetWhichEigenpairs()`.

    465:    Level: beginner

    467: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`
    468: @*/
    469: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
    470: {
    471:   PetscInt  k,nconv;
    472: #if !defined(PETSC_USE_COMPLEX)
    473:   PetscInt  k2, iquad;
    474: #endif

    476:   PetscFunctionBegin;
    478:   EPSCheckSolved(eps,1);
    479:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    480:   PetscCall(EPS_GetActualConverged(eps,&nconv));
    481:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
    482:   if (nconv==eps->nconv) {
    483:     k = eps->perm[i];
    484: #if defined(PETSC_USE_COMPLEX)
    485:     if (eigr) *eigr = eps->eigr[k];
    486:     if (eigi) *eigi = 0;
    487: #else
    488:     if (eigr) *eigr = eps->eigr[k];
    489:     if (eigi) *eigi = eps->eigi[k];
    490: #endif
    491:   } else {
    492:     PetscCheck(eps->problem_type==EPS_BSE || eps->problem_type==EPS_HAMILT || eps->problem_type==EPS_LREP,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE, Hamiltonian, or LREP");
    493:     if (eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP) {
    494:       /* BSE problem, even index is +lambda, odd index is -lambda */
    495:       k = eps->perm[i/2];
    496: #if defined(PETSC_USE_COMPLEX)
    497:       if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
    498:       if (eigi) *eigi = 0;
    499: #else
    500:       if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
    501:       if (eigi) *eigi = eps->eigi[k];
    502: #endif
    503:     } else if (eps->problem_type==EPS_HAMILT) {
    504:       /* Hamiltonian eigenproblem */
    505:       k = eps->perm[i/2];
    506: #if defined(PETSC_USE_COMPLEX)
    507:       if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
    508:       if (eigi) *eigi = 0;
    509: #else
    510:       if (eps->eigi[k]==0.0) { /* real eigenvalue */
    511:         if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
    512:         if (eigi) *eigi = 0.0;
    513:       } else if (eps->eigr[k]==0.0) { /* purely imaginary eigenvalue */
    514:         if (eigr) *eigr = 0.0;
    515:         if (eigi) *eigi = (i%2)? -eps->eigi[k]: eps->eigi[k];
    516:       } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
    517:         iquad = i%2;  /* index within the 4 values */
    518:         if (i>1) {
    519:           k2 = eps->perm[(i-2)/2];
    520:           if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
    521:         }
    522:         if (eigr) *eigr = (iquad<2)? -eps->eigr[k]: eps->eigr[k];
    523:         if (eigi) *eigi = (iquad%3)? -eps->eigi[k]: eps->eigi[k];
    524:       }
    525: #endif
    526:     }
    527:   }
    528:   PetscFunctionReturn(PETSC_SUCCESS);
    529: }

    531: /*@
    532:    EPSGetEigenvector - Gets the `i`-th right eigenvector as computed by `EPSSolve()`.

    534:    Collective

    536:    Input Parameters:
    537: +  eps - the linear eigensolver context
    538: -  i   - index of the solution

    540:    Output Parameters:
    541: +  Vr   - real part of eigenvector
    542: -  Vi   - imaginary part of eigenvector

    544:    Notes:
    545:    The caller must provide valid `Vec` objects, i.e., they must be created
    546:    by the calling program with e.g. `MatCreateVecs()`.

    548:    If the corresponding eigenvalue is real, then `Vi` is set to zero. If PETSc is
    549:    configured with complex scalars the eigenvector is stored
    550:    directly in `Vr` (`Vi` is set to zero). In any case, the user can pass `NULL` in `Vr`
    551:    or `Vi` if one of them is not required.

    553:    The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
    554:    Eigenpairs are indexed according to the ordering criterion established
    555:    with `EPSSetWhichEigenpairs()`.

    557:    The 2-norm of the eigenvector is one unless the problem is generalized
    558:    Hermitian. In this case the eigenvector is normalized with respect to the
    559:    norm defined by the $B$ matrix.

    561:    In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
    562:    [](#sec:structured-vectors).

    564:    Level: beginner

    566: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`, `EPSGetLeftEigenvector()`
    567: @*/
    568: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
    569: {
    570:   PetscInt nconv;

    572:   PetscFunctionBegin;
    577:   EPSCheckSolved(eps,1);
    578:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    579:   PetscCall(EPS_GetActualConverged(eps,&nconv));
    580:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
    581:   PetscCall(EPSComputeVectors(eps));
    582:   PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
    583:   PetscFunctionReturn(PETSC_SUCCESS);
    584: }

    586: /*@
    587:    EPSGetLeftEigenvector - Gets the `i`-th left eigenvector as computed by `EPSSolve()`.

    589:    Collective

    591:    Input Parameters:
    592: +  eps - the linear eigensolver context
    593: -  i   - index of the solution

    595:    Output Parameters:
    596: +  Wr   - real part of left eigenvector
    597: -  Wi   - imaginary part of left eigenvector

    599:    Notes:
    600:    The caller must provide valid `Vec` objects, i.e., they must be created
    601:    by the calling program with e.g. `MatCreateVecs()`.

    603:    If the corresponding eigenvalue is real, then `Wi` is set to zero. If PETSc is
    604:    configured with complex scalars the eigenvector is stored directly in `Wr`
    605:    (`Wi` is set to zero). In any case, the user can pass `NULL` in `Wr` or `Wi` if
    606:    one of them is not required.

    608:    The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
    609:    Eigensolutions are indexed according to the ordering criterion established
    610:    with `EPSSetWhichEigenpairs()`.

    612:    Left eigenvectors are available only if the `twosided` flag was set, see
    613:    `EPSSetTwoSided()`.

    615:    In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
    616:    [](#sec:structured-vectors).

    618:    Level: intermediate

    620: .seealso: [](ch:eps), `EPSGetEigenvector()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSSetTwoSided()`
    621: @*/
    622: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
    623: {
    624:   PetscInt    nconv;
    625:   PetscBool   trivial,lrepreduced=PETSC_FALSE;
    626:   Mat         H;
    627:   IS          is[2];
    628:   Vec         v0,v1,z,z0,z1;

    630:   PetscFunctionBegin;
    635:   EPSCheckSolved(eps,1);
    636:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    637:   PetscCall(EPS_GetActualConverged(eps,&nconv));
    638:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");

    640:   trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP)? PETSC_TRUE: PETSC_FALSE;
    641:   if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");

    643:   PetscCall(EPSComputeVectors(eps));
    644:   if (trivial) {
    645:     PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
    646:     PetscCall(STGetMatrix(eps->st,0,&H));
    647:     if (eps->problem_type==EPS_LREP) PetscCall(SlepcCheckMatLREPReduced(H,&lrepreduced));
    648:     if (eps->problem_type==EPS_BSE || !lrepreduced) { /* change sign of bottom part of the vector */
    649:       PetscCall(MatNestGetISs(H,is,NULL));
    650:       if (Wr) {
    651:         PetscCall(VecGetSubVector(Wr,is[1],&v1));
    652:         PetscCall(VecScale(v1,-1.0));
    653:         PetscCall(VecRestoreSubVector(Wr,is[1],&v1));
    654:       }
    655: #if !defined(PETSC_USE_COMPLEX)
    656:       if (Wi) {
    657:         PetscCall(VecGetSubVector(Wi,is[1],&v1));
    658:         PetscCall(VecScale(v1,-1.0));
    659:         PetscCall(VecRestoreSubVector(Wi,is[1],&v1));
    660:       }
    661: #endif
    662:     } else if (eps->problem_type==EPS_LREP) {   /* swap the two parts of the vector */
    663:       if (Wr) {
    664:         PetscCall(VecDuplicate(Wr,&z));
    665:         PetscCall(VecCopy(Wr,z));
    666:         PetscCall(MatNestGetISs(H,is,NULL));
    667:         PetscCall(VecGetSubVector(Wr,is[0],&v0));
    668:         PetscCall(VecGetSubVector(Wr,is[1],&v1));
    669:         PetscCall(VecGetSubVector(z,is[0],&z0));
    670:         PetscCall(VecGetSubVector(z,is[1],&z1));
    671:         PetscCall(VecCopy(z0,v1));
    672:         PetscCall(VecCopy(z1,v0));
    673:         PetscCall(VecRestoreSubVector(Wr,is[0],&v0));
    674:         PetscCall(VecRestoreSubVector(Wr,is[1],&v1));
    675:         PetscCall(VecRestoreSubVector(z,is[0],&z0));
    676:         PetscCall(VecRestoreSubVector(z,is[1],&z1));
    677:         PetscCall(VecDestroy(&z));
    678:       }
    679:     }
    680:   } else {
    681:     PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
    682:   }
    683:   PetscFunctionReturn(PETSC_SUCCESS);
    684: }

    686: /*@
    687:    EPSGetErrorEstimate - Returns the error estimate associated to the `i`-th
    688:    computed eigenpair.

    690:    Not Collective

    692:    Input Parameters:
    693: +  eps - the linear eigensolver context
    694: -  i   - index of eigenpair

    696:    Output Parameter:
    697: .  errest - the error estimate

    699:    Note:
    700:    This is the error estimate used internally by the eigensolver. The actual
    701:    error bound can be computed with `EPSComputeError()`. See discussion at
    702:    section [](#sec:errbnd).

    704:    Level: advanced

    706: .seealso: [](ch:eps), [](#sec:errbnd), `EPSComputeError()`
    707: @*/
    708: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
    709: {
    710:   PetscInt nconv;

    712:   PetscFunctionBegin;
    714:   PetscAssertPointer(errest,3);
    715:   EPSCheckSolved(eps,1);
    716:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    717:   PetscCall(EPS_GetActualConverged(eps,&nconv));
    718:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
    719:   if (nconv==eps->nconv) {
    720:     *errest = eps->errest[eps->perm[i]];
    721:   } else {
    722:     PetscCheck(eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP || eps->problem_type==EPS_HAMILT,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Wrong problem type");
    723:     /* even index is +lambda, odd index is -lambda, assume both have same error */
    724:     *errest = eps->errest[eps->perm[i/2]];
    725:   }
    726:   PetscFunctionReturn(PETSC_SUCCESS);
    727: }

    729: /*
    730:    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
    731:    associated with an eigenpair.

    733:    Input Parameters:
    734:      trans - whether A' must be used instead of A
    735:      kr,ki - eigenvalue
    736:      xr,xi - eigenvector
    737:      z     - three work vectors (the second one not referenced in complex scalars)
    738: */
    739: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
    740: {
    741:   PetscInt       nmat;
    742:   Mat            A,B;
    743:   Vec            u,w;
    744:   PetscScalar    alpha;
    745: #if !defined(PETSC_USE_COMPLEX)
    746:   Vec            v;
    747:   PetscReal      ni,nr;
    748: #endif
    749:   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;

    751:   PetscFunctionBegin;
    752:   u = z[0]; w = z[2];
    753:   PetscCall(STGetNumMatrices(eps->st,&nmat));
    754:   PetscCall(STGetMatrix(eps->st,0,&A));
    755:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));

    757: #if !defined(PETSC_USE_COMPLEX)
    758:   v = z[1];
    759:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
    760: #endif
    761:     PetscCall((*matmult)(A,xr,u));                          /* u=A*x */
    762:     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
    763:       if (nmat>1) PetscCall((*matmult)(B,xr,w));
    764:       else PetscCall(VecCopy(xr,w));                        /* w=B*x */
    765:       alpha = trans? -PetscConj(kr): -kr;
    766:       PetscCall(VecAXPY(u,alpha,w));                        /* u=A*x-k*B*x */
    767:     }
    768:     PetscCall(VecNorm(u,NORM_2,norm));
    769: #if !defined(PETSC_USE_COMPLEX)
    770:   } else {
    771:     PetscCall((*matmult)(A,xr,u));                          /* u=A*xr */
    772:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
    773:       if (nmat>1) PetscCall((*matmult)(B,xr,v));
    774:       else PetscCall(VecCopy(xr,v));                        /* v=B*xr */
    775:       PetscCall(VecAXPY(u,-kr,v));                          /* u=A*xr-kr*B*xr */
    776:       if (nmat>1) PetscCall((*matmult)(B,xi,w));
    777:       else PetscCall(VecCopy(xi,w));                        /* w=B*xi */
    778:       PetscCall(VecAXPY(u,trans?-ki:ki,w));                 /* u=A*xr-kr*B*xr+ki*B*xi */
    779:     }
    780:     PetscCall(VecNorm(u,NORM_2,&nr));
    781:     PetscCall((*matmult)(A,xi,u));                          /* u=A*xi */
    782:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
    783:       PetscCall(VecAXPY(u,-kr,w));                          /* u=A*xi-kr*B*xi */
    784:       PetscCall(VecAXPY(u,trans?ki:-ki,v));                 /* u=A*xi-kr*B*xi-ki*B*xr */
    785:     }
    786:     PetscCall(VecNorm(u,NORM_2,&ni));
    787:     *norm = SlepcAbsEigenvalue(nr,ni);
    788:   }
    789: #endif
    790:   PetscFunctionReturn(PETSC_SUCCESS);
    791: }

    793: /*@
    794:    EPSComputeError - Computes the error (based on the residual norm) associated
    795:    with the `i`-th computed eigenpair.

    797:    Collective

    799:    Input Parameters:
    800: +  eps  - the linear eigensolver context
    801: .  i    - the solution index
    802: -  type - the type of error to compute, see `EPSErrorType`

    804:    Output Parameter:
    805: .  error - the error

    807:    Notes:
    808:    The error can be computed in various ways, all of them based on the residual
    809:    norm $\|Ax-\lambda Bx\|_2$ where $(\lambda,x)$ is the approximate eigenpair.

    811:    If the computation of left eigenvectors was enabled with `EPSSetTwoSided()`,
    812:    then the error will be computed using the maximum of the value above and
    813:    the left residual norm $\|y^*A-\lambda y^*B\|_2$, where $y$ is the approximate left
    814:    eigenvector.

    816:    Level: beginner

    818: .seealso: [](ch:eps), `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`, `EPSSetTwoSided()`
    819: @*/
    820: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
    821: {
    822:   Mat            A,B;
    823:   Vec            xr,xi,w[3];
    824:   PetscReal      t,vecnorm=1.0,errorl;
    825:   PetscScalar    kr,ki;
    826:   PetscBool      flg;

    828:   PetscFunctionBegin;
    832:   PetscAssertPointer(error,4);
    833:   EPSCheckSolved(eps,1);

    835:   /* allocate work vectors */
    836: #if defined(PETSC_USE_COMPLEX)
    837:   PetscCall(EPSSetWorkVecs(eps,3));
    838:   xi   = NULL;
    839:   w[1] = NULL;
    840: #else
    841:   PetscCall(EPSSetWorkVecs(eps,5));
    842:   xi   = eps->work[3];
    843:   w[1] = eps->work[4];
    844: #endif
    845:   xr   = eps->work[0];
    846:   w[0] = eps->work[1];
    847:   w[2] = eps->work[2];

    849:   /* compute residual norm */
    850:   PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
    851:   PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));

    853:   /* compute 2-norm of eigenvector */
    854:   if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));

    856:   /* if two-sided, compute left residual norm and take the maximum */
    857:   if (eps->twosided) {
    858:     PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
    859:     PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
    860:     *error = PetscMax(*error,errorl);
    861:   }

    863:   /* compute error */
    864:   switch (type) {
    865:     case EPS_ERROR_ABSOLUTE:
    866:       break;
    867:     case EPS_ERROR_RELATIVE:
    868:       *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
    869:       break;
    870:     case EPS_ERROR_BACKWARD:
    871:       /* initialization of matrix norms */
    872:       if (!eps->nrma) {
    873:         PetscCall(STGetMatrix(eps->st,0,&A));
    874:         PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
    875:         PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
    876:         PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
    877:       }
    878:       if (eps->isgeneralized) {
    879:         if (!eps->nrmb) {
    880:           PetscCall(STGetMatrix(eps->st,1,&B));
    881:           PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
    882:           PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
    883:           PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
    884:         }
    885:       } else eps->nrmb = 1.0;
    886:       t = SlepcAbsEigenvalue(kr,ki);
    887:       *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
    888:       break;
    889:     default:
    890:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
    891:   }
    892:   PetscFunctionReturn(PETSC_SUCCESS);
    893: }

    895: /*
    896:    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
    897:    for the recurrence that builds the right subspace.

    899:    Collective

    901:    Input Parameters:
    902: +  eps - the linear eigensolver context
    903: -  i   - iteration number

    905:    Output Parameter:
    906: .  breakdown - flag indicating that a breakdown has occurred

    908:    Notes:
    909:    The start vector is computed from another vector: for the first step (i=0),
    910:    the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
    911:    vector is created. Then this vector is forced to be in the range of OP (only
    912:    for generalized definite problems) and orthonormalized with respect to all
    913:    V-vectors up to i-1. The resulting vector is placed in V[i].

    915:    The flag breakdown is set to true if either i=0 and the vector belongs to the
    916:    deflation space, or i>0 and the vector is linearly dependent with respect
    917:    to the V-vectors.
    918: */
    919: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
    920: {
    921:   PetscReal      norm;
    922:   PetscBool      lindep;
    923:   Vec            w,z;

    925:   PetscFunctionBegin;

    929:   /* For the first step, use the first initial vector, otherwise a random one */
    930:   if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));

    932:   /* Force the vector to be in the range of OP for generalized problems with B-inner product */
    933:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
    934:     PetscCall(BVCreateVec(eps->V,&w));
    935:     PetscCall(BVCopyVec(eps->V,i,w));
    936:     PetscCall(BVGetColumn(eps->V,i,&z));
    937:     PetscCall(STApply(eps->st,w,z));
    938:     PetscCall(BVRestoreColumn(eps->V,i,&z));
    939:     PetscCall(VecDestroy(&w));
    940:   }

    942:   /* Orthonormalize the vector with respect to previous vectors */
    943:   PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
    944:   if (breakdown) *breakdown = lindep;
    945:   else if (lindep || norm == 0.0) {
    946:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
    947:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
    948:   }
    949:   PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
    950:   PetscFunctionReturn(PETSC_SUCCESS);
    951: }

    953: /*
    954:    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
    955:    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
    956: */
    957: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
    958: {
    959:   PetscReal      norm;
    960:   PetscBool      lindep;
    961:   Vec            w,z;

    963:   PetscFunctionBegin;

    967:   /* For the first step, use the first initial vector, otherwise a random one */
    968:   if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));

    970:   /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
    971:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
    972:     PetscCall(BVCreateVec(eps->W,&w));
    973:     PetscCall(BVCopyVec(eps->W,i,w));
    974:     PetscCall(BVGetColumn(eps->W,i,&z));
    975:     PetscCall(STApplyHermitianTranspose(eps->st,w,z));
    976:     PetscCall(BVRestoreColumn(eps->W,i,&z));
    977:     PetscCall(VecDestroy(&w));
    978:   }

    980:   /* Orthonormalize the vector with respect to previous vectors */
    981:   PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
    982:   if (breakdown) *breakdown = lindep;
    983:   else if (lindep || norm == 0.0) {
    984:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
    985:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
    986:   }
    987:   PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
    988:   PetscFunctionReturn(PETSC_SUCCESS);
    989: }