LCOV - code coverage report
Current view: top level - eps/impls/krylov/lanczos - lanczos.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 451 482 93.6 %
Date: 2024-03-28 00:52:16 Functions: 17 18 94.4 %
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             :    SLEPc eigensolver: "lanczos"
      12             : 
      13             :    Method: Explicitly Restarted Symmetric/Hermitian Lanczos
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Lanczos method for symmetric (Hermitian) problems, with explicit
      18             :        restart and deflation. Several reorthogonalization strategies can
      19             :        be selected.
      20             : 
      21             :    References:
      22             : 
      23             :        [1] "Lanczos Methods in SLEPc", SLEPc Technical Report STR-5,
      24             :            available at https://slepc.upv.es.
      25             : */
      26             : 
      27             : #include <slepc/private/epsimpl.h>                /*I "slepceps.h" I*/
      28             : #include <slepcblaslapack.h>
      29             : 
      30             : typedef struct {
      31             :   EPSLanczosReorthogType reorthog;      /* user-provided reorthogonalization parameter */
      32             :   PetscInt               allocsize;     /* number of columns of work BV's allocated at setup */
      33             :   BV                     AV;            /* work BV used in selective reorthogonalization */
      34             : } EPS_LANCZOS;
      35             : 
      36          36 : static PetscErrorCode EPSSetUp_Lanczos(EPS eps)
      37             : {
      38          36 :   EPS_LANCZOS        *lanczos = (EPS_LANCZOS*)eps->data;
      39          36 :   BVOrthogRefineType refine;
      40          36 :   BVOrthogBlockType  btype;
      41          36 :   PetscReal          eta;
      42             : 
      43          36 :   PetscFunctionBegin;
      44          36 :   EPSCheckHermitianDefinite(eps);
      45          36 :   PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
      46          36 :   PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
      47          36 :   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
      48          36 :   if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
      49          36 :   PetscCheck(eps->which!=EPS_ALL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
      50          36 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
      51          36 :   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
      52             : 
      53          36 :   PetscCheck(lanczos->reorthog!=(EPSLanczosReorthogType)-1,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"You should explicitly provide the reorthogonalization type, e.g., -eps_lanczos_reorthog local. Note that the EPSLANCZOS solver is *NOT RECOMMENDED* for general use, because it uses explicit restart which typically has slow convergence. The recommended solver is EPSKRYLOVSCHUR (the default), which implements Lanczos with thick restart in the case of symmetric/Hermitian problems");
      54             : 
      55          36 :   PetscCall(EPSAllocateSolution(eps,1));
      56          36 :   PetscCall(EPS_SetInnerProduct(eps));
      57          36 :   if (lanczos->reorthog != EPS_LANCZOS_REORTHOG_FULL) {
      58          29 :     PetscCall(BVGetOrthogonalization(eps->V,NULL,&refine,&eta,&btype));
      59          29 :     PetscCall(BVSetOrthogonalization(eps->V,BV_ORTHOG_MGS,refine,eta,btype));
      60          29 :     PetscCall(PetscInfo(eps,"Switching to MGS orthogonalization\n"));
      61             :   }
      62          36 :   if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_SELECTIVE) {
      63           3 :     if (!lanczos->allocsize) {
      64           1 :       PetscCall(BVDuplicate(eps->V,&lanczos->AV));
      65           1 :       PetscCall(BVGetSizes(lanczos->AV,NULL,NULL,&lanczos->allocsize));
      66             :     } else { /* make sure V and AV have the same size */
      67           2 :       PetscCall(BVGetSizes(eps->V,NULL,NULL,&lanczos->allocsize));
      68           2 :       PetscCall(BVResize(lanczos->AV,lanczos->allocsize,PETSC_FALSE));
      69             :     }
      70             :   }
      71             : 
      72          36 :   PetscCall(DSSetType(eps->ds,DSHEP));
      73          36 :   PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
      74          36 :   PetscCall(DSAllocate(eps->ds,eps->ncv+1));
      75          36 :   if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) PetscCall(EPSSetWorkVecs(eps,1));
      76          36 :   PetscFunctionReturn(PETSC_SUCCESS);
      77             : }
      78             : 
      79             : /*
      80             :    EPSLocalLanczos - Local reorthogonalization.
      81             : 
      82             :    This is the simplest variant. At each Lanczos step, the corresponding Lanczos vector
      83             :    is orthogonalized with respect to the two previous Lanczos vectors, according to
      84             :    the three term Lanczos recurrence. WARNING: This variant does not track the loss of
      85             :    orthogonality that occurs in finite-precision arithmetic and, therefore, the
      86             :    generated vectors are not guaranteed to be (semi-)orthogonal.
      87             : */
      88         520 : static PetscErrorCode EPSLocalLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
      89             : {
      90         520 :   PetscInt       i,j,m = *M;
      91         520 :   Mat            Op;
      92         520 :   PetscBool      *which,lwhich[100];
      93         520 :   PetscScalar    *hwork,lhwork[100];
      94             : 
      95         520 :   PetscFunctionBegin;
      96         520 :   if (m > 100) PetscCall(PetscMalloc2(m,&which,m,&hwork));
      97             :   else {
      98         520 :     which = lwhich;
      99         520 :     hwork = lhwork;
     100             :   }
     101        1857 :   for (i=0;i<k;i++) which[i] = PETSC_TRUE;
     102             : 
     103         520 :   PetscCall(BVSetActiveColumns(eps->V,0,m));
     104         520 :   PetscCall(STGetOperator(eps->st,&Op));
     105        7578 :   for (j=k;j<m;j++) {
     106        7058 :     PetscCall(BVMatMultColumn(eps->V,Op,j));
     107        7058 :     which[j] = PETSC_TRUE;
     108        7058 :     if (j-2>=k) which[j-2] = PETSC_FALSE;
     109        7058 :     PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,beta+j,breakdown));
     110        7058 :     alpha[j] = PetscRealPart(hwork[j]);
     111        7058 :     if (PetscUnlikely(*breakdown)) {
     112           0 :       *M = j+1;
     113           0 :       break;
     114        7058 :     } else PetscCall(BVScaleColumn(eps->V,j+1,1/beta[j]));
     115             :   }
     116         520 :   PetscCall(STRestoreOperator(eps->st,&Op));
     117         520 :   if (m > 100) PetscCall(PetscFree2(which,hwork));
     118         520 :   PetscFunctionReturn(PETSC_SUCCESS);
     119             : }
     120             : 
     121             : /*
     122             :    DenseTridiagonal - Solves a real tridiagonal Hermitian Eigenvalue Problem.
     123             : 
     124             :    Input Parameters:
     125             : +  n   - dimension of the eigenproblem
     126             : .  D   - pointer to the array containing the diagonal elements
     127             : -  E   - pointer to the array containing the off-diagonal elements
     128             : 
     129             :    Output Parameters:
     130             : +  w  - pointer to the array to store the computed eigenvalues
     131             : -  V  - pointer to the array to store the eigenvectors
     132             : 
     133             :    Notes:
     134             :    If V is NULL then the eigenvectors are not computed.
     135             : 
     136             :    This routine use LAPACK routines xSTEVR.
     137             : */
     138         673 : static PetscErrorCode DenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
     139             : {
     140         673 :   PetscReal      abstol = 0.0,vl,vu,*work;
     141         673 :   PetscBLASInt   il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
     142         673 :   const char     *jobz;
     143             : #if defined(PETSC_USE_COMPLEX)
     144             :   PetscInt       i,j;
     145             :   PetscReal      *VV=NULL;
     146             : #endif
     147             : 
     148         673 :   PetscFunctionBegin;
     149         673 :   PetscCall(PetscBLASIntCast(n_,&n));
     150         673 :   PetscCall(PetscBLASIntCast(20*n_,&lwork));
     151         673 :   PetscCall(PetscBLASIntCast(10*n_,&liwork));
     152         673 :   if (V) {
     153             :     jobz = "V";
     154             : #if defined(PETSC_USE_COMPLEX)
     155             :     PetscCall(PetscMalloc1(n*n,&VV));
     156             : #endif
     157           0 :   } else jobz = "N";
     158         673 :   PetscCall(PetscMalloc3(2*n,&isuppz,lwork,&work,liwork,&iwork));
     159         673 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     160             : #if defined(PETSC_USE_COMPLEX)
     161             :   PetscCallBLAS("LAPACKstevr",LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,VV,&n,isuppz,work,&lwork,iwork,&liwork,&info));
     162             : #else
     163         673 :   PetscCallBLAS("LAPACKstevr",LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info));
     164             : #endif
     165         673 :   PetscCall(PetscFPTrapPop());
     166         673 :   SlepcCheckLapackInfo("stevr",info);
     167             : #if defined(PETSC_USE_COMPLEX)
     168             :   if (V) {
     169             :     for (i=0;i<n;i++)
     170             :       for (j=0;j<n;j++)
     171             :         V[i*n+j] = VV[i*n+j];
     172             :     PetscCall(PetscFree(VV));
     173             :   }
     174             : #endif
     175         673 :   PetscCall(PetscFree3(isuppz,work,iwork));
     176         673 :   PetscFunctionReturn(PETSC_SUCCESS);
     177             : }
     178             : 
     179             : /*
     180             :    EPSSelectiveLanczos - Selective reorthogonalization.
     181             : */
     182          38 : static PetscErrorCode EPSSelectiveLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscReal anorm)
     183             : {
     184          38 :   EPS_LANCZOS    *lanczos = (EPS_LANCZOS*)eps->data;
     185          38 :   PetscInt       i,j,m = *M,n,nritz=0,nritzo;
     186          38 :   Vec            vj1,av;
     187          38 :   Mat            Op;
     188          38 :   PetscReal      *d,*e,*ritz,norm;
     189          38 :   PetscScalar    *Y,*hwork;
     190          38 :   PetscBool      *which;
     191             : 
     192          38 :   PetscFunctionBegin;
     193          38 :   PetscCall(PetscCalloc6(m+1,&d,m,&e,m,&ritz,m*m,&Y,m,&which,m,&hwork));
     194          87 :   for (i=0;i<k;i++) which[i] = PETSC_TRUE;
     195          38 :   PetscCall(STGetOperator(eps->st,&Op));
     196             : 
     197         711 :   for (j=k;j<m;j++) {
     198         673 :     PetscCall(BVSetActiveColumns(eps->V,0,m));
     199             : 
     200             :     /* Lanczos step */
     201         673 :     PetscCall(BVMatMultColumn(eps->V,Op,j));
     202         673 :     which[j] = PETSC_TRUE;
     203         673 :     if (j-2>=k) which[j-2] = PETSC_FALSE;
     204         673 :     PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
     205         673 :     alpha[j] = PetscRealPart(hwork[j]);
     206         673 :     beta[j] = norm;
     207         673 :     if (PetscUnlikely(*breakdown)) {
     208           0 :       *M = j+1;
     209           0 :       break;
     210             :     }
     211             : 
     212             :     /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
     213         673 :     n = j-k+1;
     214        6994 :     for (i=0;i<n;i++) {
     215        6321 :       d[i] = alpha[i+k];
     216        6321 :       e[i] = beta[i+k];
     217             :     }
     218         673 :     PetscCall(DenseTridiagonal(n,d,e,ritz,Y));
     219             : 
     220             :     /* Estimate ||A|| */
     221        6994 :     for (i=0;i<n;i++)
     222        6321 :       if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
     223             : 
     224             :     /* Compute nearly converged Ritz vectors */
     225             :     nritzo = 0;
     226        6994 :     for (i=0;i<n;i++) {
     227        6321 :       if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) nritzo++;
     228             :     }
     229         673 :     if (nritzo>nritz) {
     230             :       nritz = 0;
     231         196 :       for (i=0;i<n;i++) {
     232         173 :         if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
     233          37 :           PetscCall(BVSetActiveColumns(eps->V,k,k+n));
     234          37 :           PetscCall(BVGetColumn(lanczos->AV,nritz,&av));
     235          37 :           PetscCall(BVMultVec(eps->V,1.0,0.0,av,Y+i*n));
     236          37 :           PetscCall(BVRestoreColumn(lanczos->AV,nritz,&av));
     237          37 :           nritz++;
     238             :         }
     239             :       }
     240             :     }
     241         673 :     if (nritz > 0) {
     242         247 :       PetscCall(BVGetColumn(eps->V,j+1,&vj1));
     243         247 :       PetscCall(BVSetActiveColumns(lanczos->AV,0,nritz));
     244         247 :       PetscCall(BVOrthogonalizeVec(lanczos->AV,vj1,hwork,&norm,breakdown));
     245         247 :       PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
     246         247 :       if (PetscUnlikely(*breakdown)) {
     247           0 :         *M = j+1;
     248           0 :         break;
     249             :       }
     250             :     }
     251         673 :     PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
     252             :   }
     253             : 
     254          38 :   PetscCall(STRestoreOperator(eps->st,&Op));
     255          38 :   PetscCall(PetscFree6(d,e,ritz,Y,which,hwork));
     256          38 :   PetscFunctionReturn(PETSC_SUCCESS);
     257             : }
     258             : 
     259         768 : static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
     260             : {
     261         768 :   PetscInt  k;
     262         768 :   PetscReal T,binv;
     263             : 
     264         768 :   PetscFunctionBegin;
     265             :   /* Estimate of contribution to roundoff errors from A*v
     266             :        fl(A*v) = A*v + f,
     267             :      where ||f|| \approx eps1*||A||.
     268             :      For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps */
     269         768 :   T = eps1*anorm;
     270         768 :   binv = 1.0/beta[j+1];
     271             : 
     272             :   /* Update omega(1) using omega(0)==0 */
     273         768 :   omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] - beta[j]*omega_old[0];
     274         768 :   if (omega_old[0] > 0) omega_old[0] = binv*(omega_old[0] + T);
     275         461 :   else omega_old[0] = binv*(omega_old[0] - T);
     276             : 
     277             :   /* Update remaining components */
     278        7014 :   for (k=1;k<j-1;k++) {
     279        6246 :     omega_old[k] = beta[k+1]*omega[k+1] + (alpha[k]-alpha[j])*omega[k] + beta[k]*omega[k-1] - beta[j]*omega_old[k];
     280        6246 :     if (omega_old[k] > 0) omega_old[k] = binv*(omega_old[k] + T);
     281        3192 :     else omega_old[k] = binv*(omega_old[k] - T);
     282             :   }
     283         768 :   omega_old[j-1] = binv*T;
     284             : 
     285             :   /* Swap omega and omega_old */
     286        8534 :   for (k=0;k<j;k++) {
     287        7766 :     omega[k] = omega_old[k];
     288        7766 :     omega_old[k] = omega[k];
     289             :   }
     290         768 :   omega[j] = eps1;
     291         768 :   PetscFunctionReturnVoid();
     292             : }
     293             : 
     294          48 : static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
     295             : {
     296          48 :   PetscInt  i,k,maxpos;
     297          48 :   PetscReal max;
     298          48 :   PetscBool found;
     299             : 
     300          48 :   PetscFunctionBegin;
     301             :   /* initialize which */
     302          48 :   found = PETSC_FALSE;
     303          48 :   maxpos = 0;
     304          48 :   max = 0.0;
     305         755 :   for (i=0;i<j;i++) {
     306         707 :     if (PetscAbsReal(mu[i]) >= delta) {
     307         175 :       which[i] = PETSC_TRUE;
     308         175 :       found = PETSC_TRUE;
     309         532 :     } else which[i] = PETSC_FALSE;
     310         707 :     if (PetscAbsReal(mu[i]) > max) {
     311          78 :       maxpos = i;
     312          78 :       max = PetscAbsReal(mu[i]);
     313             :     }
     314             :   }
     315          48 :   if (!found) which[maxpos] = PETSC_TRUE;
     316             : 
     317         755 :   for (i=0;i<j;i++) {
     318         707 :     if (which[i]) {
     319             :       /* find left interval */
     320         623 :       for (k=i;k>=0;k--) {
     321         623 :         if (PetscAbsReal(mu[k])<eta || which[k]) break;
     322           0 :         else which[k] = PETSC_TRUE;
     323             :       }
     324             :       /* find right interval */
     325        1071 :       for (k=i+1;k<j;k++) {
     326        1037 :         if (PetscAbsReal(mu[k])<eta || which[k]) break;
     327         448 :         else which[k] = PETSC_TRUE;
     328             :       }
     329             :     }
     330             :   }
     331          48 :   PetscFunctionReturnVoid();
     332             : }
     333             : 
     334             : /*
     335             :    EPSPartialLanczos - Partial reorthogonalization.
     336             : */
     337          76 : static PetscErrorCode EPSPartialLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscReal anorm)
     338             : {
     339          76 :   EPS_LANCZOS    *lanczos = (EPS_LANCZOS*)eps->data;
     340          76 :   PetscInt       i,j,m = *M;
     341          76 :   Mat            Op;
     342          76 :   PetscReal      norm,*omega,lomega[100],*omega_old,lomega_old[100],eps1,delta,eta;
     343          76 :   PetscBool      *which,lwhich[100],*which2,lwhich2[100];
     344          76 :   PetscBool      reorth = PETSC_FALSE,force_reorth = PETSC_FALSE;
     345          76 :   PetscBool      fro = PETSC_FALSE,estimate_anorm = PETSC_FALSE;
     346          76 :   PetscScalar    *hwork,lhwork[100];
     347             : 
     348          76 :   PetscFunctionBegin;
     349          76 :   if (m>100) PetscCall(PetscMalloc5(m,&omega,m,&omega_old,m,&which,m,&which2,m,&hwork));
     350             :   else {
     351          76 :     omega     = lomega;
     352          76 :     omega_old = lomega_old;
     353          76 :     which     = lwhich;
     354          76 :     which2    = lwhich2;
     355          76 :     hwork     = lhwork;
     356             :   }
     357             : 
     358          76 :   eps1 = PetscSqrtReal((PetscReal)eps->n)*PETSC_MACHINE_EPSILON/2;
     359          76 :   delta = PETSC_SQRT_MACHINE_EPSILON/PetscSqrtReal((PetscReal)eps->ncv);
     360          76 :   eta = PetscPowReal(PETSC_MACHINE_EPSILON,3.0/4.0)/PetscSqrtReal((PetscReal)eps->ncv);
     361          76 :   if (anorm < 0.0) {
     362           6 :     anorm = 1.0;
     363           6 :     estimate_anorm = PETSC_TRUE;
     364             :   }
     365        7676 :   for (i=0;i<PetscMax(100,m);i++) omega[i] = omega_old[i] = 0.0;
     366         174 :   for (i=0;i<k;i++) which[i] = PETSC_TRUE;
     367             : 
     368          76 :   PetscCall(BVSetActiveColumns(eps->V,0,m));
     369          76 :   PetscCall(STGetOperator(eps->st,&Op));
     370        1422 :   for (j=k;j<m;j++) {
     371        1346 :     PetscCall(BVMatMultColumn(eps->V,Op,j));
     372        1346 :     if (fro) {
     373             :       /* Lanczos step with full reorthogonalization */
     374         502 :       PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
     375         502 :       alpha[j] = PetscRealPart(hwork[j]);
     376             :     } else {
     377             :       /* Lanczos step */
     378         844 :       which[j] = PETSC_TRUE;
     379         844 :       if (j-2>=k) which[j-2] = PETSC_FALSE;
     380         844 :       PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
     381         844 :       alpha[j] = PetscRealPart(hwork[j]);
     382         844 :       beta[j] = norm;
     383             : 
     384             :       /* Estimate ||A|| if needed */
     385         844 :       if (estimate_anorm) {
     386         114 :         if (j>k) anorm = PetscMax(anorm,PetscAbsReal(alpha[j])+norm+beta[j-1]);
     387           6 :         else anorm = PetscMax(anorm,PetscAbsReal(alpha[j])+norm);
     388             :       }
     389             : 
     390             :       /* Check if reorthogonalization is needed */
     391         844 :       reorth = PETSC_FALSE;
     392         844 :       if (j>k) {
     393         768 :         update_omega(omega,omega_old,j,alpha,beta-1,eps1,anorm);
     394        8362 :         for (i=0;i<j-k;i++) {
     395        6826 :           if (PetscAbsReal(omega[i]) > delta) reorth = PETSC_TRUE;
     396             :         }
     397             :       }
     398         844 :       if (reorth || force_reorth) {
     399         290 :         for (i=0;i<k;i++) which2[i] = PETSC_FALSE;
     400        2653 :         for (i=k;i<=j;i++) which2[i] = PETSC_TRUE;
     401         159 :         if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
     402             :           /* Periodic reorthogonalization */
     403          79 :           if (force_reorth) force_reorth = PETSC_FALSE;
     404          47 :           else force_reorth = PETSC_TRUE;
     405        1238 :           for (i=0;i<j-k;i++) omega[i] = eps1;
     406             :         } else {
     407             :           /* Partial reorthogonalization */
     408          80 :           if (force_reorth) force_reorth = PETSC_FALSE;
     409             :           else {
     410          48 :             force_reorth = PETSC_TRUE;
     411          48 :             compute_int(which2+k,omega,j-k,delta,eta);
     412         803 :             for (i=0;i<j-k;i++) {
     413         707 :               if (which2[i+k]) omega[i] = eps1;
     414             :             }
     415             :           }
     416             :         }
     417         159 :         PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which2,hwork,&norm,breakdown));
     418             :       }
     419             :     }
     420             : 
     421        1346 :     if (PetscUnlikely(*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON)) {
     422           0 :       *M = j+1;
     423           0 :       break;
     424             :     }
     425        1346 :     if (!fro && norm*delta < anorm*eps1) {
     426          30 :       fro = PETSC_TRUE;
     427          30 :       PetscCall(PetscInfo(eps,"Switching to full reorthogonalization at iteration %" PetscInt_FMT "\n",eps->its));
     428             :     }
     429        1346 :     beta[j] = norm;
     430        1346 :     PetscCall(BVScaleColumn(eps->V,j+1,1.0/norm));
     431             :   }
     432             : 
     433          76 :   PetscCall(STRestoreOperator(eps->st,&Op));
     434          76 :   if (m>100) PetscCall(PetscFree5(omega,omega_old,which,which2,hwork));
     435          76 :   PetscFunctionReturn(PETSC_SUCCESS);
     436             : }
     437             : 
     438             : /*
     439             :    EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
     440             :    columns are assumed to be locked and therefore they are not modified. On
     441             :    exit, the following relation is satisfied:
     442             : 
     443             :                     OP * V - V * T = f * e_m^T
     444             : 
     445             :    where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
     446             :    f is the residual vector and e_m is the m-th vector of the canonical basis.
     447             :    The Lanczos vectors (together with vector f) are B-orthogonal (to working
     448             :    accuracy) if full reorthogonalization is being used, otherwise they are
     449             :    (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
     450             :    Lanczos vector can be computed as v_{m+1} = f / beta.
     451             : 
     452             :    This function simply calls another function which depends on the selected
     453             :    reorthogonalization strategy.
     454             : */
     455         932 : static PetscErrorCode EPSBasicLanczos(EPS eps,PetscInt k,PetscInt *m,PetscReal *betam,PetscBool *breakdown,PetscReal anorm)
     456             : {
     457         932 :   EPS_LANCZOS        *lanczos = (EPS_LANCZOS*)eps->data;
     458         932 :   PetscScalar        *T;
     459         932 :   PetscInt           i,n=*m,ld;
     460         932 :   PetscReal          *alpha,*beta;
     461         932 :   BVOrthogRefineType orthog_ref;
     462         932 :   Mat                Op,M;
     463             : 
     464         932 :   PetscFunctionBegin;
     465         932 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     466         932 :   switch (lanczos->reorthog) {
     467         520 :     case EPS_LANCZOS_REORTHOG_LOCAL:
     468         520 :       PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     469         520 :       beta = alpha + ld;
     470         520 :       PetscCall(EPSLocalLanczos(eps,alpha,beta,k,m,breakdown));
     471         520 :       *betam = beta[*m-1];
     472         520 :       PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     473             :       break;
     474         271 :     case EPS_LANCZOS_REORTHOG_FULL:
     475         271 :       PetscCall(STGetOperator(eps->st,&Op));
     476         271 :       PetscCall(DSGetMat(eps->ds,DS_MAT_T,&M));
     477         271 :       PetscCall(BVMatLanczos(eps->V,Op,M,k,m,betam,breakdown));
     478         271 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&M));
     479         271 :       PetscCall(STRestoreOperator(eps->st,&Op));
     480             :       break;
     481          38 :     case EPS_LANCZOS_REORTHOG_SELECTIVE:
     482          38 :       PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     483          38 :       beta = alpha + ld;
     484          38 :       PetscCall(EPSSelectiveLanczos(eps,alpha,beta,k,m,breakdown,anorm));
     485          38 :       *betam = beta[*m-1];
     486          38 :       PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     487             :       break;
     488          76 :     case EPS_LANCZOS_REORTHOG_PERIODIC:
     489             :     case EPS_LANCZOS_REORTHOG_PARTIAL:
     490          76 :       PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     491          76 :       beta = alpha + ld;
     492          76 :       PetscCall(EPSPartialLanczos(eps,alpha,beta,k,m,breakdown,anorm));
     493          76 :       *betam = beta[*m-1];
     494          76 :       PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     495             :       break;
     496          27 :     case EPS_LANCZOS_REORTHOG_DELAYED:
     497          27 :       PetscCall(PetscMalloc1(n*n,&T));
     498          27 :       PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
     499          27 :       if (orthog_ref == BV_ORTHOG_REFINE_NEVER) PetscCall(EPSDelayedArnoldi1(eps,T,n,k,m,betam,breakdown));
     500          27 :       else PetscCall(EPSDelayedArnoldi(eps,T,n,k,m,betam,breakdown));
     501          27 :       n = *m;
     502          27 :       PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     503          27 :       beta = alpha + ld;
     504         472 :       for (i=k;i<n-1;i++) {
     505         445 :         alpha[i] = PetscRealPart(T[n*i+i]);
     506         445 :         beta[i] = PetscRealPart(T[n*i+i+1]);
     507             :       }
     508          27 :       alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
     509          27 :       beta[n-1] = *betam;
     510          27 :       PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     511          27 :       PetscCall(PetscFree(T));
     512             :       break;
     513             :   }
     514         932 :   PetscFunctionReturn(PETSC_SUCCESS);
     515             : }
     516             : 
     517          36 : static PetscErrorCode EPSSolve_Lanczos(EPS eps)
     518             : {
     519          36 :   EPS_LANCZOS    *lanczos = (EPS_LANCZOS*)eps->data;
     520          36 :   PetscInt       nconv,i,j,k,l,x,n,*perm,restart,ncv=eps->ncv,r,ld;
     521          36 :   Vec            vi,vj,w;
     522          36 :   Mat            U;
     523          36 :   PetscScalar    *Y,*ritz,stmp;
     524          36 :   PetscReal      *bnd,anorm,beta,norm,rtmp,resnorm;
     525          36 :   PetscBool      breakdown;
     526          36 :   char           *conv,ctmp;
     527             : 
     528          36 :   PetscFunctionBegin;
     529          36 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     530          36 :   PetscCall(PetscMalloc4(ncv,&ritz,ncv,&bnd,ncv,&perm,ncv,&conv));
     531             : 
     532             :   /* The first Lanczos vector is the normalized initial vector */
     533          36 :   PetscCall(EPSGetStartVector(eps,0,NULL));
     534             : 
     535             :   anorm = -1.0;
     536             :   nconv = 0;
     537             : 
     538             :   /* Restart loop */
     539         968 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     540         932 :     eps->its++;
     541             : 
     542             :     /* Compute an ncv-step Lanczos factorization */
     543         932 :     n = PetscMin(nconv+eps->mpd,ncv);
     544         932 :     PetscCall(DSSetDimensions(eps->ds,n,nconv,PETSC_DEFAULT));
     545         932 :     PetscCall(EPSBasicLanczos(eps,nconv,&n,&beta,&breakdown,anorm));
     546         932 :     PetscCall(DSSetDimensions(eps->ds,n,nconv,0));
     547         932 :     PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
     548         932 :     PetscCall(BVSetActiveColumns(eps->V,nconv,n));
     549             : 
     550             :     /* Solve projected problem */
     551         932 :     PetscCall(DSSolve(eps->ds,ritz,NULL));
     552         932 :     PetscCall(DSSort(eps->ds,ritz,NULL,NULL,NULL,NULL));
     553         932 :     PetscCall(DSSynchronize(eps->ds,ritz,NULL));
     554             : 
     555             :     /* Estimate ||A|| */
     556       13272 :     for (i=nconv;i<n;i++)
     557       12761 :       anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
     558             : 
     559             :     /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
     560         932 :     PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
     561       13272 :     for (i=nconv;i<n;i++) {
     562       12340 :       resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
     563       12340 :       PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->convergedctx));
     564       12340 :       if (bnd[i]<eps->tol) conv[i] = 'C';
     565       12171 :       else conv[i] = 'N';
     566             :     }
     567         932 :     PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
     568             : 
     569             :     /* purge repeated ritz values */
     570         932 :     if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
     571        7058 :       for (i=nconv+1;i<n;i++) {
     572        6538 :         if (conv[i] == 'C' && PetscAbsScalar((ritz[i]-ritz[i-1])/ritz[i]) < eps->tol) conv[i] = 'R';
     573             :       }
     574             :     }
     575             : 
     576             :     /* Compute restart vector */
     577         932 :     if (breakdown) PetscCall(PetscInfo(eps,"Breakdown in Lanczos method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
     578             :     else {
     579             :       restart = nconv;
     580        1078 :       while (restart<n && conv[restart] != 'N') restart++;
     581         931 :       if (restart >= n) {
     582           0 :         breakdown = PETSC_TRUE;
     583             :       } else {
     584       12175 :         for (i=restart+1;i<n;i++) {
     585       11244 :           if (conv[i] == 'N') {
     586       11240 :             PetscCall(SlepcSCCompare(eps->sc,ritz[restart],0.0,ritz[i],0.0,&r));
     587       11240 :             if (r>0) restart = i;
     588             :           }
     589             :         }
     590         931 :         PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
     591         931 :         PetscCall(BVMultColumn(eps->V,1.0,0.0,n,Y+restart*ld+nconv));
     592         931 :         PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
     593             :       }
     594             :     }
     595             : 
     596             :     /* Count and put converged eigenvalues first */
     597       13272 :     for (i=nconv;i<n;i++) perm[i] = i;
     598        1098 :     for (k=nconv;k<n;k++) {
     599        1097 :       if (conv[perm[k]] != 'C') {
     600         934 :         j = k + 1;
     601       12181 :         while (j<n && conv[perm[j]] != 'C') j++;
     602         934 :         if (j>=n) break;
     603           3 :         l = perm[k]; perm[k] = perm[j]; perm[j] = l;
     604             :       }
     605             :     }
     606             : 
     607             :     /* Sort eigenvectors according to permutation */
     608         932 :     PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
     609        1098 :     for (i=nconv;i<k;i++) {
     610         166 :       x = perm[i];
     611         166 :       if (x != i) {
     612           3 :         j = i + 1;
     613           7 :         while (perm[j] != i) j++;
     614             :         /* swap eigenvalues i and j */
     615           3 :         stmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = stmp;
     616           3 :         rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
     617           3 :         ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
     618           3 :         perm[j] = x; perm[i] = i;
     619             :         /* swap eigenvectors i and j */
     620          60 :         for (l=0;l<n;l++) {
     621          57 :           stmp = Y[l+x*ld]; Y[l+x*ld] = Y[l+i*ld]; Y[l+i*ld] = stmp;
     622             :         }
     623             :       }
     624             :     }
     625         932 :     PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
     626             : 
     627             :     /* compute converged eigenvectors */
     628         932 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
     629         932 :     PetscCall(BVMultInPlace(eps->V,U,nconv,k));
     630         932 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
     631             : 
     632             :     /* purge spurious ritz values */
     633         932 :     if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
     634         595 :       for (i=nconv;i<k;i++) {
     635          75 :         PetscCall(BVGetColumn(eps->V,i,&vi));
     636          75 :         PetscCall(VecNorm(vi,NORM_2,&norm));
     637          75 :         PetscCall(VecScale(vi,1.0/norm));
     638          75 :         w = eps->work[0];
     639          75 :         PetscCall(STApply(eps->st,vi,w));
     640          75 :         PetscCall(VecAXPY(w,-ritz[i],vi));
     641          75 :         PetscCall(BVRestoreColumn(eps->V,i,&vi));
     642          75 :         PetscCall(VecNorm(w,NORM_2,&norm));
     643          75 :         PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],norm,&bnd[i],eps->convergedctx));
     644          75 :         if (bnd[i]>=eps->tol) conv[i] = 'S';
     645             :       }
     646         594 :       for (i=nconv;i<k;i++) {
     647          75 :         if (conv[i] != 'C') {
     648           1 :           j = i + 1;
     649           1 :           while (j<k && conv[j] != 'C') j++;
     650           1 :           if (j>=k) break;
     651             :           /* swap eigenvalues i and j */
     652           0 :           stmp = ritz[j]; ritz[j] = ritz[i]; ritz[i] = stmp;
     653           0 :           rtmp = bnd[j]; bnd[j] = bnd[i]; bnd[i] = rtmp;
     654           0 :           ctmp = conv[j]; conv[j] = conv[i]; conv[i] = ctmp;
     655             :           /* swap eigenvectors i and j */
     656           0 :           PetscCall(BVGetColumn(eps->V,i,&vi));
     657           0 :           PetscCall(BVGetColumn(eps->V,j,&vj));
     658           0 :           PetscCall(VecSwap(vi,vj));
     659           0 :           PetscCall(BVRestoreColumn(eps->V,i,&vi));
     660          74 :           PetscCall(BVRestoreColumn(eps->V,j,&vj));
     661             :         }
     662             :       }
     663             :       k = i;
     664             :     }
     665             : 
     666             :     /* store ritz values and estimated errors */
     667       13272 :     for (i=nconv;i<n;i++) {
     668       12340 :       eps->eigr[i] = ritz[i];
     669       12340 :       eps->errest[i] = bnd[i];
     670             :     }
     671         932 :     nconv = k;
     672         932 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n));
     673         932 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,nconv,eps->nev,&eps->reason,eps->stoppingctx));
     674             : 
     675         932 :     if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
     676         896 :       PetscCall(BVCopyColumn(eps->V,n,nconv));
     677         896 :       if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
     678             :         /* Reorthonormalize restart vector */
     679         502 :         PetscCall(BVOrthonormalizeColumn(eps->V,nconv,PETSC_FALSE,NULL,&breakdown));
     680             :       }
     681         896 :       if (breakdown) {
     682             :         /* Use random vector for restarting */
     683           0 :         PetscCall(PetscInfo(eps,"Using random vector for restart\n"));
     684           0 :         PetscCall(EPSGetStartVector(eps,nconv,&breakdown));
     685             :       }
     686         896 :       if (PetscUnlikely(breakdown)) { /* give up */
     687           0 :         eps->reason = EPS_DIVERGED_BREAKDOWN;
     688         968 :         PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
     689             :       }
     690             :     }
     691             :   }
     692          36 :   eps->nconv = nconv;
     693             : 
     694          36 :   PetscCall(PetscFree4(ritz,bnd,perm,conv));
     695          36 :   PetscFunctionReturn(PETSC_SUCCESS);
     696             : }
     697             : 
     698          18 : static PetscErrorCode EPSSetFromOptions_Lanczos(EPS eps,PetscOptionItems *PetscOptionsObject)
     699             : {
     700          18 :   EPS_LANCZOS            *lanczos = (EPS_LANCZOS*)eps->data;
     701          18 :   PetscBool              flg;
     702          18 :   EPSLanczosReorthogType reorthog=EPS_LANCZOS_REORTHOG_LOCAL,curval;
     703             : 
     704          18 :   PetscFunctionBegin;
     705          18 :   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Lanczos Options");
     706             : 
     707          18 :     curval = (lanczos->reorthog==(EPSLanczosReorthogType)-1)? EPS_LANCZOS_REORTHOG_LOCAL: lanczos->reorthog;
     708          18 :     PetscCall(PetscOptionsEnum("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",EPSLanczosReorthogTypes,(PetscEnum)curval,(PetscEnum*)&reorthog,&flg));
     709          18 :     if (flg) PetscCall(EPSLanczosSetReorthog(eps,reorthog));
     710             : 
     711          18 :   PetscOptionsHeadEnd();
     712          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     713             : }
     714             : 
     715          19 : static PetscErrorCode EPSLanczosSetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType reorthog)
     716             : {
     717          19 :   EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
     718             : 
     719          19 :   PetscFunctionBegin;
     720          19 :   switch (reorthog) {
     721          19 :     case EPS_LANCZOS_REORTHOG_LOCAL:
     722             :     case EPS_LANCZOS_REORTHOG_FULL:
     723             :     case EPS_LANCZOS_REORTHOG_DELAYED:
     724             :     case EPS_LANCZOS_REORTHOG_SELECTIVE:
     725             :     case EPS_LANCZOS_REORTHOG_PERIODIC:
     726             :     case EPS_LANCZOS_REORTHOG_PARTIAL:
     727          19 :       if (lanczos->reorthog != reorthog) {
     728          19 :         lanczos->reorthog = reorthog;
     729          19 :         eps->state = EPS_STATE_INITIAL;
     730             :       }
     731          19 :       break;
     732           0 :     default:
     733           0 :       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
     734             :   }
     735          19 :   PetscFunctionReturn(PETSC_SUCCESS);
     736             : }
     737             : 
     738             : /*@
     739             :    EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
     740             :    iteration.
     741             : 
     742             :    Logically Collective
     743             : 
     744             :    Input Parameters:
     745             : +  eps - the eigenproblem solver context
     746             : -  reorthog - the type of reorthogonalization
     747             : 
     748             :    Options Database Key:
     749             : .  -eps_lanczos_reorthog - Sets the reorthogonalization type (either 'local', 'selective',
     750             :                          'periodic', 'partial', 'full' or 'delayed')
     751             : 
     752             :    Level: advanced
     753             : 
     754             : .seealso: EPSLanczosGetReorthog(), EPSLanczosReorthogType
     755             : @*/
     756          19 : PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
     757             : {
     758          19 :   PetscFunctionBegin;
     759          19 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     760          76 :   PetscValidLogicalCollectiveEnum(eps,reorthog,2);
     761          19 :   PetscTryMethod(eps,"EPSLanczosSetReorthog_C",(EPS,EPSLanczosReorthogType),(eps,reorthog));
     762          19 :   PetscFunctionReturn(PETSC_SUCCESS);
     763             : }
     764             : 
     765           5 : static PetscErrorCode EPSLanczosGetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType *reorthog)
     766             : {
     767           5 :   EPS_LANCZOS *lanczos = (EPS_LANCZOS*)eps->data;
     768             : 
     769           5 :   PetscFunctionBegin;
     770           5 :   *reorthog = lanczos->reorthog;
     771           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     772             : }
     773             : 
     774             : /*@
     775             :    EPSLanczosGetReorthog - Gets the type of reorthogonalization used during
     776             :    the Lanczos iteration.
     777             : 
     778             :    Not Collective
     779             : 
     780             :    Input Parameter:
     781             : .  eps - the eigenproblem solver context
     782             : 
     783             :    Output Parameter:
     784             : .  reorthog - the type of reorthogonalization
     785             : 
     786             :    Level: advanced
     787             : 
     788             : .seealso: EPSLanczosSetReorthog(), EPSLanczosReorthogType
     789             : @*/
     790           5 : PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
     791             : {
     792           5 :   PetscFunctionBegin;
     793           5 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     794           5 :   PetscAssertPointer(reorthog,2);
     795           5 :   PetscUseMethod(eps,"EPSLanczosGetReorthog_C",(EPS,EPSLanczosReorthogType*),(eps,reorthog));
     796           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     797             : }
     798             : 
     799          22 : static PetscErrorCode EPSReset_Lanczos(EPS eps)
     800             : {
     801          22 :   EPS_LANCZOS    *lanczos = (EPS_LANCZOS*)eps->data;
     802             : 
     803          22 :   PetscFunctionBegin;
     804          22 :   PetscCall(BVDestroy(&lanczos->AV));
     805          22 :   lanczos->allocsize = 0;
     806          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     807             : }
     808             : 
     809          19 : static PetscErrorCode EPSDestroy_Lanczos(EPS eps)
     810             : {
     811          19 :   PetscFunctionBegin;
     812          19 :   PetscCall(PetscFree(eps->data));
     813          19 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",NULL));
     814          19 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",NULL));
     815          19 :   PetscFunctionReturn(PETSC_SUCCESS);
     816             : }
     817             : 
     818           0 : static PetscErrorCode EPSView_Lanczos(EPS eps,PetscViewer viewer)
     819             : {
     820           0 :   EPS_LANCZOS    *lanczos = (EPS_LANCZOS*)eps->data;
     821           0 :   PetscBool      isascii;
     822             : 
     823           0 :   PetscFunctionBegin;
     824           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     825           0 :   if (isascii) {
     826           0 :     if (lanczos->reorthog != (EPSLanczosReorthogType)-1) PetscCall(PetscViewerASCIIPrintf(viewer,"  %s reorthogonalization\n",EPSLanczosReorthogTypes[lanczos->reorthog]));
     827             :   }
     828           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     829             : }
     830             : 
     831          19 : SLEPC_EXTERN PetscErrorCode EPSCreate_Lanczos(EPS eps)
     832             : {
     833          19 :   EPS_LANCZOS    *ctx;
     834             : 
     835          19 :   PetscFunctionBegin;
     836          19 :   PetscCall(PetscNew(&ctx));
     837          19 :   eps->data = (void*)ctx;
     838          19 :   ctx->reorthog = (EPSLanczosReorthogType)-1;
     839             : 
     840          19 :   eps->useds = PETSC_TRUE;
     841             : 
     842          19 :   eps->ops->solve          = EPSSolve_Lanczos;
     843          19 :   eps->ops->setup          = EPSSetUp_Lanczos;
     844          19 :   eps->ops->setupsort      = EPSSetUpSort_Default;
     845          19 :   eps->ops->setfromoptions = EPSSetFromOptions_Lanczos;
     846          19 :   eps->ops->destroy        = EPSDestroy_Lanczos;
     847          19 :   eps->ops->reset          = EPSReset_Lanczos;
     848          19 :   eps->ops->view           = EPSView_Lanczos;
     849          19 :   eps->ops->backtransform  = EPSBackTransform_Default;
     850          19 :   eps->ops->computevectors = EPSComputeVectors_Hermitian;
     851             : 
     852          19 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",EPSLanczosSetReorthog_Lanczos));
     853          19 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",EPSLanczosGetReorthog_Lanczos));
     854          19 :   PetscFunctionReturn(PETSC_SUCCESS);
     855             : }

Generated by: LCOV version 1.14