LCOV - code coverage report
Current view: top level - pep/interface - pepsolve.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 154 158 97.5 %
Date: 2024-03-29 00:34:52 Functions: 10 10 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    PEP routines related to the solution process
      12             : 
      13             :    References:
      14             : 
      15             :        [1] C. Campos and J.E. Roman, "Parallel iterative refinement in
      16             :            polynomial eigenvalue problems", Numer. Linear Algebra Appl.
      17             :            23(4):730-745, 2016.
      18             : */
      19             : 
      20             : #include <slepc/private/pepimpl.h>       /*I "slepcpep.h" I*/
      21             : #include <slepc/private/bvimpl.h>
      22             : #include <petscdraw.h>
      23             : 
      24             : static PetscBool  cited = PETSC_FALSE;
      25             : static const char citation[] =
      26             :   "@Article{slepc-pep-refine,\n"
      27             :   "   author = \"C. Campos and J. E. Roman\",\n"
      28             :   "   title = \"Parallel iterative refinement in polynomial eigenvalue problems\",\n"
      29             :   "   journal = \"Numer. Linear Algebra Appl.\",\n"
      30             :   "   volume = \"23\",\n"
      31             :   "   number = \"4\",\n"
      32             :   "   pages = \"730--745\",\n"
      33             :   "   year = \"2016,\"\n"
      34             :   "   doi = \"https://doi.org/10.1002/nla.2052\"\n"
      35             :   "}\n";
      36             : 
      37        1616 : PetscErrorCode PEPComputeVectors(PEP pep)
      38             : {
      39        1616 :   PetscFunctionBegin;
      40        1616 :   PEPCheckSolved(pep,1);
      41        1616 :   if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,computevectors);
      42        1616 :   pep->state = PEP_STATE_EIGENVECTORS;
      43        1616 :   PetscFunctionReturn(PETSC_SUCCESS);
      44             : }
      45             : 
      46         166 : PetscErrorCode PEPExtractVectors(PEP pep)
      47             : {
      48         166 :   PetscFunctionBegin;
      49         166 :   PEPCheckSolved(pep,1);
      50         166 :   if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,extractvectors);
      51         166 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54             : /*@
      55             :    PEPSolve - Solves the polynomial eigensystem.
      56             : 
      57             :    Collective
      58             : 
      59             :    Input Parameter:
      60             : .  pep - eigensolver context obtained from PEPCreate()
      61             : 
      62             :    Options Database Keys:
      63             : +  -pep_view - print information about the solver used
      64             : .  -pep_view_matk - view the coefficient matrix Ak (replace k by an integer from 0 to nmat-1)
      65             : .  -pep_view_vectors - view the computed eigenvectors
      66             : .  -pep_view_values - view the computed eigenvalues
      67             : .  -pep_converged_reason - print reason for convergence, and number of iterations
      68             : .  -pep_error_absolute - print absolute errors of each eigenpair
      69             : .  -pep_error_relative - print relative errors of each eigenpair
      70             : -  -pep_error_backward - print backward errors of each eigenpair
      71             : 
      72             :    Notes:
      73             :    All the command-line options listed above admit an optional argument specifying
      74             :    the viewer type and options. For instance, use '-pep_view_mat0 binary:amatrix.bin'
      75             :    to save the A matrix to a binary file, '-pep_view_values draw' to draw the computed
      76             :    eigenvalues graphically, or '-pep_error_relative :myerr.m:ascii_matlab' to save
      77             :    the errors in a file that can be executed in Matlab.
      78             : 
      79             :    Level: beginner
      80             : 
      81             : .seealso: PEPCreate(), PEPSetUp(), PEPDestroy(), PEPSetTolerances()
      82             : @*/
      83         193 : PetscErrorCode PEPSolve(PEP pep)
      84             : {
      85         193 :   PetscInt       i,k;
      86         193 :   PetscBool      flg,islinear;
      87         193 :   char           str[16];
      88             : 
      89         193 :   PetscFunctionBegin;
      90         193 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
      91         193 :   if (pep->state>=PEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
      92         193 :   PetscCall(PetscLogEventBegin(PEP_Solve,pep,0,0,0));
      93             : 
      94             :   /* call setup */
      95         193 :   PetscCall(PEPSetUp(pep));
      96         193 :   pep->nconv = 0;
      97         193 :   pep->its   = 0;
      98         193 :   k = pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1);
      99       22279 :   for (i=0;i<k;i++) {
     100       22086 :     pep->eigr[i]   = 0.0;
     101       22086 :     pep->eigi[i]   = 0.0;
     102       22086 :     pep->errest[i] = 0.0;
     103       22086 :     pep->perm[i]   = i;
     104             :   }
     105         193 :   PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view_pre"));
     106         193 :   PetscCall(RGViewFromOptions(pep->rg,NULL,"-rg_view"));
     107             : 
     108             :   /* Call solver */
     109         193 :   PetscUseTypeMethod(pep,solve);
     110         193 :   PetscCheck(pep->reason,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
     111         193 :   pep->state = PEP_STATE_SOLVED;
     112             : 
     113             :   /* Only the first nconv columns contain useful information */
     114         193 :   PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
     115             : 
     116         193 :   PetscCall(PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear));
     117         193 :   if (!islinear) {
     118         158 :     PetscCall(STPostSolve(pep->st));
     119             :     /* Map eigenvalues back to the original problem */
     120         158 :     PetscCall(STGetTransform(pep->st,&flg));
     121         158 :     if (flg) PetscTryTypeMethod(pep,backtransform);
     122             :   }
     123             : 
     124             : #if !defined(PETSC_USE_COMPLEX)
     125             :   /* reorder conjugate eigenvalues (positive imaginary first) */
     126             :   for (i=0;i<pep->nconv-1;i++) {
     127             :     if (pep->eigi[i] != 0) {
     128             :       if (pep->eigi[i] < 0) {
     129             :         pep->eigi[i] = -pep->eigi[i];
     130             :         pep->eigi[i+1] = -pep->eigi[i+1];
     131             :         /* the next correction only works with eigenvectors */
     132             :         PetscCall(PEPComputeVectors(pep));
     133             :         PetscCall(BVScaleColumn(pep->V,i+1,-1.0));
     134             :       }
     135             :       i++;
     136             :     }
     137             :   }
     138             : #endif
     139             : 
     140         193 :   if (pep->refine!=PEP_REFINE_NONE) PetscCall(PetscCitationsRegister(citation,&cited));
     141             : 
     142         193 :   if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
     143          15 :     PetscCall(PEPComputeVectors(pep));
     144          15 :     PetscCall(PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv));
     145             :   }
     146             : 
     147             :   /* sort eigenvalues according to pep->which parameter */
     148         193 :   PetscCall(SlepcSortEigenvalues(pep->sc,pep->nconv,pep->eigr,pep->eigi,pep->perm));
     149         193 :   PetscCall(PetscLogEventEnd(PEP_Solve,pep,0,0,0));
     150             : 
     151             :   /* various viewers */
     152         193 :   PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view"));
     153         193 :   PetscCall(PEPConvergedReasonViewFromOptions(pep));
     154         193 :   PetscCall(PEPErrorViewFromOptions(pep));
     155         193 :   PetscCall(PEPValuesViewFromOptions(pep));
     156         193 :   PetscCall(PEPVectorsViewFromOptions(pep));
     157         928 :   for (i=0;i<pep->nmat;i++) {
     158         735 :     PetscCall(PetscSNPrintf(str,sizeof(str),"-pep_view_mat%" PetscInt_FMT,i));
     159         735 :     PetscCall(MatViewFromOptions(pep->A[i],(PetscObject)pep,str));
     160             :   }
     161             : 
     162             :   /* Remove the initial subspace */
     163         193 :   pep->nini = 0;
     164         193 :   PetscFunctionReturn(PETSC_SUCCESS);
     165             : }
     166             : 
     167             : /*@
     168             :    PEPGetIterationNumber - Gets the current iteration number. If the
     169             :    call to PEPSolve() is complete, then it returns the number of iterations
     170             :    carried out by the solution method.
     171             : 
     172             :    Not Collective
     173             : 
     174             :    Input Parameter:
     175             : .  pep - the polynomial eigensolver context
     176             : 
     177             :    Output Parameter:
     178             : .  its - number of iterations
     179             : 
     180             :    Note:
     181             :    During the i-th iteration this call returns i-1. If PEPSolve() is
     182             :    complete, then parameter "its" contains either the iteration number at
     183             :    which convergence was successfully reached, or failure was detected.
     184             :    Call PEPGetConvergedReason() to determine if the solver converged or
     185             :    failed and why.
     186             : 
     187             :    Level: intermediate
     188             : 
     189             : .seealso: PEPGetConvergedReason(), PEPSetTolerances()
     190             : @*/
     191          29 : PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
     192             : {
     193          29 :   PetscFunctionBegin;
     194          29 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     195          29 :   PetscAssertPointer(its,2);
     196          29 :   *its = pep->its;
     197          29 :   PetscFunctionReturn(PETSC_SUCCESS);
     198             : }
     199             : 
     200             : /*@
     201             :    PEPGetConverged - Gets the number of converged eigenpairs.
     202             : 
     203             :    Not Collective
     204             : 
     205             :    Input Parameter:
     206             : .  pep - the polynomial eigensolver context
     207             : 
     208             :    Output Parameter:
     209             : .  nconv - number of converged eigenpairs
     210             : 
     211             :    Note:
     212             :    This function should be called after PEPSolve() has finished.
     213             : 
     214             :    Level: beginner
     215             : 
     216             : .seealso: PEPSetDimensions(), PEPSolve(), PEPGetEigenpair()
     217             : @*/
     218          89 : PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
     219             : {
     220          89 :   PetscFunctionBegin;
     221          89 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     222          89 :   PetscAssertPointer(nconv,2);
     223          89 :   PEPCheckSolved(pep,1);
     224          89 :   *nconv = pep->nconv;
     225          89 :   PetscFunctionReturn(PETSC_SUCCESS);
     226             : }
     227             : 
     228             : /*@
     229             :    PEPGetConvergedReason - Gets the reason why the PEPSolve() iteration was
     230             :    stopped.
     231             : 
     232             :    Not Collective
     233             : 
     234             :    Input Parameter:
     235             : .  pep - the polynomial eigensolver context
     236             : 
     237             :    Output Parameter:
     238             : .  reason - negative value indicates diverged, positive value converged
     239             : 
     240             :    Options Database Key:
     241             : .  -pep_converged_reason - print the reason to a viewer
     242             : 
     243             :    Notes:
     244             :    Possible values for reason are
     245             : +  PEP_CONVERGED_TOL - converged up to tolerance
     246             : .  PEP_CONVERGED_USER - converged due to a user-defined condition
     247             : .  PEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
     248             : .  PEP_DIVERGED_BREAKDOWN - generic breakdown in method
     249             : -  PEP_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
     250             : 
     251             :    Can only be called after the call to PEPSolve() is complete.
     252             : 
     253             :    Level: intermediate
     254             : 
     255             : .seealso: PEPSetTolerances(), PEPSolve(), PEPConvergedReason
     256             : @*/
     257          26 : PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
     258             : {
     259          26 :   PetscFunctionBegin;
     260          26 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     261          26 :   PetscAssertPointer(reason,2);
     262          26 :   PEPCheckSolved(pep,1);
     263          26 :   *reason = pep->reason;
     264          26 :   PetscFunctionReturn(PETSC_SUCCESS);
     265             : }
     266             : 
     267             : /*@C
     268             :    PEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
     269             :    PEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
     270             : 
     271             :    Collective
     272             : 
     273             :    Input Parameters:
     274             : +  pep - polynomial eigensolver context
     275             : -  i   - index of the solution
     276             : 
     277             :    Output Parameters:
     278             : +  eigr - real part of eigenvalue
     279             : .  eigi - imaginary part of eigenvalue
     280             : .  Vr   - real part of eigenvector
     281             : -  Vi   - imaginary part of eigenvector
     282             : 
     283             :    Notes:
     284             :    It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
     285             :    required. Otherwise, the caller must provide valid Vec objects, i.e.,
     286             :    they must be created by the calling program with e.g. MatCreateVecs().
     287             : 
     288             :    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
     289             :    configured with complex scalars the eigenvalue is stored
     290             :    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
     291             :    set to zero). In any case, the user can pass NULL in Vr or Vi if one of
     292             :    them is not required.
     293             : 
     294             :    The index i should be a value between 0 and nconv-1 (see PEPGetConverged()).
     295             :    Eigenpairs are indexed according to the ordering criterion established
     296             :    with PEPSetWhichEigenpairs().
     297             : 
     298             :    Level: beginner
     299             : 
     300             : .seealso: PEPSolve(), PEPGetConverged(), PEPSetWhichEigenpairs()
     301             : @*/
     302        1600 : PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
     303             : {
     304        1600 :   PetscInt       k;
     305             : 
     306        1600 :   PetscFunctionBegin;
     307        1600 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     308        6400 :   PetscValidLogicalCollectiveInt(pep,i,2);
     309        1600 :   if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,5); PetscCheckSameComm(pep,1,Vr,5); }
     310        1600 :   if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,6); PetscCheckSameComm(pep,1,Vi,6); }
     311        1600 :   PEPCheckSolved(pep,1);
     312        1600 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     313        1600 :   PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
     314             : 
     315        1600 :   PetscCall(PEPComputeVectors(pep));
     316        1600 :   k = pep->perm[i];
     317             : 
     318             :   /* eigenvalue */
     319             : #if defined(PETSC_USE_COMPLEX)
     320        1600 :   if (eigr) *eigr = pep->eigr[k];
     321        1600 :   if (eigi) *eigi = 0;
     322             : #else
     323             :   if (eigr) *eigr = pep->eigr[k];
     324             :   if (eigi) *eigi = pep->eigi[k];
     325             : #endif
     326             : 
     327             :   /* eigenvector */
     328        1600 :   PetscCall(BV_GetEigenvector(pep->V,k,pep->eigi[k],Vr,Vi));
     329        1600 :   PetscFunctionReturn(PETSC_SUCCESS);
     330             : }
     331             : 
     332             : /*@
     333             :    PEPGetErrorEstimate - Returns the error estimate associated to the i-th
     334             :    computed eigenpair.
     335             : 
     336             :    Not Collective
     337             : 
     338             :    Input Parameters:
     339             : +  pep - polynomial eigensolver context
     340             : -  i   - index of eigenpair
     341             : 
     342             :    Output Parameter:
     343             : .  errest - the error estimate
     344             : 
     345             :    Notes:
     346             :    This is the error estimate used internally by the eigensolver. The actual
     347             :    error bound can be computed with PEPComputeError(). See also the users
     348             :    manual for details.
     349             : 
     350             :    Level: advanced
     351             : 
     352             : .seealso: PEPComputeError()
     353             : @*/
     354           4 : PetscErrorCode PEPGetErrorEstimate(PEP pep,PetscInt i,PetscReal *errest)
     355             : {
     356           4 :   PetscFunctionBegin;
     357           4 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     358           4 :   PetscAssertPointer(errest,3);
     359           4 :   PEPCheckSolved(pep,1);
     360           4 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     361           4 :   PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
     362           4 :   *errest = pep->errest[pep->perm[i]];
     363           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     364             : }
     365             : 
     366             : /*
     367             :    PEPComputeResidualNorm_Private - Computes the norm of the residual vector
     368             :    associated with an eigenpair.
     369             : 
     370             :    Input Parameters:
     371             :      kr,ki - eigenvalue
     372             :      xr,xi - eigenvector
     373             :      z     - array of 4 work vectors (z[2],z[3] not referenced in complex scalars)
     374             : */
     375        1235 : PetscErrorCode PEPComputeResidualNorm_Private(PEP pep,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
     376             : {
     377        1235 :   Mat            *A=pep->A;
     378        1235 :   PetscInt       i,nmat=pep->nmat;
     379        1235 :   PetscScalar    t[20],*vals=t,*ivals=NULL;
     380        1235 :   Vec            u,w;
     381             : #if !defined(PETSC_USE_COMPLEX)
     382             :   Vec            ui,wi;
     383             :   PetscReal      ni;
     384             :   PetscBool      imag;
     385             :   PetscScalar    it[20];
     386             : #endif
     387             : 
     388        1235 :   PetscFunctionBegin;
     389        1235 :   u = z[0]; w = z[1];
     390        1235 :   PetscCall(VecSet(u,0.0));
     391             : #if !defined(PETSC_USE_COMPLEX)
     392             :   ui = z[2]; wi = z[3];
     393             :   ivals = it;
     394             : #endif
     395        1235 :   if (nmat>20) {
     396           0 :     PetscCall(PetscMalloc1(nmat,&vals));
     397             : #if !defined(PETSC_USE_COMPLEX)
     398             :     PetscCall(PetscMalloc1(nmat,&ivals));
     399             : #endif
     400             :   }
     401        1235 :   PetscCall(PEPEvaluateBasis(pep,kr,ki,vals,ivals));
     402             : #if !defined(PETSC_USE_COMPLEX)
     403             :   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
     404             :     imag = PETSC_FALSE;
     405             :   else {
     406             :     imag = PETSC_TRUE;
     407             :     PetscCall(VecSet(ui,0.0));
     408             :   }
     409             : #endif
     410        5349 :   for (i=0;i<nmat;i++) {
     411        4114 :     if (vals[i]!=0.0) {
     412        4114 :       PetscCall(MatMult(A[i],xr,w));
     413        4114 :       PetscCall(VecAXPY(u,vals[i],w));
     414             :     }
     415             : #if !defined(PETSC_USE_COMPLEX)
     416             :     if (imag) {
     417             :       if (ivals[i]!=0 || vals[i]!=0) {
     418             :         PetscCall(MatMult(A[i],xi,wi));
     419             :         if (vals[i]==0) PetscCall(MatMult(A[i],xr,w));
     420             :       }
     421             :       if (ivals[i]!=0) {
     422             :         PetscCall(VecAXPY(u,-ivals[i],wi));
     423             :         PetscCall(VecAXPY(ui,ivals[i],w));
     424             :       }
     425             :       if (vals[i]!=0) PetscCall(VecAXPY(ui,vals[i],wi));
     426             :     }
     427             : #endif
     428             :   }
     429        1235 :   PetscCall(VecNorm(u,NORM_2,norm));
     430             : #if !defined(PETSC_USE_COMPLEX)
     431             :   if (imag) {
     432             :     PetscCall(VecNorm(ui,NORM_2,&ni));
     433             :     *norm = SlepcAbsEigenvalue(*norm,ni);
     434             :   }
     435             : #endif
     436        1235 :   if (nmat>20) {
     437           0 :     PetscCall(PetscFree(vals));
     438             : #if !defined(PETSC_USE_COMPLEX)
     439             :     PetscCall(PetscFree(ivals));
     440             : #endif
     441             :   }
     442        1235 :   PetscFunctionReturn(PETSC_SUCCESS);
     443             : }
     444             : 
     445             : /*@
     446             :    PEPComputeError - Computes the error (based on the residual norm) associated
     447             :    with the i-th computed eigenpair.
     448             : 
     449             :    Collective
     450             : 
     451             :    Input Parameters:
     452             : +  pep  - the polynomial eigensolver context
     453             : .  i    - the solution index
     454             : -  type - the type of error to compute
     455             : 
     456             :    Output Parameter:
     457             : .  error - the error
     458             : 
     459             :    Notes:
     460             :    The error can be computed in various ways, all of them based on the residual
     461             :    norm ||P(l)x||_2 where l is the eigenvalue and x is the eigenvector.
     462             :    See the users guide for additional details.
     463             : 
     464             :    Level: beginner
     465             : 
     466             : .seealso: PEPErrorType, PEPSolve(), PEPGetErrorEstimate()
     467             : @*/
     468        1057 : PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
     469             : {
     470        1057 :   Vec            xr,xi,w[4];
     471        1057 :   PetscScalar    kr,ki;
     472        1057 :   PetscReal      t,z=0.0;
     473        1057 :   PetscInt       j;
     474        1057 :   PetscBool      flg;
     475             : 
     476        1057 :   PetscFunctionBegin;
     477        1057 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     478        4228 :   PetscValidLogicalCollectiveInt(pep,i,2);
     479        4228 :   PetscValidLogicalCollectiveEnum(pep,type,3);
     480        1057 :   PetscAssertPointer(error,4);
     481        1057 :   PEPCheckSolved(pep,1);
     482             : 
     483             :   /* allocate work vectors */
     484             : #if defined(PETSC_USE_COMPLEX)
     485        1057 :   PetscCall(PEPSetWorkVecs(pep,3));
     486        1057 :   xi   = NULL;
     487        1057 :   w[2] = NULL;
     488        1057 :   w[3] = NULL;
     489             : #else
     490             :   PetscCall(PEPSetWorkVecs(pep,6));
     491             :   xi   = pep->work[3];
     492             :   w[2] = pep->work[4];
     493             :   w[3] = pep->work[5];
     494             : #endif
     495        1057 :   xr   = pep->work[0];
     496        1057 :   w[0] = pep->work[1];
     497        1057 :   w[1] = pep->work[2];
     498             : 
     499             :   /* compute residual norms */
     500        1057 :   PetscCall(PEPGetEigenpair(pep,i,&kr,&ki,xr,xi));
     501        1057 :   PetscCall(PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error));
     502             : 
     503             :   /* compute error */
     504        1057 :   switch (type) {
     505             :     case PEP_ERROR_ABSOLUTE:
     506             :       break;
     507          26 :     case PEP_ERROR_RELATIVE:
     508          26 :       *error /= SlepcAbsEigenvalue(kr,ki);
     509          26 :       break;
     510        1019 :     case PEP_ERROR_BACKWARD:
     511             :       /* initialization of matrix norms */
     512        1019 :       if (!pep->nrma[pep->nmat-1]) {
     513         642 :         for (j=0;j<pep->nmat;j++) {
     514         489 :           PetscCall(MatHasOperation(pep->A[j],MATOP_NORM,&flg));
     515         489 :           PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
     516         489 :           PetscCall(MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]));
     517             :         }
     518             :       }
     519        1019 :       t = SlepcAbsEigenvalue(kr,ki);
     520        4304 :       for (j=pep->nmat-1;j>=0;j--) {
     521        3285 :         z = z*t+pep->nrma[j];
     522             :       }
     523        1019 :       *error /= z;
     524        1019 :       break;
     525           0 :     default:
     526           0 :       SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
     527             :   }
     528        1057 :   PetscFunctionReturn(PETSC_SUCCESS);
     529             : }

Generated by: LCOV version 1.14