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.

    123:    Level: beginner

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    235:    Not Collective

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

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

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

    250:    Level: intermediate

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

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

    266:    Not Collective

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

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

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

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

    280:    Level: beginner

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

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

    298:    Not Collective

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

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

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

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

    314:    Level: intermediate

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

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

    332:    Collective

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

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

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

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

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

    353:    Level: intermediate

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

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

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

    389:    Collective

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

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

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

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

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

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

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

    422:    Level: beginner

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

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

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

    445:    Not Collective

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

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

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

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

    464:    Level: beginner

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

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

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

    533:    Collective

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

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

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

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

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

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

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

    563:    Level: beginner

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

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

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

    588:    Collective

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

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

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

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

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

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

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

    617:    Level: intermediate

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

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

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

    642:   PetscCall(EPSComputeVectors(eps));
    643:   if (trivial) {
    644:     PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
    645:     if (eps->problem_type==EPS_BSE) {   /* change sign of bottom part of the vector */
    646:       PetscCall(STGetMatrix(eps->st,0,&H));
    647:       PetscCall(MatNestGetISs(H,is,NULL));
    648:       if (Wr) {
    649:         PetscCall(VecGetSubVector(Wr,is[1],&v));
    650:         PetscCall(VecScale(v,-1.0));
    651:         PetscCall(VecRestoreSubVector(Wr,is[1],&v));
    652:       }
    653: #if !defined(PETSC_USE_COMPLEX)
    654:       if (Wi) {
    655:         PetscCall(VecGetSubVector(Wi,is[1],&v));
    656:         PetscCall(VecScale(v,-1.0));
    657:         PetscCall(VecRestoreSubVector(Wi,is[1],&v));
    658:       }
    659: #endif
    660:     }
    661:   } else {
    662:     PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
    663:   }
    664:   PetscFunctionReturn(PETSC_SUCCESS);
    665: }

    667: /*@
    668:    EPSGetErrorEstimate - Returns the error estimate associated to the `i`-th
    669:    computed eigenpair.

    671:    Not Collective

    673:    Input Parameters:
    674: +  eps - the linear eigensolver context
    675: -  i   - index of eigenpair

    677:    Output Parameter:
    678: .  errest - the error estimate

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

    685:    Level: advanced

    687: .seealso: [](ch:eps), [](#sec:errbnd), `EPSComputeError()`
    688: @*/
    689: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
    690: {
    691:   PetscInt nconv;

    693:   PetscFunctionBegin;
    695:   PetscAssertPointer(errest,3);
    696:   EPSCheckSolved(eps,1);
    697:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    698:   PetscCall(EPS_GetActualConverged(eps,&nconv));
    699:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
    700:   if (nconv==eps->nconv) {
    701:     *errest = eps->errest[eps->perm[i]];
    702:   } else {
    703:     PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
    704:     /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
    705:     *errest = eps->errest[eps->perm[i/2]];
    706:   }
    707:   PetscFunctionReturn(PETSC_SUCCESS);
    708: }

    710: /*
    711:    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
    712:    associated with an eigenpair.

    714:    Input Parameters:
    715:      trans - whether A' must be used instead of A
    716:      kr,ki - eigenvalue
    717:      xr,xi - eigenvector
    718:      z     - three work vectors (the second one not referenced in complex scalars)
    719: */
    720: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
    721: {
    722:   PetscInt       nmat;
    723:   Mat            A,B;
    724:   Vec            u,w;
    725:   PetscScalar    alpha;
    726: #if !defined(PETSC_USE_COMPLEX)
    727:   Vec            v;
    728:   PetscReal      ni,nr;
    729: #endif
    730:   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;

    732:   PetscFunctionBegin;
    733:   u = z[0]; w = z[2];
    734:   PetscCall(STGetNumMatrices(eps->st,&nmat));
    735:   PetscCall(STGetMatrix(eps->st,0,&A));
    736:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));

    738: #if !defined(PETSC_USE_COMPLEX)
    739:   v = z[1];
    740:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
    741: #endif
    742:     PetscCall((*matmult)(A,xr,u));                          /* u=A*x */
    743:     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
    744:       if (nmat>1) PetscCall((*matmult)(B,xr,w));
    745:       else PetscCall(VecCopy(xr,w));                        /* w=B*x */
    746:       alpha = trans? -PetscConj(kr): -kr;
    747:       PetscCall(VecAXPY(u,alpha,w));                        /* u=A*x-k*B*x */
    748:     }
    749:     PetscCall(VecNorm(u,NORM_2,norm));
    750: #if !defined(PETSC_USE_COMPLEX)
    751:   } else {
    752:     PetscCall((*matmult)(A,xr,u));                          /* u=A*xr */
    753:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
    754:       if (nmat>1) PetscCall((*matmult)(B,xr,v));
    755:       else PetscCall(VecCopy(xr,v));                        /* v=B*xr */
    756:       PetscCall(VecAXPY(u,-kr,v));                          /* u=A*xr-kr*B*xr */
    757:       if (nmat>1) PetscCall((*matmult)(B,xi,w));
    758:       else PetscCall(VecCopy(xi,w));                        /* w=B*xi */
    759:       PetscCall(VecAXPY(u,trans?-ki:ki,w));                 /* u=A*xr-kr*B*xr+ki*B*xi */
    760:     }
    761:     PetscCall(VecNorm(u,NORM_2,&nr));
    762:     PetscCall((*matmult)(A,xi,u));                          /* u=A*xi */
    763:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
    764:       PetscCall(VecAXPY(u,-kr,w));                          /* u=A*xi-kr*B*xi */
    765:       PetscCall(VecAXPY(u,trans?ki:-ki,v));                 /* u=A*xi-kr*B*xi-ki*B*xr */
    766:     }
    767:     PetscCall(VecNorm(u,NORM_2,&ni));
    768:     *norm = SlepcAbsEigenvalue(nr,ni);
    769:   }
    770: #endif
    771:   PetscFunctionReturn(PETSC_SUCCESS);
    772: }

    774: /*@
    775:    EPSComputeError - Computes the error (based on the residual norm) associated
    776:    with the `i`-th computed eigenpair.

    778:    Collective

    780:    Input Parameters:
    781: +  eps  - the linear eigensolver context
    782: .  i    - the solution index
    783: -  type - the type of error to compute, see `EPSErrorType`

    785:    Output Parameter:
    786: .  error - the error

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

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

    797:    Level: beginner

    799: .seealso: [](ch:eps), `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`, `EPSSetTwoSided()`
    800: @*/
    801: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
    802: {
    803:   Mat            A,B;
    804:   Vec            xr,xi,w[3];
    805:   PetscReal      t,vecnorm=1.0,errorl;
    806:   PetscScalar    kr,ki;
    807:   PetscBool      flg;

    809:   PetscFunctionBegin;
    813:   PetscAssertPointer(error,4);
    814:   EPSCheckSolved(eps,1);

    816:   /* allocate work vectors */
    817: #if defined(PETSC_USE_COMPLEX)
    818:   PetscCall(EPSSetWorkVecs(eps,3));
    819:   xi   = NULL;
    820:   w[1] = NULL;
    821: #else
    822:   PetscCall(EPSSetWorkVecs(eps,5));
    823:   xi   = eps->work[3];
    824:   w[1] = eps->work[4];
    825: #endif
    826:   xr   = eps->work[0];
    827:   w[0] = eps->work[1];
    828:   w[2] = eps->work[2];

    830:   /* compute residual norm */
    831:   PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
    832:   PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));

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

    837:   /* if two-sided, compute left residual norm and take the maximum */
    838:   if (eps->twosided) {
    839:     PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
    840:     PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
    841:     *error = PetscMax(*error,errorl);
    842:   }

    844:   /* compute error */
    845:   switch (type) {
    846:     case EPS_ERROR_ABSOLUTE:
    847:       break;
    848:     case EPS_ERROR_RELATIVE:
    849:       *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
    850:       break;
    851:     case EPS_ERROR_BACKWARD:
    852:       /* initialization of matrix norms */
    853:       if (!eps->nrma) {
    854:         PetscCall(STGetMatrix(eps->st,0,&A));
    855:         PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
    856:         PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
    857:         PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
    858:       }
    859:       if (eps->isgeneralized) {
    860:         if (!eps->nrmb) {
    861:           PetscCall(STGetMatrix(eps->st,1,&B));
    862:           PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
    863:           PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
    864:           PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
    865:         }
    866:       } else eps->nrmb = 1.0;
    867:       t = SlepcAbsEigenvalue(kr,ki);
    868:       *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
    869:       break;
    870:     default:
    871:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
    872:   }
    873:   PetscFunctionReturn(PETSC_SUCCESS);
    874: }

    876: /*
    877:    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
    878:    for the recurrence that builds the right subspace.

    880:    Collective

    882:    Input Parameters:
    883: +  eps - the linear eigensolver context
    884: -  i   - iteration number

    886:    Output Parameter:
    887: .  breakdown - flag indicating that a breakdown has occurred

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

    896:    The flag breakdown is set to true if either i=0 and the vector belongs to the
    897:    deflation space, or i>0 and the vector is linearly dependent with respect
    898:    to the V-vectors.
    899: */
    900: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
    901: {
    902:   PetscReal      norm;
    903:   PetscBool      lindep;
    904:   Vec            w,z;

    906:   PetscFunctionBegin;

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

    913:   /* Force the vector to be in the range of OP for generalized problems with B-inner product */
    914:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
    915:     PetscCall(BVCreateVec(eps->V,&w));
    916:     PetscCall(BVCopyVec(eps->V,i,w));
    917:     PetscCall(BVGetColumn(eps->V,i,&z));
    918:     PetscCall(STApply(eps->st,w,z));
    919:     PetscCall(BVRestoreColumn(eps->V,i,&z));
    920:     PetscCall(VecDestroy(&w));
    921:   }

    923:   /* Orthonormalize the vector with respect to previous vectors */
    924:   PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
    925:   if (breakdown) *breakdown = lindep;
    926:   else if (lindep || norm == 0.0) {
    927:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
    928:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
    929:   }
    930:   PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
    931:   PetscFunctionReturn(PETSC_SUCCESS);
    932: }

    934: /*
    935:    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
    936:    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
    937: */
    938: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
    939: {
    940:   PetscReal      norm;
    941:   PetscBool      lindep;
    942:   Vec            w,z;

    944:   PetscFunctionBegin;

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

    951:   /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
    952:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
    953:     PetscCall(BVCreateVec(eps->W,&w));
    954:     PetscCall(BVCopyVec(eps->W,i,w));
    955:     PetscCall(BVGetColumn(eps->W,i,&z));
    956:     PetscCall(STApplyHermitianTranspose(eps->st,w,z));
    957:     PetscCall(BVRestoreColumn(eps->W,i,&z));
    958:     PetscCall(VecDestroy(&w));
    959:   }

    961:   /* Orthonormalize the vector with respect to previous vectors */
    962:   PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
    963:   if (breakdown) *breakdown = lindep;
    964:   else if (lindep || norm == 0.0) {
    965:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
    966:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
    967:   }
    968:   PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
    969:   PetscFunctionReturn(PETSC_SUCCESS);
    970: }