LCOV - code coverage report
Current view: top level - eps/impls/krylov/lanczos - lanczos.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 460 490 93.9 %
Date: 2024-03-28 00:28:38 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         514 : static PetscErrorCode EPSLocalLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
      89             : {
      90         514 :   PetscInt       i,j,m = *M;
      91         514 :   Mat            Op;
      92         514 :   PetscBool      *which,lwhich[100];
      93         514 :   PetscScalar    *hwork,lhwork[100];
      94             : 
      95         514 :   PetscFunctionBegin;
      96         514 :   if (m > 100) PetscCall(PetscMalloc2(m,&which,m,&hwork));
      97             :   else {
      98         514 :     which = lwhich;
      99         514 :     hwork = lhwork;
     100             :   }
     101        1784 :   for (i=0;i<k;i++) which[i] = PETSC_TRUE;
     102             : 
     103         514 :   PetscCall(BVSetActiveColumns(eps->V,0,m));
     104         514 :   PetscCall(STGetOperator(eps->st,&Op));
     105        7480 :   for (j=k;j<m;j++) {
     106        6966 :     PetscCall(BVMatMultColumn(eps->V,Op,j));
     107        6966 :     which[j] = PETSC_TRUE;
     108        6966 :     if (j-2>=k) which[j-2] = PETSC_FALSE;
     109        6966 :     PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,beta+j,breakdown));
     110        6966 :     alpha[j] = PetscRealPart(hwork[j]);
     111        6966 :     if (PetscUnlikely(*breakdown)) {
     112           0 :       *M = j+1;
     113           0 :       break;
     114        6966 :     } else PetscCall(BVScaleColumn(eps->V,j+1,1/beta[j]));
     115             :   }
     116         514 :   PetscCall(STRestoreOperator(eps->st,&Op));
     117         514 :   if (m > 100) PetscCall(PetscFree2(which,hwork));
     118         514 :   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         676 : static PetscErrorCode DenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
     139             : {
     140         676 :   PetscReal      abstol = 0.0,vl,vu,*work;
     141         676 :   PetscBLASInt   il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
     142         676 :   const char     *jobz;
     143             : #if defined(PETSC_USE_COMPLEX)
     144         676 :   PetscInt       i,j;
     145         676 :   PetscReal      *VV=NULL;
     146             : #endif
     147             : 
     148         676 :   PetscFunctionBegin;
     149         676 :   PetscCall(PetscBLASIntCast(n_,&n));
     150         676 :   PetscCall(PetscBLASIntCast(20*n_,&lwork));
     151         676 :   PetscCall(PetscBLASIntCast(10*n_,&liwork));
     152         676 :   if (V) {
     153         676 :     jobz = "V";
     154             : #if defined(PETSC_USE_COMPLEX)
     155         676 :     PetscCall(PetscMalloc1(n*n,&VV));
     156             : #endif
     157             :   } else jobz = "N";
     158         676 :   PetscCall(PetscMalloc3(2*n,&isuppz,lwork,&work,liwork,&iwork));
     159         676 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     160             : #if defined(PETSC_USE_COMPLEX)
     161         676 :   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             :   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         676 :   PetscCall(PetscFPTrapPop());
     166         676 :   SlepcCheckLapackInfo("stevr",info);
     167             : #if defined(PETSC_USE_COMPLEX)
     168         676 :   if (V) {
     169        7053 :     for (i=0;i<n;i++)
     170       84762 :       for (j=0;j<n;j++)
     171       78385 :         V[i*n+j] = VV[i*n+j];
     172         676 :     PetscCall(PetscFree(VV));
     173             :   }
     174             : #endif
     175         676 :   PetscCall(PetscFree3(isuppz,work,iwork));
     176         676 :   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          84 :   for (i=0;i<k;i++) which[i] = PETSC_TRUE;
     195          38 :   PetscCall(STGetOperator(eps->st,&Op));
     196             : 
     197         714 :   for (j=k;j<m;j++) {
     198         676 :     PetscCall(BVSetActiveColumns(eps->V,0,m));
     199             : 
     200             :     /* Lanczos step */
     201         676 :     PetscCall(BVMatMultColumn(eps->V,Op,j));
     202         676 :     which[j] = PETSC_TRUE;
     203         676 :     if (j-2>=k) which[j-2] = PETSC_FALSE;
     204         676 :     PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
     205         676 :     alpha[j] = PetscRealPart(hwork[j]);
     206         676 :     beta[j] = norm;
     207         676 :     if (PetscUnlikely(*breakdown)) {
     208           0 :       *M = j+1;
     209           0 :       break;
     210             :     }
     211             : 
     212             :     /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
     213         676 :     n = j-k+1;
     214        7053 :     for (i=0;i<n;i++) {
     215        6377 :       d[i] = alpha[i+k];
     216        6377 :       e[i] = beta[i+k];
     217             :     }
     218         676 :     PetscCall(DenseTridiagonal(n,d,e,ritz,Y));
     219             : 
     220             :     /* Estimate ||A|| */
     221        7053 :     for (i=0;i<n;i++)
     222        6377 :       if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
     223             : 
     224             :     /* Compute nearly converged Ritz vectors */
     225             :     nritzo = 0;
     226        7053 :     for (i=0;i<n;i++) {
     227        6377 :       if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) nritzo++;
     228             :     }
     229         676 :     if (nritzo>nritz) {
     230             :       nritz = 0;
     231         160 :       for (i=0;i<n;i++) {
     232         138 :         if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
     233          32 :           PetscCall(BVSetActiveColumns(eps->V,k,k+n));
     234          32 :           PetscCall(BVGetColumn(lanczos->AV,nritz,&av));
     235          32 :           PetscCall(BVMultVec(eps->V,1.0,0.0,av,Y+i*n));
     236          32 :           PetscCall(BVRestoreColumn(lanczos->AV,nritz,&av));
     237          32 :           nritz++;
     238             :         }
     239             :       }
     240             :     }
     241         676 :     if (nritz > 0) {
     242         255 :       PetscCall(BVGetColumn(eps->V,j+1,&vj1));
     243         255 :       PetscCall(BVSetActiveColumns(lanczos->AV,0,nritz));
     244         255 :       PetscCall(BVOrthogonalizeVec(lanczos->AV,vj1,hwork,&norm,breakdown));
     245         255 :       PetscCall(BVRestoreColumn(eps->V,j+1,&vj1));
     246         255 :       if (PetscUnlikely(*breakdown)) {
     247           0 :         *M = j+1;
     248           0 :         break;
     249             :       }
     250             :     }
     251         676 :     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         842 : static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
     260             : {
     261         842 :   PetscInt  k;
     262         842 :   PetscReal T,binv;
     263             : 
     264         842 :   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         842 :   T = eps1*anorm;
     270         842 :   binv = 1.0/beta[j+1];
     271             : 
     272             :   /* Update omega(1) using omega(0)==0 */
     273         842 :   omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] - beta[j]*omega_old[0];
     274         842 :   if (omega_old[0] > 0) omega_old[0] = binv*(omega_old[0] + T);
     275         506 :   else omega_old[0] = binv*(omega_old[0] - T);
     276             : 
     277             :   /* Update remaining components */
     278        7632 :   for (k=1;k<j-1;k++) {
     279        6790 :     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        6790 :     if (omega_old[k] > 0) omega_old[k] = binv*(omega_old[k] + T);
     281        3594 :     else omega_old[k] = binv*(omega_old[k] - T);
     282             :   }
     283         842 :   omega_old[j-1] = binv*T;
     284             : 
     285             :   /* Swap omega and omega_old */
     286        9296 :   for (k=0;k<j;k++) {
     287        8454 :     omega[k] = omega_old[k];
     288        8454 :     omega_old[k] = omega[k];
     289             :   }
     290         842 :   omega[j] = eps1;
     291         842 :   PetscFunctionReturnVoid();
     292             : }
     293             : 
     294          53 : static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
     295             : {
     296          53 :   PetscInt  i,k,maxpos;
     297          53 :   PetscReal max;
     298          53 :   PetscBool found;
     299             : 
     300          53 :   PetscFunctionBegin;
     301             :   /* initialize which */
     302          53 :   found = PETSC_FALSE;
     303          53 :   maxpos = 0;
     304          53 :   max = 0.0;
     305         843 :   for (i=0;i<j;i++) {
     306         790 :     if (PetscAbsReal(mu[i]) >= delta) {
     307         189 :       which[i] = PETSC_TRUE;
     308         189 :       found = PETSC_TRUE;
     309         601 :     } else which[i] = PETSC_FALSE;
     310         790 :     if (PetscAbsReal(mu[i]) > max) {
     311          88 :       maxpos = i;
     312          88 :       max = PetscAbsReal(mu[i]);
     313             :     }
     314             :   }
     315          53 :   if (!found) which[maxpos] = PETSC_TRUE;
     316             : 
     317         843 :   for (i=0;i<j;i++) {
     318         790 :     if (which[i]) {
     319             :       /* find left interval */
     320         689 :       for (k=i;k>=0;k--) {
     321         689 :         if (PetscAbsReal(mu[k])<eta || which[k]) break;
     322           0 :         else which[k] = PETSC_TRUE;
     323             :       }
     324             :       /* find right interval */
     325        1189 :       for (k=i+1;k<j;k++) {
     326        1159 :         if (PetscAbsReal(mu[k])<eta || which[k]) break;
     327         500 :         else which[k] = PETSC_TRUE;
     328             :       }
     329             :     }
     330             :   }
     331          53 :   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         168 :   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        1428 :   for (j=k;j<m;j++) {
     371        1352 :     PetscCall(BVMatMultColumn(eps->V,Op,j));
     372        1352 :     if (fro) {
     373             :       /* Lanczos step with full reorthogonalization */
     374         434 :       PetscCall(BVOrthogonalizeColumn(eps->V,j+1,hwork,&norm,breakdown));
     375         434 :       alpha[j] = PetscRealPart(hwork[j]);
     376             :     } else {
     377             :       /* Lanczos step */
     378         918 :       which[j] = PETSC_TRUE;
     379         918 :       if (j-2>=k) which[j-2] = PETSC_FALSE;
     380         918 :       PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which,hwork,&norm,breakdown));
     381         918 :       alpha[j] = PetscRealPart(hwork[j]);
     382         918 :       beta[j] = norm;
     383             : 
     384             :       /* Estimate ||A|| if needed */
     385         918 :       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         918 :       reorth = PETSC_FALSE;
     392         918 :       if (j>k) {
     393         842 :         update_omega(omega,omega_old,j,alpha,beta-1,eps1,anorm);
     394        9228 :         for (i=0;i<j-k;i++) {
     395        7544 :           if (PetscAbsReal(omega[i]) > delta) reorth = PETSC_TRUE;
     396             :         }
     397             :       }
     398         918 :       if (reorth || force_reorth) {
     399         298 :         for (i=0;i<k;i++) which2[i] = PETSC_FALSE;
     400        2905 :         for (i=k;i<=j;i++) which2[i] = PETSC_TRUE;
     401         173 :         if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
     402             :           /* Periodic reorthogonalization */
     403          85 :           if (force_reorth) force_reorth = PETSC_FALSE;
     404          50 :           else force_reorth = PETSC_TRUE;
     405        1338 :           for (i=0;i<j-k;i++) omega[i] = eps1;
     406             :         } else {
     407             :           /* Partial reorthogonalization */
     408          88 :           if (force_reorth) force_reorth = PETSC_FALSE;
     409             :           else {
     410          53 :             force_reorth = PETSC_TRUE;
     411          53 :             compute_int(which2+k,omega,j-k,delta,eta);
     412         896 :             for (i=0;i<j-k;i++) {
     413         790 :               if (which2[i+k]) omega[i] = eps1;
     414             :             }
     415             :           }
     416             :         }
     417         173 :         PetscCall(BVOrthogonalizeSomeColumn(eps->V,j+1,which2,hwork,&norm,breakdown));
     418             :       }
     419             :     }
     420             : 
     421        1352 :     if (PetscUnlikely(*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON)) {
     422           0 :       *M = j+1;
     423           0 :       break;
     424             :     }
     425        1352 :     if (!fro && norm*delta < anorm*eps1) {
     426          26 :       fro = PETSC_TRUE;
     427          26 :       PetscCall(PetscInfo(eps,"Switching to full reorthogonalization at iteration %" PetscInt_FMT "\n",eps->its));
     428             :     }
     429        1352 :     beta[j] = norm;
     430        1352 :     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         938 : static PetscErrorCode EPSBasicLanczos(EPS eps,PetscInt k,PetscInt *m,PetscReal *betam,PetscBool *breakdown,PetscReal anorm)
     456             : {
     457         938 :   EPS_LANCZOS        *lanczos = (EPS_LANCZOS*)eps->data;
     458         938 :   PetscScalar        *T;
     459         938 :   PetscInt           i,n=*m,ld;
     460         938 :   PetscReal          *alpha,*beta;
     461         938 :   BVOrthogRefineType orthog_ref;
     462         938 :   Mat                Op,M;
     463             : 
     464         938 :   PetscFunctionBegin;
     465         938 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     466         938 :   switch (lanczos->reorthog) {
     467         514 :     case EPS_LANCZOS_REORTHOG_LOCAL:
     468         514 :       PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     469         514 :       beta = alpha + ld;
     470         514 :       PetscCall(EPSLocalLanczos(eps,alpha,beta,k,m,breakdown));
     471         514 :       *betam = beta[*m-1];
     472         514 :       PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     473             :       break;
     474         287 :     case EPS_LANCZOS_REORTHOG_FULL:
     475         287 :       PetscCall(STGetOperator(eps->st,&Op));
     476         287 :       PetscCall(DSGetMat(eps->ds,DS_MAT_T,&M));
     477         287 :       PetscCall(BVMatLanczos(eps->V,Op,M,k,m,betam,breakdown));
     478         287 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&M));
     479         287 :       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          23 :     case EPS_LANCZOS_REORTHOG_DELAYED:
     497          23 :       PetscCall(PetscMalloc1(n*n,&T));
     498          23 :       PetscCall(BVGetOrthogonalization(eps->V,NULL,&orthog_ref,NULL,NULL));
     499          23 :       if (orthog_ref == BV_ORTHOG_REFINE_NEVER) PetscCall(EPSDelayedArnoldi1(eps,T,n,k,m,betam,breakdown));
     500          23 :       else PetscCall(EPSDelayedArnoldi(eps,T,n,k,m,betam,breakdown));
     501          23 :       n = *m;
     502          23 :       PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     503          23 :       beta = alpha + ld;
     504         403 :       for (i=k;i<n-1;i++) {
     505         380 :         alpha[i] = PetscRealPart(T[n*i+i]);
     506         380 :         beta[i] = PetscRealPart(T[n*i+i+1]);
     507             :       }
     508          23 :       alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
     509          23 :       beta[n-1] = *betam;
     510          23 :       PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     511          23 :       PetscCall(PetscFree(T));
     512             :       break;
     513             :   }
     514         938 :   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         974 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     540         938 :     eps->its++;
     541             : 
     542             :     /* Compute an ncv-step Lanczos factorization */
     543         938 :     n = PetscMin(nconv+eps->mpd,ncv);
     544         938 :     PetscCall(DSSetDimensions(eps->ds,n,nconv,PETSC_DEFAULT));
     545         938 :     PetscCall(EPSBasicLanczos(eps,nconv,&n,&beta,&breakdown,anorm));
     546         938 :     PetscCall(DSSetDimensions(eps->ds,n,nconv,0));
     547         938 :     PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
     548         938 :     PetscCall(BVSetActiveColumns(eps->V,nconv,n));
     549             : 
     550             :     /* Solve projected problem */
     551         938 :     PetscCall(DSSolve(eps->ds,ritz,NULL));
     552         938 :     PetscCall(DSSort(eps->ds,ritz,NULL,NULL,NULL,NULL));
     553         938 :     PetscCall(DSSynchronize(eps->ds,ritz,NULL));
     554             : 
     555             :     /* Estimate ||A|| */
     556       13337 :     for (i=nconv;i<n;i++)
     557       12830 :       anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
     558             : 
     559             :     /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
     560         938 :     PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
     561       13337 :     for (i=nconv;i<n;i++) {
     562       12399 :       resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
     563       12399 :       PetscCall((*eps->converged)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->convergedctx));
     564       12399 :       if (bnd[i]<eps->tol) conv[i] = 'C';
     565       12230 :       else conv[i] = 'N';
     566             :     }
     567         938 :     PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
     568             : 
     569             :     /* purge repeated ritz values */
     570         938 :     if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
     571        6966 :       for (i=nconv+1;i<n;i++) {
     572        6452 :         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         938 :     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        1084 :       while (restart<n && conv[restart] != 'N') restart++;
     581         937 :       if (restart >= n) {
     582           0 :         breakdown = PETSC_TRUE;
     583             :       } else {
     584       12234 :         for (i=restart+1;i<n;i++) {
     585       11297 :           if (conv[i] == 'N') {
     586       11293 :             PetscCall(SlepcSCCompare(eps->sc,ritz[restart],0.0,ritz[i],0.0,&r));
     587       11293 :             if (r>0) restart = i;
     588             :           }
     589             :         }
     590         937 :         PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
     591         937 :         PetscCall(BVMultColumn(eps->V,1.0,0.0,n,Y+restart*ld+nconv));
     592         937 :         PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
     593             :       }
     594             :     }
     595             : 
     596             :     /* Count and put converged eigenvalues first */
     597       13337 :     for (i=nconv;i<n;i++) perm[i] = i;
     598        1104 :     for (k=nconv;k<n;k++) {
     599        1103 :       if (conv[perm[k]] != 'C') {
     600         940 :         j = k + 1;
     601       12240 :         while (j<n && conv[perm[j]] != 'C') j++;
     602         940 :         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         938 :     PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Y));
     609        1104 :     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         938 :     PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Y));
     626             : 
     627             :     /* compute converged eigenvectors */
     628         938 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
     629         938 :     PetscCall(BVMultInPlace(eps->V,U,nconv,k));
     630         938 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
     631             : 
     632             :     /* purge spurious ritz values */
     633         938 :     if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
     634         589 :       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         588 :       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       13337 :     for (i=nconv;i<n;i++) {
     668       12399 :       eps->eigr[i] = ritz[i];
     669       12399 :       eps->errest[i] = bnd[i];
     670             :     }
     671         938 :     nconv = k;
     672         938 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n));
     673         938 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,nconv,eps->nev,&eps->reason,eps->stoppingctx));
     674             : 
     675         938 :     if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
     676         902 :       PetscCall(BVCopyColumn(eps->V,n,nconv));
     677         902 :       if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
     678             :         /* Reorthonormalize restart vector */
     679         496 :         PetscCall(BVOrthonormalizeColumn(eps->V,nconv,PETSC_FALSE,NULL,&breakdown));
     680             :       }
     681         902 :       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         902 :       if (PetscUnlikely(breakdown)) { /* give up */
     687           0 :         eps->reason = EPS_DIVERGED_BREAKDOWN;
     688         974 :         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