LCOV - code coverage report
Current view: top level - nep/interface - nepsolve.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 287 291 98.6 %
Date: 2021-08-02 00:35:43 Functions: 16 16 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-2021, 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             :    NEP routines related to the solution process
      12             : 
      13             :    References:
      14             : 
      15             :        [1] C. Campos and J.E. Roman, "NEP: a module for the parallel solution
      16             :            of nonlinear eigenvalue problems in SLEPc", ACM Trans. Math. Soft.
      17             :            (to appear), arXiv:1910.11712, 2021.
      18             : */
      19             : 
      20             : #include <slepc/private/nepimpl.h>       /*I "slepcnep.h" I*/
      21             : #include <slepc/private/bvimpl.h>        /*I "slepcbv.h" I*/
      22             : #include <petscdraw.h>
      23             : 
      24             : static PetscBool  cited = PETSC_FALSE;
      25             : static const char citation[] =
      26             :   "@Article{slepc-nep,\n"
      27             :   "   author = \"C. Campos and J. E. Roman\",\n"
      28             :   "   title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
      29             :   "   journal = \"{ACM} Trans. Math. Software\",\n"
      30             :   "   volume = \"47\",\n"
      31             :   "   number = \"3\",\n"
      32             :   "   pages = \"23:1--23:29\",\n"
      33             :   "   year = \"2021\",\n"
      34             :   "   doi = \"10.1145/3447544\"\n"
      35             :   "}\n";
      36             : 
      37         694 : PetscErrorCode NEPComputeVectors(NEP nep)
      38             : {
      39         694 :   PetscErrorCode ierr;
      40             : 
      41         694 :   PetscFunctionBegin;
      42         694 :   NEPCheckSolved(nep,1);
      43         694 :   if (nep->state==NEP_STATE_SOLVED && nep->ops->computevectors) {
      44          67 :     ierr = (*nep->ops->computevectors)(nep);CHKERRQ(ierr);
      45             :   }
      46         694 :   nep->state = NEP_STATE_EIGENVECTORS;
      47         694 :   PetscFunctionReturn(0);
      48             : }
      49             : 
      50             : /*@
      51             :    NEPSolve - Solves the nonlinear eigensystem.
      52             : 
      53             :    Collective on nep
      54             : 
      55             :    Input Parameter:
      56             : .  nep - eigensolver context obtained from NEPCreate()
      57             : 
      58             :    Options Database Keys:
      59             : +  -nep_view - print information about the solver used
      60             : .  -nep_view_vectors binary - save the computed eigenvectors to the default binary viewer
      61             : .  -nep_view_values - print computed eigenvalues
      62             : .  -nep_converged_reason - print reason for convergence, and number of iterations
      63             : .  -nep_error_absolute - print absolute errors of each eigenpair
      64             : -  -nep_error_relative - print relative errors of each eigenpair
      65             : 
      66             :    Level: beginner
      67             : 
      68             : .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
      69             : @*/
      70         123 : PetscErrorCode NEPSolve(NEP nep)
      71             : {
      72         123 :   PetscErrorCode ierr;
      73         123 :   PetscInt       i;
      74             : 
      75         123 :   PetscFunctionBegin;
      76         123 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
      77         123 :   if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(0);
      78         123 :   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
      79         123 :   ierr = PetscLogEventBegin(NEP_Solve,nep,0,0,0);CHKERRQ(ierr);
      80             : 
      81             :   /* call setup */
      82         123 :   ierr = NEPSetUp(nep);CHKERRQ(ierr);
      83         123 :   nep->nconv = 0;
      84         123 :   nep->its = 0;
      85        5884 :   for (i=0;i<nep->ncv;i++) {
      86        5761 :     nep->eigr[i]   = 0.0;
      87        5761 :     nep->eigi[i]   = 0.0;
      88        5761 :     nep->errest[i] = 0.0;
      89        5761 :     nep->perm[i]   = i;
      90             :   }
      91         123 :   ierr = NEPViewFromOptions(nep,NULL,"-nep_view_pre");CHKERRQ(ierr);
      92         123 :   ierr = RGViewFromOptions(nep->rg,NULL,"-rg_view");CHKERRQ(ierr);
      93             : 
      94             :   /* call solver */
      95         123 :   ierr = (*nep->ops->solve)(nep);CHKERRQ(ierr);
      96         123 :   if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
      97         123 :   nep->state = NEP_STATE_SOLVED;
      98             : 
      99             :   /* Only the first nconv columns contain useful information */
     100         123 :   ierr = BVSetActiveColumns(nep->V,0,nep->nconv);CHKERRQ(ierr);
     101         123 :   if (nep->twosided) { ierr = BVSetActiveColumns(nep->W,0,nep->nconv);CHKERRQ(ierr); }
     102             : 
     103         123 :   if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
     104          13 :     ierr = NEPComputeVectors(nep);CHKERRQ(ierr);
     105          13 :     ierr = NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);CHKERRQ(ierr);
     106          13 :     nep->state = NEP_STATE_EIGENVECTORS;
     107             :   }
     108             : 
     109             :   /* sort eigenvalues according to nep->which parameter */
     110         123 :   ierr = SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);CHKERRQ(ierr);
     111         123 :   ierr = PetscLogEventEnd(NEP_Solve,nep,0,0,0);CHKERRQ(ierr);
     112             : 
     113             :   /* various viewers */
     114         123 :   ierr = NEPViewFromOptions(nep,NULL,"-nep_view");CHKERRQ(ierr);
     115         123 :   ierr = NEPConvergedReasonViewFromOptions(nep);CHKERRQ(ierr);
     116         123 :   ierr = NEPErrorViewFromOptions(nep);CHKERRQ(ierr);
     117         123 :   ierr = NEPValuesViewFromOptions(nep);CHKERRQ(ierr);
     118         123 :   ierr = NEPVectorsViewFromOptions(nep);CHKERRQ(ierr);
     119             : 
     120             :   /* Remove the initial subspace */
     121         123 :   nep->nini = 0;
     122             : 
     123             :   /* Reset resolvent information */
     124         123 :   ierr = MatDestroy(&nep->resolvent);CHKERRQ(ierr);
     125         123 :   PetscFunctionReturn(0);
     126             : }
     127             : 
     128             : /*@
     129             :    NEPProjectOperator - Computes the projection of the nonlinear operator.
     130             : 
     131             :    Collective on nep
     132             : 
     133             :    Input Parameters:
     134             : +  nep - the nonlinear eigensolver context
     135             : .  j0  - initial index
     136             : -  j1  - final index
     137             : 
     138             :    Notes:
     139             :    This is available for split operator only.
     140             : 
     141             :    The nonlinear operator T(lambda) is projected onto span(V), where V is
     142             :    an orthonormal basis built internally by the solver. The projected
     143             :    operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
     144             :    computes all matrices Ei = V'*A_i*V, and stores them in the extra
     145             :    matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
     146             :    the previous ones are assumed to be available already.
     147             : 
     148             :    Level: developer
     149             : 
     150             : .seealso: NEPSetSplitOperator()
     151             : @*/
     152           1 : PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
     153             : {
     154           1 :   PetscErrorCode ierr;
     155           1 :   PetscInt       k;
     156           1 :   Mat            G;
     157             : 
     158           1 :   PetscFunctionBegin;
     159           1 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     160           2 :   PetscValidLogicalCollectiveInt(nep,j0,2);
     161           2 :   PetscValidLogicalCollectiveInt(nep,j1,3);
     162           1 :   NEPCheckProblem(nep,1);
     163           1 :   NEPCheckSplit(nep,1);
     164           1 :   ierr = BVSetActiveColumns(nep->V,j0,j1);CHKERRQ(ierr);
     165           4 :   for (k=0;k<nep->nt;k++) {
     166           3 :     ierr = DSGetMat(nep->ds,DSMatExtra[k],&G);CHKERRQ(ierr);
     167           3 :     ierr = BVMatProject(nep->V,nep->A[k],nep->V,G);CHKERRQ(ierr);
     168           3 :     ierr = DSRestoreMat(nep->ds,DSMatExtra[k],&G);CHKERRQ(ierr);
     169             :   }
     170           1 :   PetscFunctionReturn(0);
     171             : }
     172             : 
     173             : /*@
     174             :    NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
     175             : 
     176             :    Collective on nep
     177             : 
     178             :    Input Parameters:
     179             : +  nep    - the nonlinear eigensolver context
     180             : .  lambda - scalar argument
     181             : .  x      - vector to be multiplied against
     182             : -  v      - workspace vector (used only in the case of split form)
     183             : 
     184             :    Output Parameters:
     185             : +  y   - result vector
     186             : .  A   - Function matrix
     187             : -  B   - optional preconditioning matrix
     188             : 
     189             :    Note:
     190             :    If the nonlinear operator is represented in split form, the result
     191             :    y = T(lambda)*x is computed without building T(lambda) explicitly. In
     192             :    that case, parameters A and B are not used. Otherwise, the matrix
     193             :    T(lambda) is built and the effect is the same as a call to
     194             :    NEPComputeFunction() followed by a MatMult().
     195             : 
     196             :    Level: developer
     197             : 
     198             : .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
     199             : @*/
     200         700 : PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
     201             : {
     202         700 :   PetscErrorCode ierr;
     203         700 :   PetscInt       i;
     204         700 :   PetscScalar    alpha;
     205             : 
     206         700 :   PetscFunctionBegin;
     207         700 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     208        1400 :   PetscValidLogicalCollectiveScalar(nep,lambda,2);
     209         700 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     210         700 :   if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
     211         700 :   PetscValidHeaderSpecific(y,VEC_CLASSID,5);
     212         700 :   if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
     213         700 :   if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,7);
     214             : 
     215         700 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     216         644 :     ierr = VecSet(y,0.0);CHKERRQ(ierr);
     217        2534 :     for (i=0;i<nep->nt;i++) {
     218        1890 :       ierr = FNEvaluateFunction(nep->f[i],lambda,&alpha);CHKERRQ(ierr);
     219        1890 :       ierr = MatMult(nep->A[i],x,v);CHKERRQ(ierr);
     220        1890 :       ierr = VecAXPY(y,alpha,v);CHKERRQ(ierr);
     221             :     }
     222             :   } else {
     223          56 :     ierr = NEPComputeFunction(nep,lambda,A,B);CHKERRQ(ierr);
     224          56 :     ierr = MatMult(A,x,y);CHKERRQ(ierr);
     225             :   }
     226         700 :   PetscFunctionReturn(0);
     227             : }
     228             : 
     229             : /*@
     230             :    NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
     231             : 
     232             :    Collective on nep
     233             : 
     234             :    Input Parameters:
     235             : +  nep    - the nonlinear eigensolver context
     236             : .  lambda - scalar argument
     237             : .  x      - vector to be multiplied against
     238             : -  v      - workspace vector (used only in the case of split form)
     239             : 
     240             :    Output Parameters:
     241             : +  y   - result vector
     242             : .  A   - Function matrix
     243             : -  B   - optional preconditioning matrix
     244             : 
     245             :    Level: developer
     246             : 
     247             : .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
     248             : @*/
     249          17 : PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
     250             : {
     251          17 :   PetscErrorCode ierr;
     252          17 :   PetscInt       i;
     253          17 :   PetscScalar    alpha;
     254             : 
     255          17 :   PetscFunctionBegin;
     256          17 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     257          34 :   PetscValidLogicalCollectiveScalar(nep,lambda,2);
     258          17 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     259          17 :   if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
     260          17 :   PetscValidHeaderSpecific(y,VEC_CLASSID,5);
     261          17 :   if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
     262          17 :   if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,7);
     263             : 
     264          17 :   ierr = VecConjugate(x);CHKERRQ(ierr);
     265          17 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     266          11 :     ierr = VecSet(y,0.0);CHKERRQ(ierr);
     267          44 :     for (i=0;i<nep->nt;i++) {
     268          33 :       ierr = FNEvaluateFunction(nep->f[i],lambda,&alpha);CHKERRQ(ierr);
     269          33 :       ierr = MatMultTranspose(nep->A[i],x,v);CHKERRQ(ierr);
     270          33 :       ierr = VecAXPY(y,alpha,v);CHKERRQ(ierr);
     271             :     }
     272             :   } else {
     273           6 :     ierr = NEPComputeFunction(nep,lambda,A,B);CHKERRQ(ierr);
     274           6 :     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
     275             :   }
     276          17 :   ierr = VecConjugate(x);CHKERRQ(ierr);
     277          17 :   ierr = VecConjugate(y);CHKERRQ(ierr);
     278          17 :   PetscFunctionReturn(0);
     279             : }
     280             : 
     281             : /*@
     282             :    NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
     283             : 
     284             :    Collective on nep
     285             : 
     286             :    Input Parameters:
     287             : +  nep    - the nonlinear eigensolver context
     288             : .  lambda - scalar argument
     289             : .  x      - vector to be multiplied against
     290             : -  v      - workspace vector (used only in the case of split form)
     291             : 
     292             :    Output Parameters:
     293             : +  y   - result vector
     294             : -  A   - Jacobian matrix
     295             : 
     296             :    Note:
     297             :    If the nonlinear operator is represented in split form, the result
     298             :    y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
     299             :    that case, parameter A is not used. Otherwise, the matrix
     300             :    T'(lambda) is built and the effect is the same as a call to
     301             :    NEPComputeJacobian() followed by a MatMult().
     302             : 
     303             :    Level: developer
     304             : 
     305             : .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
     306             : @*/
     307          32 : PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
     308             : {
     309          32 :   PetscErrorCode ierr;
     310          32 :   PetscInt       i;
     311          32 :   PetscScalar    alpha;
     312             : 
     313          32 :   PetscFunctionBegin;
     314          32 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     315          64 :   PetscValidLogicalCollectiveScalar(nep,lambda,2);
     316          32 :   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
     317          32 :   if (v) PetscValidHeaderSpecific(v,VEC_CLASSID,4);
     318          32 :   PetscValidHeaderSpecific(y,VEC_CLASSID,5);
     319          32 :   if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,6);
     320             : 
     321          32 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     322          32 :     ierr = VecSet(y,0.0);CHKERRQ(ierr);
     323         128 :     for (i=0;i<nep->nt;i++) {
     324          96 :       ierr = FNEvaluateDerivative(nep->f[i],lambda,&alpha);CHKERRQ(ierr);
     325          96 :       ierr = MatMult(nep->A[i],x,v);CHKERRQ(ierr);
     326          96 :       ierr = VecAXPY(y,alpha,v);CHKERRQ(ierr);
     327             :     }
     328             :   } else {
     329           0 :     ierr = NEPComputeJacobian(nep,lambda,A);CHKERRQ(ierr);
     330           0 :     ierr = MatMult(A,x,y);CHKERRQ(ierr);
     331             :   }
     332          32 :   PetscFunctionReturn(0);
     333             : }
     334             : 
     335             : /*@
     336             :    NEPGetIterationNumber - Gets the current iteration number. If the
     337             :    call to NEPSolve() is complete, then it returns the number of iterations
     338             :    carried out by the solution method.
     339             : 
     340             :    Not Collective
     341             : 
     342             :    Input Parameter:
     343             : .  nep - the nonlinear eigensolver context
     344             : 
     345             :    Output Parameter:
     346             : .  its - number of iterations
     347             : 
     348             :    Note:
     349             :    During the i-th iteration this call returns i-1. If NEPSolve() is
     350             :    complete, then parameter "its" contains either the iteration number at
     351             :    which convergence was successfully reached, or failure was detected.
     352             :    Call NEPGetConvergedReason() to determine if the solver converged or
     353             :    failed and why.
     354             : 
     355             :    Level: intermediate
     356             : 
     357             : .seealso: NEPGetConvergedReason(), NEPSetTolerances()
     358             : @*/
     359           9 : PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
     360             : {
     361           9 :   PetscFunctionBegin;
     362           9 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     363           9 :   PetscValidIntPointer(its,2);
     364           9 :   *its = nep->its;
     365           9 :   PetscFunctionReturn(0);
     366             : }
     367             : 
     368             : /*@
     369             :    NEPGetConverged - Gets the number of converged eigenpairs.
     370             : 
     371             :    Not Collective
     372             : 
     373             :    Input Parameter:
     374             : .  nep - the nonlinear eigensolver context
     375             : 
     376             :    Output Parameter:
     377             : .  nconv - number of converged eigenpairs
     378             : 
     379             :    Note:
     380             :    This function should be called after NEPSolve() has finished.
     381             : 
     382             :    Level: beginner
     383             : 
     384             : .seealso: NEPSetDimensions(), NEPSolve()
     385             : @*/
     386          38 : PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
     387             : {
     388          38 :   PetscFunctionBegin;
     389          38 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     390          38 :   PetscValidIntPointer(nconv,2);
     391          38 :   NEPCheckSolved(nep,1);
     392          38 :   *nconv = nep->nconv;
     393          38 :   PetscFunctionReturn(0);
     394             : }
     395             : 
     396             : /*@
     397             :    NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
     398             :    stopped.
     399             : 
     400             :    Not Collective
     401             : 
     402             :    Input Parameter:
     403             : .  nep - the nonlinear eigensolver context
     404             : 
     405             :    Output Parameter:
     406             : .  reason - negative value indicates diverged, positive value converged
     407             : 
     408             :    Notes:
     409             : 
     410             :    Possible values for reason are
     411             : +  NEP_CONVERGED_TOL - converged up to tolerance
     412             : .  NEP_CONVERGED_USER - converged due to a user-defined condition
     413             : .  NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
     414             : .  NEP_DIVERGED_BREAKDOWN - generic breakdown in method
     415             : .  NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
     416             : -  NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
     417             :    unrestarted solver
     418             : 
     419             :    Can only be called after the call to NEPSolve() is complete.
     420             : 
     421             :    Level: intermediate
     422             : 
     423             : .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
     424             : @*/
     425           1 : PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
     426             : {
     427           1 :   PetscFunctionBegin;
     428           1 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     429           1 :   PetscValidPointer(reason,2);
     430           1 :   NEPCheckSolved(nep,1);
     431           1 :   *reason = nep->reason;
     432           1 :   PetscFunctionReturn(0);
     433             : }
     434             : 
     435             : /*@C
     436             :    NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
     437             :    NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
     438             : 
     439             :    Logically Collective on nep
     440             : 
     441             :    Input Parameters:
     442             : +  nep - nonlinear eigensolver context
     443             : -  i   - index of the solution
     444             : 
     445             :    Output Parameters:
     446             : +  eigr - real part of eigenvalue
     447             : .  eigi - imaginary part of eigenvalue
     448             : .  Vr   - real part of eigenvector
     449             : -  Vi   - imaginary part of eigenvector
     450             : 
     451             :    Notes:
     452             :    It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
     453             :    required. Otherwise, the caller must provide valid Vec objects, i.e.,
     454             :    they must be created by the calling program with e.g. MatCreateVecs().
     455             : 
     456             :    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
     457             :    configured with complex scalars the eigenvalue is stored
     458             :    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
     459             :    set to zero). In any case, the user can pass NULL in Vr or Vi if one of
     460             :    them is not required.
     461             : 
     462             :    The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
     463             :    Eigenpairs are indexed according to the ordering criterion established
     464             :    with NEPSetWhichEigenpairs().
     465             : 
     466             :    Level: beginner
     467             : 
     468             : .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
     469             : @*/
     470         657 : PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
     471             : {
     472         657 :   PetscInt       k;
     473         657 :   PetscErrorCode ierr;
     474             : 
     475         657 :   PetscFunctionBegin;
     476         657 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     477        1314 :   PetscValidLogicalCollectiveInt(nep,i,2);
     478         657 :   if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,5); PetscCheckSameComm(nep,1,Vr,5); }
     479         657 :   if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,6); PetscCheckSameComm(nep,1,Vi,6); }
     480         657 :   NEPCheckSolved(nep,1);
     481         657 :   if (i<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     482         657 :   if (i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
     483             : 
     484         657 :   ierr = NEPComputeVectors(nep);CHKERRQ(ierr);
     485         657 :   k = nep->perm[i];
     486             : 
     487             :   /* eigenvalue */
     488             : #if defined(PETSC_USE_COMPLEX)
     489         657 :   if (eigr) *eigr = nep->eigr[k];
     490         657 :   if (eigi) *eigi = 0;
     491             : #else
     492             :   if (eigr) *eigr = nep->eigr[k];
     493             :   if (eigi) *eigi = nep->eigi[k];
     494             : #endif
     495             : 
     496             :   /* eigenvector */
     497         657 :   ierr = BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);CHKERRQ(ierr);
     498         657 :   PetscFunctionReturn(0);
     499             : }
     500             : 
     501             : /*@
     502             :    NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
     503             : 
     504             :    Logically Collective on nep
     505             : 
     506             :    Input Parameters:
     507             : +  nep - eigensolver context
     508             : -  i   - index of the solution
     509             : 
     510             :    Output Parameters:
     511             : +  Wr   - real part of left eigenvector
     512             : -  Wi   - imaginary part of left eigenvector
     513             : 
     514             :    Notes:
     515             :    The caller must provide valid Vec objects, i.e., they must be created
     516             :    by the calling program with e.g. MatCreateVecs().
     517             : 
     518             :    If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
     519             :    configured with complex scalars the eigenvector is stored directly in Wr
     520             :    (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
     521             :    them is not required.
     522             : 
     523             :    The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
     524             :    Eigensolutions are indexed according to the ordering criterion established
     525             :    with NEPSetWhichEigenpairs().
     526             : 
     527             :    Left eigenvectors are available only if the twosided flag was set, see
     528             :    NEPSetTwoSided().
     529             : 
     530             :    Level: intermediate
     531             : 
     532             : .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
     533             : @*/
     534          17 : PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
     535             : {
     536          17 :   PetscErrorCode ierr;
     537          17 :   PetscInt       k;
     538             : 
     539          17 :   PetscFunctionBegin;
     540          17 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     541          34 :   PetscValidLogicalCollectiveInt(nep,i,2);
     542          17 :   if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(nep,1,Wr,3); }
     543          17 :   if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(nep,1,Wi,4); }
     544          17 :   NEPCheckSolved(nep,1);
     545          17 :   if (!nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
     546          17 :   if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
     547          17 :   ierr = NEPComputeVectors(nep);CHKERRQ(ierr);
     548          17 :   k = nep->perm[i];
     549          17 :   ierr = BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);CHKERRQ(ierr);
     550          17 :   PetscFunctionReturn(0);
     551             : }
     552             : 
     553             : /*@
     554             :    NEPGetErrorEstimate - Returns the error estimate associated to the i-th
     555             :    computed eigenpair.
     556             : 
     557             :    Not Collective
     558             : 
     559             :    Input Parameters:
     560             : +  nep - nonlinear eigensolver context
     561             : -  i   - index of eigenpair
     562             : 
     563             :    Output Parameter:
     564             : .  errest - the error estimate
     565             : 
     566             :    Notes:
     567             :    This is the error estimate used internally by the eigensolver. The actual
     568             :    error bound can be computed with NEPComputeError().
     569             : 
     570             :    Level: advanced
     571             : 
     572             : .seealso: NEPComputeError()
     573             : @*/
     574           4 : PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
     575             : {
     576           4 :   PetscFunctionBegin;
     577           4 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     578           4 :   PetscValidRealPointer(errest,3);
     579           4 :   NEPCheckSolved(nep,1);
     580           4 :   if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
     581           4 :   *errest = nep->errest[nep->perm[i]];
     582           4 :   PetscFunctionReturn(0);
     583             : }
     584             : 
     585             : /*
     586             :    NEPComputeResidualNorm_Private - Computes the norm of the residual vector
     587             :    associated with an eigenpair.
     588             : 
     589             :    Input Parameters:
     590             :      adj    - whether the adjoint T^* must be used instead of T
     591             :      lambda - eigenvalue
     592             :      x      - eigenvector
     593             :      w      - array of work vectors (two vectors in split form, one vector otherwise)
     594             : */
     595         709 : PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
     596             : {
     597         709 :   PetscErrorCode ierr;
     598         709 :   Vec            y,z=NULL;
     599             : 
     600         709 :   PetscFunctionBegin;
     601         709 :   y = w[0];
     602         709 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
     603         709 :   if (adj) {
     604           9 :     ierr = NEPApplyAdjoint(nep,lambda,x,z,y,nep->function,nep->function_pre);CHKERRQ(ierr);
     605             :   } else {
     606         700 :     ierr = NEPApplyFunction(nep,lambda,x,z,y,nep->function,nep->function_pre);CHKERRQ(ierr);
     607             :   }
     608         709 :   ierr = VecNorm(y,NORM_2,norm);CHKERRQ(ierr);
     609         709 :   PetscFunctionReturn(0);
     610             : }
     611             : 
     612             : /*@
     613             :    NEPComputeError - Computes the error (based on the residual norm) associated
     614             :    with the i-th computed eigenpair.
     615             : 
     616             :    Collective on nep
     617             : 
     618             :    Input Parameters:
     619             : +  nep  - the nonlinear eigensolver context
     620             : .  i    - the solution index
     621             : -  type - the type of error to compute
     622             : 
     623             :    Output Parameter:
     624             : .  error - the error
     625             : 
     626             :    Notes:
     627             :    The error can be computed in various ways, all of them based on the residual
     628             :    norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
     629             :    eigenvector.
     630             : 
     631             :    Level: beginner
     632             : 
     633             : .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
     634             : @*/
     635         593 : PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
     636             : {
     637         593 :   PetscErrorCode ierr;
     638         593 :   Vec            xr,xi=NULL;
     639         593 :   PetscInt       j,nwork,issplit=0;
     640         593 :   PetscScalar    kr,ki,s;
     641         593 :   PetscReal      er,z=0.0,errorl,nrm;
     642         593 :   PetscBool      flg;
     643             : 
     644         593 :   PetscFunctionBegin;
     645         593 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     646        1186 :   PetscValidLogicalCollectiveInt(nep,i,2);
     647        1186 :   PetscValidLogicalCollectiveEnum(nep,type,3);
     648         593 :   PetscValidRealPointer(error,4);
     649         593 :   NEPCheckSolved(nep,1);
     650             : 
     651             :   /* allocate work vectors */
     652             : #if defined(PETSC_USE_COMPLEX)
     653         593 :   nwork = 2;
     654             : #else
     655             :   nwork = 3;
     656             : #endif
     657         593 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     658         549 :     issplit = 1;
     659         549 :     nwork++;  /* need an extra work vector for NEPComputeResidualNorm_Private */
     660             :   }
     661         593 :   ierr = NEPSetWorkVecs(nep,nwork);CHKERRQ(ierr);
     662         593 :   xr = nep->work[issplit+1];
     663             : #if !defined(PETSC_USE_COMPLEX)
     664             :   xi = nep->work[issplit+2];
     665             : #endif
     666             : 
     667             :   /* compute residual norms */
     668         593 :   ierr = NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
     669             : #if !defined(PETSC_USE_COMPLEX)
     670             :   if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
     671             : #endif
     672         593 :   ierr = NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);CHKERRQ(ierr);
     673         593 :   ierr = VecNorm(xr,NORM_2,&er);CHKERRQ(ierr);
     674             : 
     675             :   /* if two-sided, compute left residual norm and take the maximum */
     676         593 :   if (nep->twosided) {
     677           9 :     ierr = NEPGetLeftEigenvector(nep,i,xr,xi);CHKERRQ(ierr);
     678           9 :     ierr = NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl);CHKERRQ(ierr);
     679          15 :     *error = PetscMax(*error,errorl);
     680             :   }
     681             : 
     682             :   /* compute error */
     683         593 :   switch (type) {
     684             :     case NEP_ERROR_ABSOLUTE:
     685             :       break;
     686         524 :     case NEP_ERROR_RELATIVE:
     687         524 :       *error /= PetscAbsScalar(kr)*er;
     688         524 :       break;
     689          68 :     case NEP_ERROR_BACKWARD:
     690          68 :       if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
     691          18 :         ierr = NEPComputeFunction(nep,kr,nep->function,nep->function);CHKERRQ(ierr);
     692          18 :         ierr = MatHasOperation(nep->function,MATOP_NORM,&flg);CHKERRQ(ierr);
     693          18 :         if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
     694          18 :         ierr = MatNorm(nep->function,NORM_INFINITY,&nrm);CHKERRQ(ierr);
     695          18 :         *error /= nrm*er;
     696          18 :         break;
     697             :       }
     698             :       /* initialization of matrix norms */
     699          50 :       if (!nep->nrma[0]) {
     700          46 :         for (j=0;j<nep->nt;j++) {
     701          32 :           ierr = MatHasOperation(nep->A[j],MATOP_NORM,&flg);CHKERRQ(ierr);
     702          32 :           if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
     703          32 :           ierr = MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);CHKERRQ(ierr);
     704             :         }
     705             :       }
     706         161 :       for (j=0;j<nep->nt;j++) {
     707         111 :         ierr = FNEvaluateFunction(nep->f[j],kr,&s);CHKERRQ(ierr);
     708         111 :         z = z + nep->nrma[j]*PetscAbsScalar(s);
     709             :       }
     710          50 :       *error /= z*er;
     711          50 :       break;
     712           0 :     default:
     713           0 :       SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
     714             :   }
     715         593 :   PetscFunctionReturn(0);
     716             : }
     717             : 
     718             : /*@
     719             :    NEPComputeFunction - Computes the function matrix T(lambda) that has been
     720             :    set with NEPSetFunction().
     721             : 
     722             :    Collective on nep
     723             : 
     724             :    Input Parameters:
     725             : +  nep    - the NEP context
     726             : -  lambda - the scalar argument
     727             : 
     728             :    Output Parameters:
     729             : +  A   - Function matrix
     730             : -  B   - optional preconditioning matrix
     731             : 
     732             :    Notes:
     733             :    NEPComputeFunction() is typically used within nonlinear eigensolvers
     734             :    implementations, so most users would not generally call this routine
     735             :    themselves.
     736             : 
     737             :    Level: developer
     738             : 
     739             : .seealso: NEPSetFunction(), NEPGetFunction()
     740             : @*/
     741        2052 : PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
     742             : {
     743        2052 :   PetscErrorCode ierr;
     744        2052 :   PetscInt       i;
     745        2052 :   PetscScalar    alpha;
     746             : 
     747        2052 :   PetscFunctionBegin;
     748        2052 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     749        2052 :   NEPCheckProblem(nep,1);
     750        2052 :   switch (nep->fui) {
     751         799 :   case NEP_USER_INTERFACE_CALLBACK:
     752         799 :     if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
     753         799 :     ierr = PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);CHKERRQ(ierr);
     754         799 :     PetscStackPush("NEP user Function function");
     755         799 :     ierr = (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);CHKERRQ(ierr);
     756         799 :     PetscStackPop;
     757         799 :     ierr = PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);CHKERRQ(ierr);
     758             :     break;
     759        1253 :   case NEP_USER_INTERFACE_SPLIT:
     760        1253 :     ierr = MatZeroEntries(A);CHKERRQ(ierr);
     761        4942 :     for (i=0;i<nep->nt;i++) {
     762        3689 :       ierr = FNEvaluateFunction(nep->f[i],lambda,&alpha);CHKERRQ(ierr);
     763        3689 :       ierr = MatAXPY(A,alpha,nep->A[i],nep->mstr);CHKERRQ(ierr);
     764             :     }
     765        1253 :     if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented");
     766             :     break;
     767             :   }
     768        2052 :   PetscFunctionReturn(0);
     769             : }
     770             : 
     771             : /*@
     772             :    NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
     773             :    set with NEPSetJacobian().
     774             : 
     775             :    Collective on nep
     776             : 
     777             :    Input Parameters:
     778             : +  nep    - the NEP context
     779             : -  lambda - the scalar argument
     780             : 
     781             :    Output Parameters:
     782             : .  A   - Jacobian matrix
     783             : 
     784             :    Notes:
     785             :    Most users should not need to explicitly call this routine, as it
     786             :    is used internally within the nonlinear eigensolvers.
     787             : 
     788             :    Level: developer
     789             : 
     790             : .seealso: NEPSetJacobian(), NEPGetJacobian()
     791             : @*/
     792        1239 : PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
     793             : {
     794        1239 :   PetscErrorCode ierr;
     795        1239 :   PetscInt       i;
     796        1239 :   PetscScalar    alpha;
     797             : 
     798        1239 :   PetscFunctionBegin;
     799        1239 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     800        1239 :   NEPCheckProblem(nep,1);
     801        1239 :   switch (nep->fui) {
     802         255 :   case NEP_USER_INTERFACE_CALLBACK:
     803         255 :     if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
     804         255 :     ierr = PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);CHKERRQ(ierr);
     805         255 :     PetscStackPush("NEP user Jacobian function");
     806         255 :     ierr = (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);CHKERRQ(ierr);
     807         255 :     PetscStackPop;
     808         255 :     ierr = PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);CHKERRQ(ierr);
     809             :     break;
     810         984 :   case NEP_USER_INTERFACE_SPLIT:
     811         984 :     ierr = MatZeroEntries(A);CHKERRQ(ierr);
     812        3868 :     for (i=0;i<nep->nt;i++) {
     813        2884 :       ierr = FNEvaluateDerivative(nep->f[i],lambda,&alpha);CHKERRQ(ierr);
     814        2884 :       ierr = MatAXPY(A,alpha,nep->A[i],nep->mstr);CHKERRQ(ierr);
     815             :     }
     816             :     break;
     817             :   }
     818        1239 :   PetscFunctionReturn(0);
     819             : }
     820             : 

Generated by: LCOV version 1.14