LCOV - code coverage report
Current view: top level - eps/interface - epssolve.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 312 327 95.4 %
Date: 2024-03-29 00:55:19 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-, 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             :    EPS routines related to the solution process
      12             : */
      13             : 
      14             : #include <slepc/private/epsimpl.h>   /*I "slepceps.h" I*/
      15             : #include <slepc/private/bvimpl.h>
      16             : #include <petscdraw.h>
      17             : 
      18        6695 : PetscErrorCode EPSComputeVectors(EPS eps)
      19             : {
      20        6695 :   PetscFunctionBegin;
      21        6695 :   EPSCheckSolved(eps,1);
      22        6695 :   if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
      23        6695 :   eps->state = EPS_STATE_EIGENVECTORS;
      24        6695 :   PetscFunctionReturn(PETSC_SUCCESS);
      25             : }
      26             : 
      27             : #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
      28             : 
      29         855 : static PetscErrorCode EPSComputeValues(EPS eps)
      30             : {
      31         855 :   PetscBool      injective,iscomp,isfilter;
      32         855 :   PetscInt       i,n,aux,nconv0;
      33         855 :   Mat            A,B=NULL,G,Z;
      34             : 
      35         855 :   PetscFunctionBegin;
      36         855 :   switch (eps->categ) {
      37         639 :     case EPS_CATEGORY_KRYLOV:
      38             :     case EPS_CATEGORY_OTHER:
      39         639 :       PetscCall(STIsInjective(eps->st,&injective));
      40         639 :       if (injective) {
      41             :         /* one-to-one mapping: backtransform eigenvalues */
      42         629 :         PetscUseTypeMethod(eps,backtransform);
      43             :       } else {
      44             :         /* compute eigenvalues from Rayleigh quotient */
      45          10 :         PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
      46          10 :         if (!n) break;
      47          10 :         PetscCall(EPSGetOperators(eps,&A,&B));
      48          10 :         PetscCall(BVSetActiveColumns(eps->V,0,n));
      49          10 :         PetscCall(DSGetCompact(eps->ds,&iscomp));
      50          10 :         PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
      51          10 :         PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
      52          10 :         PetscCall(BVMatProject(eps->V,A,eps->V,G));
      53          10 :         PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
      54          10 :         if (B) {
      55           0 :           PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
      56           0 :           PetscCall(BVMatProject(eps->V,B,eps->V,G));
      57           0 :           PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
      58             :         }
      59          10 :         PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
      60          10 :         PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
      61          10 :         PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
      62          10 :         PetscCall(DSSetCompact(eps->ds,iscomp));
      63          10 :         if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
      64          10 :           PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
      65          10 :           PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
      66          10 :           PetscCall(BVMultInPlace(eps->V,Z,0,n));
      67          10 :           PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
      68             :         }
      69             :         /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
      70          10 :         PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
      71          10 :         if (isfilter) {
      72           9 :           nconv0 = eps->nconv;
      73          82 :           for (i=0;i<eps->nconv;i++) {
      74          73 :             if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
      75           0 :               eps->nconv--;
      76           0 :               if (i<eps->nconv) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
      77             :             }
      78             :           }
      79           9 :           if (nconv0>eps->nconv) PetscCall(PetscInfo(eps,"Discarded %" PetscInt_FMT " computed eigenvalues lying outside the interval\n",nconv0-eps->nconv));
      80             :         }
      81             :       }
      82             :       break;
      83             :     case EPS_CATEGORY_PRECOND:
      84             :     case EPS_CATEGORY_CONTOUR:
      85             :       /* eigenvalues already available as an output of the solver */
      86             :       break;
      87             :   }
      88         855 :   PetscFunctionReturn(PETSC_SUCCESS);
      89             : }
      90             : 
      91             : /*@
      92             :    EPSSolve - Solves the eigensystem.
      93             : 
      94             :    Collective
      95             : 
      96             :    Input Parameter:
      97             : .  eps - eigensolver context obtained from EPSCreate()
      98             : 
      99             :    Options Database Keys:
     100             : +  -eps_view - print information about the solver used
     101             : .  -eps_view_mat0 - view the first matrix (A)
     102             : .  -eps_view_mat1 - view the second matrix (B)
     103             : .  -eps_view_vectors - view the computed eigenvectors
     104             : .  -eps_view_values - view the computed eigenvalues
     105             : .  -eps_converged_reason - print reason for convergence, and number of iterations
     106             : .  -eps_error_absolute - print absolute errors of each eigenpair
     107             : .  -eps_error_relative - print relative errors of each eigenpair
     108             : -  -eps_error_backward - print backward errors of each eigenpair
     109             : 
     110             :    Notes:
     111             :    All the command-line options listed above admit an optional argument specifying
     112             :    the viewer type and options. For instance, use '-eps_view_mat0 binary:amatrix.bin'
     113             :    to save the A matrix to a binary file, '-eps_view_values draw' to draw the computed
     114             :    eigenvalues graphically, or '-eps_error_relative :myerr.m:ascii_matlab' to save
     115             :    the errors in a file that can be executed in Matlab.
     116             : 
     117             :    Level: beginner
     118             : 
     119             : .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
     120             : @*/
     121         855 : PetscErrorCode EPSSolve(EPS eps)
     122             : {
     123         855 :   PetscInt       i;
     124         855 :   PetscBool      hasname;
     125         855 :   STMatMode      matmode;
     126         855 :   Mat            A,B;
     127             : 
     128         855 :   PetscFunctionBegin;
     129         855 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     130         855 :   if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
     131         855 :   PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
     132             : 
     133             :   /* Call setup */
     134         855 :   PetscCall(EPSSetUp(eps));
     135         855 :   eps->nconv = 0;
     136         855 :   eps->its   = 0;
     137       27537 :   for (i=0;i<eps->ncv;i++) {
     138       26682 :     eps->eigr[i]   = 0.0;
     139       26682 :     eps->eigi[i]   = 0.0;
     140       26682 :     eps->errest[i] = 0.0;
     141       26682 :     eps->perm[i]   = i;
     142             :   }
     143         855 :   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
     144         855 :   PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
     145             : 
     146             :   /* Call solver */
     147         855 :   PetscUseTypeMethod(eps,solve);
     148         855 :   PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
     149         855 :   eps->state = EPS_STATE_SOLVED;
     150             : 
     151             :   /* Only the first nconv columns contain useful information (except in CISS) */
     152         855 :   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
     153         855 :   if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
     154             : 
     155             :   /* If inplace, purify eigenvectors before reverting operator */
     156         855 :   PetscCall(STGetMatMode(eps->st,&matmode));
     157         855 :   if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
     158         855 :   PetscCall(STPostSolve(eps->st));
     159             : 
     160             :   /* Map eigenvalues back to the original problem if appropriate */
     161         855 :   PetscCall(EPSComputeValues(eps));
     162             : 
     163             : #if !defined(PETSC_USE_COMPLEX)
     164             :   /* Reorder conjugate eigenvalues (positive imaginary first) */
     165        9087 :   for (i=0;i<eps->nconv-1;i++) {
     166        8232 :     if (eps->eigi[i] != 0) {
     167         121 :       if (eps->eigi[i] < 0) {
     168          21 :         eps->eigi[i] = -eps->eigi[i];
     169          21 :         eps->eigi[i+1] = -eps->eigi[i+1];
     170             :         /* the next correction only works with eigenvectors */
     171          21 :         PetscCall(EPSComputeVectors(eps));
     172          21 :         PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
     173             :       }
     174         121 :       i++;
     175             :     }
     176             :   }
     177             : #endif
     178             : 
     179             :   /* Sort eigenvalues according to eps->which parameter */
     180         855 :   PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
     181         855 :   PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
     182             : 
     183             :   /* Various viewers */
     184         855 :   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
     185         855 :   PetscCall(EPSConvergedReasonViewFromOptions(eps));
     186         855 :   PetscCall(EPSErrorViewFromOptions(eps));
     187         855 :   PetscCall(EPSValuesViewFromOptions(eps));
     188         855 :   PetscCall(EPSVectorsViewFromOptions(eps));
     189             : 
     190         855 :   PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
     191         855 :   if (hasname) {
     192           0 :     PetscCall(EPSGetOperators(eps,&A,NULL));
     193           0 :     PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
     194             :   }
     195         855 :   if (eps->isgeneralized) {
     196         210 :     PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
     197         210 :     if (hasname) {
     198           0 :       PetscCall(EPSGetOperators(eps,NULL,&B));
     199           0 :       PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
     200             :     }
     201             :   }
     202             : 
     203             :   /* Remove deflation and initial subspaces */
     204         855 :   if (eps->nds) {
     205          14 :     PetscCall(BVSetNumConstraints(eps->V,0));
     206          14 :     eps->nds = 0;
     207             :   }
     208         855 :   eps->nini = 0;
     209         855 :   PetscFunctionReturn(PETSC_SUCCESS);
     210             : }
     211             : 
     212             : /*@
     213             :    EPSGetIterationNumber - Gets the current iteration number. If the
     214             :    call to EPSSolve() is complete, then it returns the number of iterations
     215             :    carried out by the solution method.
     216             : 
     217             :    Not Collective
     218             : 
     219             :    Input Parameter:
     220             : .  eps - the eigensolver context
     221             : 
     222             :    Output Parameter:
     223             : .  its - number of iterations
     224             : 
     225             :    Note:
     226             :    During the i-th iteration this call returns i-1. If EPSSolve() is
     227             :    complete, then parameter "its" contains either the iteration number at
     228             :    which convergence was successfully reached, or failure was detected.
     229             :    Call EPSGetConvergedReason() to determine if the solver converged or
     230             :    failed and why.
     231             : 
     232             :    Level: intermediate
     233             : 
     234             : .seealso: EPSGetConvergedReason(), EPSSetTolerances()
     235             : @*/
     236         146 : PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
     237             : {
     238         146 :   PetscFunctionBegin;
     239         146 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     240         146 :   PetscAssertPointer(its,2);
     241         146 :   *its = eps->its;
     242         146 :   PetscFunctionReturn(PETSC_SUCCESS);
     243             : }
     244             : 
     245             : /*@
     246             :    EPSGetConverged - Gets the number of converged eigenpairs.
     247             : 
     248             :    Not Collective
     249             : 
     250             :    Input Parameter:
     251             : .  eps - the eigensolver context
     252             : 
     253             :    Output Parameter:
     254             : .  nconv - number of converged eigenpairs
     255             : 
     256             :    Note:
     257             :    This function should be called after EPSSolve() has finished.
     258             : 
     259             :    Level: beginner
     260             : 
     261             : .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
     262             : @*/
     263         501 : PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
     264             : {
     265         501 :   PetscFunctionBegin;
     266         501 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     267         501 :   PetscAssertPointer(nconv,2);
     268         501 :   EPSCheckSolved(eps,1);
     269         501 :   *nconv = eps->nconv;
     270         501 :   PetscFunctionReturn(PETSC_SUCCESS);
     271             : }
     272             : 
     273             : /*@
     274             :    EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
     275             :    stopped.
     276             : 
     277             :    Not Collective
     278             : 
     279             :    Input Parameter:
     280             : .  eps - the eigensolver context
     281             : 
     282             :    Output Parameter:
     283             : .  reason - negative value indicates diverged, positive value converged
     284             : 
     285             :    Options Database Key:
     286             : .  -eps_converged_reason - print the reason to a viewer
     287             : 
     288             :    Notes:
     289             :    Possible values for reason are
     290             : +  EPS_CONVERGED_TOL - converged up to tolerance
     291             : .  EPS_CONVERGED_USER - converged due to a user-defined condition
     292             : .  EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
     293             : .  EPS_DIVERGED_BREAKDOWN - generic breakdown in method
     294             : -  EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
     295             : 
     296             :    Can only be called after the call to EPSSolve() is complete.
     297             : 
     298             :    Level: intermediate
     299             : 
     300             : .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
     301             : @*/
     302         143 : PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
     303             : {
     304         143 :   PetscFunctionBegin;
     305         143 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     306         143 :   PetscAssertPointer(reason,2);
     307         143 :   EPSCheckSolved(eps,1);
     308         143 :   *reason = eps->reason;
     309         143 :   PetscFunctionReturn(PETSC_SUCCESS);
     310             : }
     311             : 
     312             : /*@
     313             :    EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
     314             :    subspace.
     315             : 
     316             :    Collective
     317             : 
     318             :    Input Parameter:
     319             : .  eps - the eigensolver context
     320             : 
     321             :    Output Parameter:
     322             : .  v - an array of vectors
     323             : 
     324             :    Notes:
     325             :    This function should be called after EPSSolve() has finished.
     326             : 
     327             :    The user should provide in v an array of nconv vectors, where nconv is
     328             :    the value returned by EPSGetConverged().
     329             : 
     330             :    The first k vectors returned in v span an invariant subspace associated
     331             :    with the first k computed eigenvalues (note that this is not true if the
     332             :    k-th eigenvalue is complex and matrix A is real; in this case the first
     333             :    k+1 vectors should be used). An invariant subspace X of A satisfies Ax
     334             :    in X for all x in X (a similar definition applies for generalized
     335             :    eigenproblems).
     336             : 
     337             :    Level: intermediate
     338             : 
     339             : .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
     340             : @*/
     341           8 : PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
     342             : {
     343           8 :   PetscInt       i;
     344           8 :   BV             V=eps->V;
     345           8 :   Vec            w;
     346             : 
     347           8 :   PetscFunctionBegin;
     348           8 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     349           8 :   PetscAssertPointer(v,2);
     350           8 :   PetscValidHeaderSpecific(*v,VEC_CLASSID,2);
     351           8 :   EPSCheckSolved(eps,1);
     352           8 :   PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
     353           8 :   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
     354           1 :     PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
     355           1 :     PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
     356           1 :     PetscCall(BVCopy(eps->V,V));
     357           5 :     for (i=0;i<eps->nconv;i++) {
     358           4 :       PetscCall(BVGetColumn(V,i,&w));
     359           4 :       PetscCall(VecPointwiseDivide(w,w,eps->D));
     360           4 :       PetscCall(BVRestoreColumn(V,i,&w));
     361             :     }
     362           1 :     PetscCall(BVOrthogonalize(V,NULL));
     363             :   }
     364          60 :   for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
     365           8 :   if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
     366           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     367             : }
     368             : 
     369             : /*@C
     370             :    EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
     371             :    EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
     372             : 
     373             :    Collective
     374             : 
     375             :    Input Parameters:
     376             : +  eps - eigensolver context
     377             : -  i   - index of the solution
     378             : 
     379             :    Output Parameters:
     380             : +  eigr - real part of eigenvalue
     381             : .  eigi - imaginary part of eigenvalue
     382             : .  Vr   - real part of eigenvector
     383             : -  Vi   - imaginary part of eigenvector
     384             : 
     385             :    Notes:
     386             :    It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
     387             :    required. Otherwise, the caller must provide valid Vec objects, i.e.,
     388             :    they must be created by the calling program with e.g. MatCreateVecs().
     389             : 
     390             :    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
     391             :    configured with complex scalars the eigenvalue is stored
     392             :    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
     393             :    set to zero). In both cases, the user can pass NULL in eigi and Vi.
     394             : 
     395             :    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
     396             :    Eigenpairs are indexed according to the ordering criterion established
     397             :    with EPSSetWhichEigenpairs().
     398             : 
     399             :    The 2-norm of the eigenvector is one unless the problem is generalized
     400             :    Hermitian. In this case the eigenvector is normalized with respect to the
     401             :    norm defined by the B matrix.
     402             : 
     403             :    Level: beginner
     404             : 
     405             : .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
     406             :           EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
     407             : @*/
     408        5100 : PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
     409             : {
     410        5100 :   PetscFunctionBegin;
     411        5100 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     412       20400 :   PetscValidLogicalCollectiveInt(eps,i,2);
     413        5100 :   EPSCheckSolved(eps,1);
     414        5100 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     415        5100 :   PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
     416        5100 :   PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
     417        5100 :   if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
     418        5100 :   PetscFunctionReturn(PETSC_SUCCESS);
     419             : }
     420             : 
     421             : /*@C
     422             :    EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
     423             : 
     424             :    Not Collective
     425             : 
     426             :    Input Parameters:
     427             : +  eps - eigensolver context
     428             : -  i   - index of the solution
     429             : 
     430             :    Output Parameters:
     431             : +  eigr - real part of eigenvalue
     432             : -  eigi - imaginary part of eigenvalue
     433             : 
     434             :    Notes:
     435             :    If the eigenvalue is real, then eigi is set to zero. If PETSc is
     436             :    configured with complex scalars the eigenvalue is stored
     437             :    directly in eigr (eigi is set to zero).
     438             : 
     439             :    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
     440             :    Eigenpairs are indexed according to the ordering criterion established
     441             :    with EPSSetWhichEigenpairs().
     442             : 
     443             :    Level: beginner
     444             : 
     445             : .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
     446             : @*/
     447        6705 : PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
     448             : {
     449        6705 :   PetscInt k;
     450             : 
     451        6705 :   PetscFunctionBegin;
     452        6705 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     453        6705 :   EPSCheckSolved(eps,1);
     454        6705 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     455        6705 :   PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
     456        6705 :   k = eps->perm[i];
     457             : #if defined(PETSC_USE_COMPLEX)
     458             :   if (eigr) *eigr = eps->eigr[k];
     459             :   if (eigi) *eigi = 0;
     460             : #else
     461        6705 :   if (eigr) *eigr = eps->eigr[k];
     462        6705 :   if (eigi) *eigi = eps->eigi[k];
     463             : #endif
     464        6705 :   PetscFunctionReturn(PETSC_SUCCESS);
     465             : }
     466             : 
     467             : /*@
     468             :    EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
     469             : 
     470             :    Collective
     471             : 
     472             :    Input Parameters:
     473             : +  eps - eigensolver context
     474             : -  i   - index of the solution
     475             : 
     476             :    Output Parameters:
     477             : +  Vr   - real part of eigenvector
     478             : -  Vi   - imaginary part of eigenvector
     479             : 
     480             :    Notes:
     481             :    The caller must provide valid Vec objects, i.e., they must be created
     482             :    by the calling program with e.g. MatCreateVecs().
     483             : 
     484             :    If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
     485             :    configured with complex scalars the eigenvector is stored
     486             :    directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
     487             :    or Vi if one of them is not required.
     488             : 
     489             :    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
     490             :    Eigenpairs are indexed according to the ordering criterion established
     491             :    with EPSSetWhichEigenpairs().
     492             : 
     493             :    The 2-norm of the eigenvector is one unless the problem is generalized
     494             :    Hermitian. In this case the eigenvector is normalized with respect to the
     495             :    norm defined by the B matrix.
     496             : 
     497             :    Level: beginner
     498             : 
     499             : .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
     500             : @*/
     501        6570 : PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
     502             : {
     503        6570 :   PetscInt       k;
     504             : 
     505        6570 :   PetscFunctionBegin;
     506        6570 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     507       26280 :   PetscValidLogicalCollectiveInt(eps,i,2);
     508        6570 :   if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Vr,3); }
     509        6570 :   if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Vi,4); }
     510        6570 :   EPSCheckSolved(eps,1);
     511        6570 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     512        6570 :   PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
     513        6570 :   PetscCall(EPSComputeVectors(eps));
     514        6570 :   k = eps->perm[i];
     515        6570 :   PetscCall(BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi));
     516        6570 :   PetscFunctionReturn(PETSC_SUCCESS);
     517             : }
     518             : 
     519             : /*@
     520             :    EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
     521             : 
     522             :    Collective
     523             : 
     524             :    Input Parameters:
     525             : +  eps - eigensolver context
     526             : -  i   - index of the solution
     527             : 
     528             :    Output Parameters:
     529             : +  Wr   - real part of left eigenvector
     530             : -  Wi   - imaginary part of left eigenvector
     531             : 
     532             :    Notes:
     533             :    The caller must provide valid Vec objects, i.e., they must be created
     534             :    by the calling program with e.g. MatCreateVecs().
     535             : 
     536             :    If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
     537             :    configured with complex scalars the eigenvector is stored directly in Wr
     538             :    (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
     539             :    one of them is not required.
     540             : 
     541             :    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
     542             :    Eigensolutions are indexed according to the ordering criterion established
     543             :    with EPSSetWhichEigenpairs().
     544             : 
     545             :    Left eigenvectors are available only if the twosided flag was set, see
     546             :    EPSSetTwoSided().
     547             : 
     548             :    Level: intermediate
     549             : 
     550             : .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
     551             : @*/
     552          72 : PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
     553             : {
     554          72 :   PetscInt       k;
     555             : 
     556          72 :   PetscFunctionBegin;
     557          72 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     558         288 :   PetscValidLogicalCollectiveInt(eps,i,2);
     559          72 :   if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Wr,3); }
     560          72 :   if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Wi,4); }
     561          72 :   EPSCheckSolved(eps,1);
     562          72 :   PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
     563          72 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     564          72 :   PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
     565          72 :   PetscCall(EPSComputeVectors(eps));
     566          72 :   k = eps->perm[i];
     567          72 :   PetscCall(BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi));
     568          72 :   PetscFunctionReturn(PETSC_SUCCESS);
     569             : }
     570             : 
     571             : /*@
     572             :    EPSGetErrorEstimate - Returns the error estimate associated to the i-th
     573             :    computed eigenpair.
     574             : 
     575             :    Not Collective
     576             : 
     577             :    Input Parameters:
     578             : +  eps - eigensolver context
     579             : -  i   - index of eigenpair
     580             : 
     581             :    Output Parameter:
     582             : .  errest - the error estimate
     583             : 
     584             :    Notes:
     585             :    This is the error estimate used internally by the eigensolver. The actual
     586             :    error bound can be computed with EPSComputeError(). See also the users
     587             :    manual for details.
     588             : 
     589             :    Level: advanced
     590             : 
     591             : .seealso: EPSComputeError()
     592             : @*/
     593          46 : PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
     594             : {
     595          46 :   PetscFunctionBegin;
     596          46 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     597          46 :   PetscAssertPointer(errest,3);
     598          46 :   EPSCheckSolved(eps,1);
     599          46 :   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
     600          46 :   PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
     601          46 :   *errest = eps->errest[eps->perm[i]];
     602          46 :   PetscFunctionReturn(PETSC_SUCCESS);
     603             : }
     604             : 
     605             : /*
     606             :    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
     607             :    associated with an eigenpair.
     608             : 
     609             :    Input Parameters:
     610             :      trans - whether A' must be used instead of A
     611             :      kr,ki - eigenvalue
     612             :      xr,xi - eigenvector
     613             :      z     - three work vectors (the second one not referenced in complex scalars)
     614             : */
     615        3960 : PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
     616             : {
     617        3960 :   PetscInt       nmat;
     618        3960 :   Mat            A,B;
     619        3960 :   Vec            u,w;
     620        3960 :   PetscScalar    alpha;
     621             : #if !defined(PETSC_USE_COMPLEX)
     622        3960 :   Vec            v;
     623        3960 :   PetscReal      ni,nr;
     624             : #endif
     625        3960 :   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
     626             : 
     627        3960 :   PetscFunctionBegin;
     628        3960 :   u = z[0]; w = z[2];
     629        3960 :   PetscCall(STGetNumMatrices(eps->st,&nmat));
     630        3960 :   PetscCall(STGetMatrix(eps->st,0,&A));
     631        3960 :   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
     632             : 
     633             : #if !defined(PETSC_USE_COMPLEX)
     634        3960 :   v = z[1];
     635        3960 :   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
     636             : #endif
     637        3830 :     PetscCall((*matmult)(A,xr,u));                          /* u=A*x */
     638        3830 :     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
     639        3830 :       if (nmat>1) PetscCall((*matmult)(B,xr,w));
     640        1725 :       else PetscCall(VecCopy(xr,w));                        /* w=B*x */
     641        3830 :       alpha = trans? -PetscConj(kr): -kr;
     642        3830 :       PetscCall(VecAXPY(u,alpha,w));                        /* u=A*x-k*B*x */
     643             :     }
     644        3830 :     PetscCall(VecNorm(u,NORM_2,norm));
     645             : #if !defined(PETSC_USE_COMPLEX)
     646             :   } else {
     647         130 :     PetscCall((*matmult)(A,xr,u));                          /* u=A*xr */
     648         130 :     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
     649         130 :       if (nmat>1) PetscCall((*matmult)(B,xr,v));
     650         126 :       else PetscCall(VecCopy(xr,v));                        /* v=B*xr */
     651         130 :       PetscCall(VecAXPY(u,-kr,v));                          /* u=A*xr-kr*B*xr */
     652         130 :       if (nmat>1) PetscCall((*matmult)(B,xi,w));
     653         126 :       else PetscCall(VecCopy(xi,w));                        /* w=B*xi */
     654         130 :       PetscCall(VecAXPY(u,trans?-ki:ki,w));                 /* u=A*xr-kr*B*xr+ki*B*xi */
     655             :     }
     656         130 :     PetscCall(VecNorm(u,NORM_2,&nr));
     657         130 :     PetscCall((*matmult)(A,xi,u));                          /* u=A*xi */
     658         130 :     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
     659         130 :       PetscCall(VecAXPY(u,-kr,w));                          /* u=A*xi-kr*B*xi */
     660         130 :       PetscCall(VecAXPY(u,trans?ki:-ki,v));                 /* u=A*xi-kr*B*xi-ki*B*xr */
     661             :     }
     662         130 :     PetscCall(VecNorm(u,NORM_2,&ni));
     663         130 :     *norm = SlepcAbsEigenvalue(nr,ni);
     664             :   }
     665             : #endif
     666        3960 :   PetscFunctionReturn(PETSC_SUCCESS);
     667             : }
     668             : 
     669             : /*@
     670             :    EPSComputeError - Computes the error (based on the residual norm) associated
     671             :    with the i-th computed eigenpair.
     672             : 
     673             :    Collective
     674             : 
     675             :    Input Parameters:
     676             : +  eps  - the eigensolver context
     677             : .  i    - the solution index
     678             : -  type - the type of error to compute
     679             : 
     680             :    Output Parameter:
     681             : .  error - the error
     682             : 
     683             :    Notes:
     684             :    The error can be computed in various ways, all of them based on the residual
     685             :    norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
     686             : 
     687             :    Level: beginner
     688             : 
     689             : .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
     690             : @*/
     691        3663 : PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
     692             : {
     693        3663 :   Mat            A,B;
     694        3663 :   Vec            xr,xi,w[3];
     695        3663 :   PetscReal      t,vecnorm=1.0,errorl;
     696        3663 :   PetscScalar    kr,ki;
     697        3663 :   PetscBool      flg;
     698             : 
     699        3663 :   PetscFunctionBegin;
     700        3663 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     701       14652 :   PetscValidLogicalCollectiveInt(eps,i,2);
     702       14652 :   PetscValidLogicalCollectiveEnum(eps,type,3);
     703        3663 :   PetscAssertPointer(error,4);
     704        3663 :   EPSCheckSolved(eps,1);
     705             : 
     706             :   /* allocate work vectors */
     707             : #if defined(PETSC_USE_COMPLEX)
     708             :   PetscCall(EPSSetWorkVecs(eps,3));
     709             :   xi   = NULL;
     710             :   w[1] = NULL;
     711             : #else
     712        3663 :   PetscCall(EPSSetWorkVecs(eps,5));
     713        3663 :   xi   = eps->work[3];
     714        3663 :   w[1] = eps->work[4];
     715             : #endif
     716        3663 :   xr   = eps->work[0];
     717        3663 :   w[0] = eps->work[1];
     718        3663 :   w[2] = eps->work[2];
     719             : 
     720             :   /* compute residual norm */
     721        3663 :   PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
     722        3663 :   PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
     723             : 
     724             :   /* compute 2-norm of eigenvector */
     725        3663 :   if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
     726             : 
     727             :   /* if two-sided, compute left residual norm and take the maximum */
     728        3663 :   if (eps->twosided) {
     729          36 :     PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
     730          36 :     PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
     731          51 :     *error = PetscMax(*error,errorl);
     732             :   }
     733             : 
     734             :   /* compute error */
     735        3663 :   switch (type) {
     736             :     case EPS_ERROR_ABSOLUTE:
     737             :       break;
     738        1991 :     case EPS_ERROR_RELATIVE:
     739        1991 :       *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
     740        1991 :       break;
     741        1638 :     case EPS_ERROR_BACKWARD:
     742             :       /* initialization of matrix norms */
     743        1638 :       if (!eps->nrma) {
     744          34 :         PetscCall(STGetMatrix(eps->st,0,&A));
     745          34 :         PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
     746          34 :         PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
     747          34 :         PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
     748             :       }
     749        1638 :       if (eps->isgeneralized) {
     750        1632 :         if (!eps->nrmb) {
     751          33 :           PetscCall(STGetMatrix(eps->st,1,&B));
     752          33 :           PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
     753          33 :           PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
     754          33 :           PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
     755             :         }
     756           6 :       } else eps->nrmb = 1.0;
     757        1638 :       t = SlepcAbsEigenvalue(kr,ki);
     758        1638 :       *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
     759        1638 :       break;
     760           0 :     default:
     761           0 :       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
     762             :   }
     763        3663 :   PetscFunctionReturn(PETSC_SUCCESS);
     764             : }
     765             : 
     766             : /*
     767             :    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
     768             :    for the recurrence that builds the right subspace.
     769             : 
     770             :    Collective
     771             : 
     772             :    Input Parameters:
     773             : +  eps - the eigensolver context
     774             : -  i   - iteration number
     775             : 
     776             :    Output Parameters:
     777             : .  breakdown - flag indicating that a breakdown has occurred
     778             : 
     779             :    Notes:
     780             :    The start vector is computed from another vector: for the first step (i=0),
     781             :    the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
     782             :    vector is created. Then this vector is forced to be in the range of OP (only
     783             :    for generalized definite problems) and orthonormalized with respect to all
     784             :    V-vectors up to i-1. The resulting vector is placed in V[i].
     785             : 
     786             :    The flag breakdown is set to true if either i=0 and the vector belongs to the
     787             :    deflation space, or i>0 and the vector is linearly dependent with respect
     788             :    to the V-vectors.
     789             : */
     790         636 : PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
     791             : {
     792         636 :   PetscReal      norm;
     793         636 :   PetscBool      lindep;
     794         636 :   Vec            w,z;
     795             : 
     796         636 :   PetscFunctionBegin;
     797         636 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     798        2544 :   PetscValidLogicalCollectiveInt(eps,i,2);
     799             : 
     800             :   /* For the first step, use the first initial vector, otherwise a random one */
     801         636 :   if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
     802             : 
     803             :   /* Force the vector to be in the range of OP for definite generalized problems */
     804         636 :   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
     805         112 :     PetscCall(BVCreateVec(eps->V,&w));
     806         112 :     PetscCall(BVCopyVec(eps->V,i,w));
     807         112 :     PetscCall(BVGetColumn(eps->V,i,&z));
     808         112 :     PetscCall(STApply(eps->st,w,z));
     809         112 :     PetscCall(BVRestoreColumn(eps->V,i,&z));
     810         112 :     PetscCall(VecDestroy(&w));
     811             :   }
     812             : 
     813             :   /* Orthonormalize the vector with respect to previous vectors */
     814         636 :   PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
     815         636 :   if (breakdown) *breakdown = lindep;
     816         599 :   else if (lindep || norm == 0.0) {
     817           0 :     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
     818           0 :     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
     819             :   }
     820         636 :   PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
     821         636 :   PetscFunctionReturn(PETSC_SUCCESS);
     822             : }
     823             : 
     824             : /*
     825             :    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
     826             :    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
     827             : */
     828          24 : PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
     829             : {
     830          24 :   PetscReal      norm;
     831          24 :   PetscBool      lindep;
     832             : 
     833          24 :   PetscFunctionBegin;
     834          24 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     835          96 :   PetscValidLogicalCollectiveInt(eps,i,2);
     836             : 
     837             :   /* For the first step, use the first initial vector, otherwise a random one */
     838          24 :   if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
     839             : 
     840             :   /* Orthonormalize the vector with respect to previous vectors */
     841          24 :   PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
     842          24 :   if (breakdown) *breakdown = lindep;
     843          18 :   else if (lindep || norm == 0.0) {
     844           0 :     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
     845           0 :     SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
     846             :   }
     847          24 :   PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
     848          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     849             : }

Generated by: LCOV version 1.14