LCOV - code coverage report
Current view: top level - nep/impls/nleigs - nleigs.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 1180 1213 97.3 %
Date: 2024-03-28 00:28:38 Functions: 58 58 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    SLEPc nonlinear eigensolver: "nleigs"
      12             : 
      13             :    Method: NLEIGS
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Fully rational Krylov method for nonlinear eigenvalue problems.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] S. Guttel et al., "NLEIGS: A class of robust fully rational Krylov
      22             :            method for nonlinear eigenvalue problems", SIAM J. Sci. Comput.
      23             :            36(6):A2842-A2864, 2014.
      24             : */
      25             : 
      26             : #include <slepc/private/nepimpl.h>         /*I "slepcnep.h" I*/
      27             : #include <slepcblaslapack.h>
      28             : #include "nleigs.h"
      29             : 
      30       19380 : PetscErrorCode NEPNLEIGSBackTransform(PetscObject ob,PetscInt n,PetscScalar *valr,PetscScalar *vali)
      31             : {
      32       19380 :   NEP         nep;
      33       19380 :   PetscInt    j;
      34             : #if !defined(PETSC_USE_COMPLEX)
      35             :   PetscScalar t;
      36             : #endif
      37             : 
      38       19380 :   PetscFunctionBegin;
      39       19380 :   nep = (NEP)ob;
      40             : #if !defined(PETSC_USE_COMPLEX)
      41             :   for (j=0;j<n;j++) {
      42             :     if (vali[j] == 0) valr[j] = 1.0 / valr[j] + nep->target;
      43             :     else {
      44             :       t = valr[j] * valr[j] + vali[j] * vali[j];
      45             :       valr[j] = valr[j] / t + nep->target;
      46             :       vali[j] = - vali[j] / t;
      47             :     }
      48             :   }
      49             : #else
      50       59427 :   for (j=0;j<n;j++) {
      51       40047 :     valr[j] = 1.0 / valr[j] + nep->target;
      52             :   }
      53             : #endif
      54       19380 :   PetscFunctionReturn(PETSC_SUCCESS);
      55             : }
      56             : 
      57             : /* Computes the roots of a polynomial */
      58           3 : static PetscErrorCode NEPNLEIGSAuxiliarPRootFinder(PetscInt deg,PetscScalar *polcoeffs,PetscScalar *wr,PetscScalar *wi,PetscBool *avail)
      59             : {
      60           3 :   PetscScalar    *C;
      61           3 :   PetscBLASInt   n_,lwork;
      62           3 :   PetscInt       i;
      63             : #if defined(PETSC_USE_COMPLEX)
      64           3 :   PetscReal      *rwork=NULL;
      65             : #endif
      66           3 :   PetscScalar    *work;
      67           3 :   PetscBLASInt   info;
      68             : 
      69           3 :   PetscFunctionBegin;
      70           3 :   *avail = PETSC_TRUE;
      71           3 :   if (deg>0) {
      72           3 :     PetscCall(PetscCalloc1(deg*deg,&C));
      73           3 :     PetscCall(PetscBLASIntCast(deg,&n_));
      74           3 :     for (i=0;i<deg-1;i++) {
      75           0 :       C[(deg+1)*i+1]   = 1.0;
      76           0 :       C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
      77             :     }
      78           3 :     C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
      79           3 :     PetscCall(PetscBLASIntCast(3*deg,&lwork));
      80             : 
      81           3 :     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
      82             : #if !defined(PETSC_USE_COMPLEX)
      83             :     PetscCall(PetscMalloc1(lwork,&work));
      84             :     PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
      85             :     if (info) *avail = PETSC_FALSE;
      86             :     PetscCall(PetscFree(work));
      87             : #else
      88           3 :     PetscCall(PetscMalloc2(2*deg,&rwork,lwork,&work));
      89           3 :     PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
      90           3 :     if (info) *avail = PETSC_FALSE;
      91           3 :     PetscCall(PetscFree2(rwork,work));
      92             : #endif
      93           3 :     PetscCall(PetscFPTrapPop());
      94           3 :     PetscCall(PetscFree(C));
      95             :   }
      96           3 :   PetscFunctionReturn(PETSC_SUCCESS);
      97             : }
      98             : 
      99          15 : static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
     100             : {
     101          15 :   PetscInt i,j;
     102             : 
     103          15 :   PetscFunctionBegin;
     104          45 :   for (i=0;i<nin;i++) {
     105          30 :     if (max && *nout>=max) break;
     106          30 :     pout[(*nout)++] = pin[i];
     107         125 :     for (j=0;j<*nout-1;j++)
     108          95 :       if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
     109           0 :         (*nout)--;
     110           0 :         break;
     111             :       }
     112             :   }
     113          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     114             : }
     115             : 
     116          11 : static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
     117             : {
     118          11 :   FNCombineType  ctype;
     119          11 :   FN             f1,f2;
     120          11 :   PetscInt       i,nq,nisol1,nisol2;
     121          11 :   PetscScalar    *qcoeff,*wr,*wi,*isol1,*isol2;
     122          11 :   PetscBool      flg,avail,rat1,rat2;
     123             : 
     124          11 :   PetscFunctionBegin;
     125          11 :   *rational = PETSC_FALSE;
     126          11 :   PetscCall(PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg));
     127          11 :   if (flg) {
     128          10 :     *rational = PETSC_TRUE;
     129          10 :     PetscCall(FNRationalGetDenominator(f,&nq,&qcoeff));
     130          10 :     if (nq>1) {
     131           3 :       PetscCall(PetscMalloc2(nq-1,&wr,nq-1,&wi));
     132           3 :       PetscCall(NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail));
     133           3 :       if (avail) {
     134           3 :         PetscCall(PetscCalloc1(nq-1,isol));
     135           3 :         *nisol = 0;
     136           6 :         for (i=0;i<nq-1;i++)
     137             : #if !defined(PETSC_USE_COMPLEX)
     138             :           if (wi[i]==0)
     139             : #endif
     140           3 :             (*isol)[(*nisol)++] = wr[i];
     141           3 :         nq = *nisol; *nisol = 0;
     142           6 :         for (i=0;i<nq;i++) wr[i] = (*isol)[i];
     143           3 :         PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol,0));
     144           3 :         PetscCall(PetscFree2(wr,wi));
     145           0 :       } else { *nisol=0; *isol = NULL; }
     146           7 :     } else { *nisol = 0; *isol = NULL; }
     147          10 :     PetscCall(PetscFree(qcoeff));
     148             :   }
     149          11 :   PetscCall(PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg));
     150          11 :   if (flg) {
     151           1 :     PetscCall(FNCombineGetChildren(f,&ctype,&f1,&f2));
     152           1 :     if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
     153           1 :       PetscCall(NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1));
     154           1 :       PetscCall(NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2));
     155           1 :       if (nisol1+nisol2>0) {
     156           1 :         PetscCall(PetscCalloc1(nisol1+nisol2,isol));
     157           1 :         *nisol = 0;
     158           1 :         PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol,0));
     159           1 :         PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol,0));
     160             :       }
     161           1 :       *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
     162           1 :       PetscCall(PetscFree(isol1));
     163           1 :       PetscCall(PetscFree(isol2));
     164             :     }
     165             :   }
     166          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     167             : }
     168             : 
     169           3 : static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
     170             : {
     171           3 :   PetscInt       nt,i,nisol;
     172           3 :   FN             f;
     173           3 :   PetscScalar    *isol;
     174           3 :   PetscBool      rat;
     175             : 
     176           3 :   PetscFunctionBegin;
     177           3 :   *rational = PETSC_TRUE;
     178           3 :   *ndptx = 0;
     179           3 :   PetscCall(NEPGetSplitOperatorInfo(nep,&nt,NULL));
     180          12 :   for (i=0;i<nt;i++) {
     181           9 :     PetscCall(NEPGetSplitOperatorTerm(nep,i,NULL,&f));
     182           9 :     PetscCall(NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat));
     183           9 :     if (nisol) {
     184           3 :       PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0));
     185           3 :       PetscCall(PetscFree(isol));
     186             :     }
     187           9 :     *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
     188             :   }
     189           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     190             : }
     191             : 
     192             : #if defined(SLEPC_MISSING_LAPACK_GGEV3)
     193             : #define LAPGEEV "ggev"
     194             : #else
     195             : #define LAPGEEV "ggev3"
     196             : #endif
     197             : 
     198             : /*  Adaptive Anderson-Antoulas algorithm */
     199          25 : static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
     200             : {
     201          25 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     202          25 :   PetscScalar    mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
     203          25 :   PetscScalar    *N,*D;
     204          25 :   PetscReal      *S,norm,err,*R;
     205          25 :   PetscInt       i,k,j,idx=0,cont;
     206          25 :   PetscBLASInt   n_,m_,lda_,lwork,info,one=1;
     207             : #if defined(PETSC_USE_COMPLEX)
     208          25 :   PetscReal      *rwork;
     209             : #endif
     210             : 
     211          25 :   PetscFunctionBegin;
     212          25 :   PetscCall(PetscBLASIntCast(8*ndpt,&lwork));
     213          25 :   PetscCall(PetscMalloc5(ndpt,&R,ndpt,&z,ndpt,&f,ndpt*ndpt,&C,ndpt,&ww));
     214          25 :   PetscCall(PetscMalloc6(ndpt*ndpt,&A,ndpt,&S,ndpt*ndpt,&VT,lwork,&work,ndpt,&D,ndpt,&N));
     215             : #if defined(PETSC_USE_COMPLEX)
     216          25 :   PetscCall(PetscMalloc1(8*ndpt,&rwork));
     217             : #endif
     218          25 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     219             :   norm = 0.0;
     220        2525 :   for (i=0;i<ndpt;i++) {
     221        2500 :     mean += F[i];
     222        3613 :     norm = PetscMax(PetscAbsScalar(F[i]),norm);
     223             :   }
     224          25 :   mean /= ndpt;
     225          25 :   PetscCall(PetscBLASIntCast(ndpt,&lda_));
     226        2525 :   for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
     227             :   /* next support point */
     228             :   err = 0.0;
     229        2525 :   for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
     230          64 :   for (k=0;k<ndpt-1;k++) {
     231          64 :     z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
     232             :     /* next column of Cauchy matrix */
     233        6464 :     for (i=0;i<ndpt;i++) {
     234        6400 :       C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
     235             :     }
     236             : 
     237          64 :     PetscCall(PetscArrayzero(A,ndpt*ndpt));
     238             :     cont = 0;
     239        6464 :     for (i=0;i<ndpt;i++) {
     240        6400 :       if (R[i]!=-1.0) {
     241       21058 :         for (j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
     242        6246 :         cont++;
     243             :       }
     244             :     }
     245          64 :     PetscCall(PetscBLASIntCast(cont,&m_));
     246          64 :     PetscCall(PetscBLASIntCast(k+1,&n_));
     247             : #if defined(PETSC_USE_COMPLEX)
     248          64 :     PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
     249             : #else
     250             :     PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
     251             : #endif
     252          64 :     SlepcCheckLapackInfo("gesvd",info);
     253         218 :     for (i=0;i<=k;i++) {
     254         154 :       ww[i] = PetscConj(VT[i*ndpt+k]);
     255         154 :       D[i] = ww[i]*f[i];
     256             :     }
     257          64 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
     258          64 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
     259        6464 :     for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
     260             :     /* next support point */
     261             :     err = 0.0;
     262        6464 :     for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
     263          64 :     if (err <= ctx->ddtol*norm) break;
     264             :   }
     265             : 
     266          25 :   PetscCheck(k<ndpt-1,PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in general problem");
     267             :   /* poles */
     268          25 :   PetscCall(PetscArrayzero(C,ndpt*ndpt));
     269          25 :   PetscCall(PetscArrayzero(A,ndpt*ndpt));
     270          89 :   for (i=0;i<=k;i++) {
     271          64 :     C[i+ndpt*i] = 1.0;
     272          64 :     A[(i+1)*ndpt] = ww[i];
     273          64 :     A[i+1] = 1.0;
     274          64 :     A[i+1+(i+1)*ndpt] = z[i];
     275             :   }
     276          25 :   C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
     277          25 :   n_++;
     278             : #if defined(PETSC_USE_COMPLEX)
     279          25 :   PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
     280             : #else
     281             :   PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
     282             : #endif
     283          25 :   SlepcCheckLapackInfo(LAPGEEV,info);
     284             :   cont = 0.0;
     285         114 :   for (i=0;i<n_;i++) if (N[i]!=0.0) {
     286          33 :     dxi[cont++] = D[i]/N[i];
     287             :   }
     288          25 :   *ndptx = cont;
     289          25 :   PetscCall(PetscFPTrapPop());
     290          25 :   PetscCall(PetscFree5(R,z,f,C,ww));
     291          25 :   PetscCall(PetscFree6(A,S,VT,work,D,N));
     292             : #if defined(PETSC_USE_COMPLEX)
     293          25 :   PetscCall(PetscFree(rwork));
     294             : #endif
     295          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     296             : }
     297             : 
     298             : /*  Singularities using Adaptive Anderson-Antoulas algorithm */
     299          12 : static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
     300             : {
     301          12 :   Vec            u,v,w;
     302          12 :   PetscRandom    rand=NULL;
     303          12 :   PetscScalar    *F,*isol;
     304          12 :   PetscInt       i,k,nisol,nt;
     305          12 :   Mat            T;
     306          12 :   FN             f;
     307             : 
     308          12 :   PetscFunctionBegin;
     309          12 :   PetscCall(PetscMalloc1(ndpt,&F));
     310          12 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     311           6 :     PetscCall(PetscMalloc1(ndpt,&isol));
     312           6 :     *ndptx = 0;
     313           6 :     PetscCall(NEPGetSplitOperatorInfo(nep,&nt,NULL));
     314           6 :     nisol = *ndptx;
     315          25 :     for (k=0;k<nt;k++) {
     316          19 :       PetscCall(NEPGetSplitOperatorTerm(nep,k,NULL,&f));
     317        1919 :       for (i=0;i<ndpt;i++) PetscCall(FNEvaluateFunction(f,ds[i],&F[i]));
     318          19 :       PetscCall(NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol));
     319          19 :       if (nisol) PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt));
     320             :     }
     321           6 :     PetscCall(PetscFree(isol));
     322             :   } else {
     323           6 :     PetscCall(MatCreateVecs(nep->function,&u,NULL));
     324           6 :     PetscCall(VecDuplicate(u,&v));
     325           6 :     PetscCall(VecDuplicate(u,&w));
     326           6 :     if (nep->V) PetscCall(BVGetRandomContext(nep->V,&rand));
     327           6 :     PetscCall(VecSetRandom(u,rand));
     328           6 :     PetscCall(VecNormalize(u,NULL));
     329           6 :     PetscCall(VecSetRandom(v,rand));
     330           6 :     PetscCall(VecNormalize(v,NULL));
     331           6 :     T = nep->function;
     332         606 :     for (i=0;i<ndpt;i++) {
     333         600 :       PetscCall(NEPComputeFunction(nep,ds[i],T,T));
     334         600 :       PetscCall(MatMult(T,v,w));
     335         600 :       PetscCall(VecDot(w,u,&F[i]));
     336             :     }
     337           6 :     PetscCall(NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi));
     338           6 :     PetscCall(VecDestroy(&u));
     339           6 :     PetscCall(VecDestroy(&v));
     340           6 :     PetscCall(VecDestroy(&w));
     341             :   }
     342          12 :   PetscCall(PetscFree(F));
     343          12 :   PetscFunctionReturn(PETSC_SUCCESS);
     344             : }
     345             : 
     346          27 : static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
     347             : {
     348          27 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     349          27 :   PetscInt       i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
     350          27 :   PetscScalar    *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
     351          27 :   PetscReal      maxnrs,minnrxi;
     352          27 :   PetscBool      rational;
     353             : #if !defined(PETSC_USE_COMPLEX)
     354             :   PetscReal      a,b,h;
     355             : #endif
     356             : 
     357          27 :   PetscFunctionBegin;
     358          27 :   if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
     359          27 :   PetscCall(PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi));
     360             : 
     361             :   /* Discretize the target region boundary */
     362          27 :   PetscCall(RGComputeContour(nep->rg,ndpt,ds,dsi));
     363             : #if !defined(PETSC_USE_COMPLEX)
     364             :   for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
     365             :   if (i<ndpt) {
     366             :     PetscCheck(nep->problem_type==NEP_RATIONAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NLEIGS with real arithmetic requires the target set to be included in the real axis");
     367             :     /* Select a segment in the real axis */
     368             :     PetscCall(RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL));
     369             :     PetscCheck(a>-PETSC_MAX_REAL && b<PETSC_MAX_REAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"NLEIGS requires a bounded target set");
     370             :     h = (b-a)/ndpt;
     371             :     for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
     372             :   }
     373             : #endif
     374             :   /* Discretize the singularity region */
     375          27 :   if (ctx->computesingularities) PetscCall(ctx->computesingularities(nep,&ndptx,dxi,ctx->singularitiesctx));
     376             :   else {
     377          15 :     if (nep->problem_type==NEP_RATIONAL) {
     378           3 :       PetscCall(NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational));
     379           3 :       PetscCheck(rational,PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
     380             :     } else {
     381             :       /* AAA algorithm */
     382          12 :       PetscCall(NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi));
     383             :     }
     384             :   }
     385             :   /* Look for Leja-Bagby points in the discretization sets */
     386          27 :   s[0]  = ds[0];
     387          27 :   xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
     388          27 :   PetscCheck(PetscAbsScalar(xi[0])>=10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"Singularity point 0 is nearly zero: %g; consider removing the singularity or shifting the problem",(double)PetscAbsScalar(xi[0]));
     389          27 :   beta[0] = 1.0; /* scaling factors are also computed here */
     390      151227 :   for (i=0;i<ndpt;i++) {
     391      151200 :     nrs[i] = 1.0;
     392      151200 :     nrxi[i] = 1.0;
     393             :   }
     394        2325 :   for (k=1;k<ctx->ddmaxit;k++) {
     395             :     maxnrs = 0.0;
     396    12062598 :     minnrxi = PETSC_MAX_REAL;
     397    12062598 :     for (i=0;i<ndpt;i++) {
     398    12060300 :       nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
     399    12060300 :       if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
     400             :     }
     401        2298 :     if (ndptx>k) {
     402     8980198 :       for (i=1;i<ndptx;i++) {
     403     8979279 :         nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
     404     8979279 :         if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
     405             :       }
     406         919 :       PetscCheck(PetscAbsScalar(xi[k])>=10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"Singularity point %" PetscInt_FMT " is nearly zero: %g; consider removing the singularity or shifting the problem",k,(double)PetscAbsScalar(xi[k]));
     407        1379 :     } else xi[k] = PETSC_INFINITY;
     408        2298 :     beta[k] = maxnrs;
     409             :   }
     410          27 :   PetscCall(PetscFree5(ds,dsi,dxi,nrs,nrxi));
     411          27 :   PetscFunctionReturn(PETSC_SUCCESS);
     412             : }
     413             : 
     414        1240 : PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
     415             : {
     416        1240 :   NEP_NLEIGS  *ctx=(NEP_NLEIGS*)nep->data;
     417        1240 :   PetscInt    i;
     418        1240 :   PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
     419             : 
     420        1240 :   PetscFunctionBegin;
     421        1240 :   b[0] = 1.0/beta[0];
     422       12584 :   for (i=0;i<k;i++) {
     423       11344 :     b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
     424             :   }
     425        1240 :   PetscFunctionReturn(PETSC_SUCCESS);
     426             : }
     427             : 
     428       18105 : static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
     429             : {
     430       18105 :   NEP_NLEIGS_MATSHELL *ctx;
     431       18105 :   PetscInt            i;
     432             : 
     433       18105 :   PetscFunctionBeginUser;
     434       18105 :   PetscCall(MatShellGetContext(A,&ctx));
     435       18105 :   PetscCall(MatMult(ctx->A[0],x,y));
     436       18105 :   if (ctx->coeff[0]!=1.0) PetscCall(VecScale(y,ctx->coeff[0]));
     437      211277 :   for (i=1;i<ctx->nmat;i++) {
     438      193172 :     PetscCall(MatMult(ctx->A[i],x,ctx->t));
     439      193172 :     PetscCall(VecAXPY(y,ctx->coeff[i],ctx->t));
     440             :   }
     441       18105 :   PetscFunctionReturn(PETSC_SUCCESS);
     442             : }
     443             : 
     444       16401 : static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
     445             : {
     446       16401 :   NEP_NLEIGS_MATSHELL *ctx;
     447       16401 :   PetscInt            i;
     448             : 
     449       16401 :   PetscFunctionBeginUser;
     450       16401 :   PetscCall(MatShellGetContext(A,&ctx));
     451       16401 :   PetscCall(MatMultTranspose(ctx->A[0],x,y));
     452       16401 :   if (ctx->coeff[0]!=1.0) PetscCall(VecScale(y,ctx->coeff[0]));
     453      192422 :   for (i=1;i<ctx->nmat;i++) {
     454      176021 :     PetscCall(MatMultTranspose(ctx->A[i],x,ctx->t));
     455      176021 :     PetscCall(VecAXPY(y,ctx->coeff[i],ctx->t));
     456             :   }
     457       16401 :   PetscFunctionReturn(PETSC_SUCCESS);
     458             : }
     459             : 
     460           6 : static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
     461             : {
     462           6 :   NEP_NLEIGS_MATSHELL *ctx;
     463           6 :   PetscInt            i;
     464             : 
     465           6 :   PetscFunctionBeginUser;
     466           6 :   PetscCall(MatShellGetContext(A,&ctx));
     467           6 :   PetscCall(MatGetDiagonal(ctx->A[0],diag));
     468           6 :   if (ctx->coeff[0]!=1.0) PetscCall(VecScale(diag,ctx->coeff[0]));
     469          72 :   for (i=1;i<ctx->nmat;i++) {
     470          66 :     PetscCall(MatGetDiagonal(ctx->A[i],ctx->t));
     471          66 :     PetscCall(VecAXPY(diag,ctx->coeff[i],ctx->t));
     472             :   }
     473           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     474             : }
     475             : 
     476           3 : static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
     477             : {
     478           3 :   PetscInt            m,n,M,N,i;
     479           3 :   NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
     480           3 :   void                (*fun)(void);
     481             : 
     482           3 :   PetscFunctionBeginUser;
     483           3 :   PetscCall(MatShellGetContext(A,&ctx));
     484           3 :   PetscCall(PetscNew(&ctxnew));
     485           3 :   ctxnew->nmat = ctx->nmat;
     486           3 :   ctxnew->maxnmat = ctx->maxnmat;
     487           3 :   PetscCall(PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff));
     488           6 :   for (i=0;i<ctx->nmat;i++) {
     489           3 :     PetscCall(PetscObjectReference((PetscObject)ctx->A[i]));
     490           3 :     ctxnew->A[i] = ctx->A[i];
     491           3 :     ctxnew->coeff[i] = ctx->coeff[i];
     492             :   }
     493           3 :   PetscCall(MatGetSize(ctx->A[0],&M,&N));
     494           3 :   PetscCall(MatGetLocalSize(ctx->A[0],&m,&n));
     495           3 :   PetscCall(VecDuplicate(ctx->t,&ctxnew->t));
     496           3 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctxnew,B));
     497           3 :   PetscCall(MatShellSetManageScalingShifts(*B));
     498           3 :   PetscCall(MatShellGetOperation(A,MATOP_MULT,&fun));
     499           3 :   PetscCall(MatShellSetOperation(*B,MATOP_MULT,fun));
     500           3 :   PetscCall(MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun));
     501           3 :   PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun));
     502           3 :   PetscCall(MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun));
     503           3 :   PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun));
     504           3 :   PetscCall(MatShellGetOperation(A,MATOP_DUPLICATE,&fun));
     505           3 :   PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,fun));
     506           3 :   PetscCall(MatShellGetOperation(A,MATOP_DESTROY,&fun));
     507           3 :   PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,fun));
     508           3 :   PetscCall(MatShellGetOperation(A,MATOP_AXPY,&fun));
     509           3 :   PetscCall(MatShellSetOperation(*B,MATOP_AXPY,fun));
     510           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     511             : }
     512             : 
     513          75 : static PetscErrorCode MatDestroy_Fun(Mat A)
     514             : {
     515          75 :   NEP_NLEIGS_MATSHELL *ctx;
     516          75 :   PetscInt            i;
     517             : 
     518          75 :   PetscFunctionBeginUser;
     519          75 :   if (A) {
     520          75 :     PetscCall(MatShellGetContext(A,&ctx));
     521         909 :     for (i=0;i<ctx->nmat;i++) PetscCall(MatDestroy(&ctx->A[i]));
     522          75 :     PetscCall(VecDestroy(&ctx->t));
     523          75 :     PetscCall(PetscFree2(ctx->A,ctx->coeff));
     524          75 :     PetscCall(PetscFree(ctx));
     525             :   }
     526          75 :   PetscFunctionReturn(PETSC_SUCCESS);
     527             : }
     528             : 
     529         759 : static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
     530             : {
     531         759 :   NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
     532         759 :   PetscInt            i,j;
     533         759 :   PetscBool           found;
     534             : 
     535         759 :   PetscFunctionBeginUser;
     536         759 :   PetscCall(MatShellGetContext(Y,&ctxY));
     537         759 :   PetscCall(MatShellGetContext(X,&ctxX));
     538        6831 :   for (i=0;i<ctxX->nmat;i++) {
     539             :     found = PETSC_FALSE;
     540       48579 :     for (j=0;!found&&j<ctxY->nmat;j++) {
     541       42507 :       if (ctxX->A[i]==ctxY->A[j]) {
     542        5313 :         found = PETSC_TRUE;
     543        5313 :         ctxY->coeff[j] += a*ctxX->coeff[i];
     544             :       }
     545             :     }
     546        6072 :     if (!found) {
     547         759 :       ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
     548         759 :       ctxY->A[ctxY->nmat++] = ctxX->A[i];
     549        6072 :       PetscCall(PetscObjectReference((PetscObject)ctxX->A[i]));
     550             :     }
     551             :   }
     552         759 :   PetscFunctionReturn(PETSC_SUCCESS);
     553             : }
     554             : 
     555          63 : static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
     556             : {
     557          63 :   NEP_NLEIGS_MATSHELL *ctx;
     558          63 :   PetscInt            i;
     559             : 
     560          63 :   PetscFunctionBeginUser;
     561          63 :   PetscCall(MatShellGetContext(M,&ctx));
     562         819 :   for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
     563          63 :   PetscFunctionReturn(PETSC_SUCCESS);
     564             : }
     565             : 
     566          72 : static PetscErrorCode NLEIGSMatToMatShellArray(Mat A,Mat *Ms,PetscInt maxnmat)
     567             : {
     568          72 :   NEP_NLEIGS_MATSHELL *ctx;
     569          72 :   PetscInt            m,n,M,N;
     570          72 :   PetscBool           has;
     571             : 
     572          72 :   PetscFunctionBegin;
     573          72 :   PetscCall(MatHasOperation(A,MATOP_DUPLICATE,&has));
     574          72 :   PetscCheck(has,PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"MatDuplicate operation required");
     575          72 :   PetscCall(PetscNew(&ctx));
     576          72 :   ctx->maxnmat = maxnmat;
     577          72 :   PetscCall(PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff));
     578          72 :   PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&ctx->A[0]));
     579          72 :   ctx->nmat = 1;
     580          72 :   ctx->coeff[0] = 1.0;
     581          72 :   PetscCall(MatCreateVecs(A,&ctx->t,NULL));
     582          72 :   PetscCall(MatGetSize(A,&M,&N));
     583          72 :   PetscCall(MatGetLocalSize(A,&m,&n));
     584          72 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctx,Ms));
     585          72 :   PetscCall(MatShellSetManageScalingShifts(*Ms));
     586          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun));
     587          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun));
     588          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun));
     589          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun));
     590          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun));
     591          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun));
     592          72 :   PetscCall(MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun));
     593          72 :   PetscFunctionReturn(PETSC_SUCCESS);
     594             : }
     595             : 
     596             : /*
     597             :    MatIsShellAny - returns true if any of the n matrices is a shell matrix
     598             :  */
     599          16 : static PetscErrorCode MatIsShellAny(Mat *A,PetscInt n,PetscBool *shell)
     600             : {
     601          16 :   PetscInt       i;
     602          16 :   PetscBool      flg;
     603             : 
     604          16 :   PetscFunctionBegin;
     605          16 :   *shell = PETSC_FALSE;
     606          52 :   for (i=0;i<n;i++) {
     607          39 :     PetscCall(MatIsShell(A[i],&flg));
     608          39 :     if (flg) { *shell = PETSC_TRUE; break; }
     609             :   }
     610          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     611             : }
     612             : 
     613          16 : static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
     614             : {
     615          16 :   PetscErrorCode ierr;
     616          16 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     617          16 :   PetscInt       k,j,i,maxnmat,nmax;
     618          16 :   PetscReal      norm0,norm,*matnorm;
     619          16 :   PetscScalar    *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
     620          16 :   Mat            T,P,Ts,K,H;
     621          16 :   PetscBool      shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
     622          16 :   PetscBLASInt   n_;
     623             : 
     624          16 :   PetscFunctionBegin;
     625          16 :   nmax = ctx->ddmaxit;
     626          16 :   PetscCall(PetscMalloc1(nep->nt*nmax,&ctx->coeffD));
     627          16 :   PetscCall(PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm));
     628          52 :   for (j=0;j<nep->nt;j++) {
     629          39 :     PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm));
     630          39 :     if (!hasmnorm) break;
     631          36 :     PetscCall(MatNorm(nep->A[j],NORM_INFINITY,matnorm+j));
     632             :   }
     633             :   /* Try matrix functions scheme */
     634          16 :   PetscCall(PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH));
     635        1315 :   for (i=0;i<nmax-1;i++) {
     636        1299 :     pK[(nmax+1)*i]   = 1.0;
     637        1299 :     pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
     638        1299 :     pH[(nmax+1)*i]   = s[i];
     639        1299 :     pH[(nmax+1)*i+1] = beta[i+1];
     640             :   }
     641          16 :   pH[nmax*nmax-1] = s[nmax-1];
     642          16 :   pK[nmax*nmax-1] = 1.0;
     643          16 :   PetscCall(PetscBLASIntCast(nmax,&n_));
     644          16 :   PetscCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
     645             :   /* The matrix to be used is in H. K will be a work-space matrix */
     646          16 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H));
     647          16 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K));
     648          58 :   for (j=0;matrix&&j<nep->nt;j++) {
     649          42 :     PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler,NULL));
     650          42 :     ierr = FNEvaluateFunctionMat(nep->f[j],H,K);
     651          42 :     PetscCall(PetscPopErrorHandler());
     652          42 :     if (!ierr) {
     653        3587 :       for (i=0;i<nmax;i++) ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0];
     654             :     } else {
     655           0 :       matrix = PETSC_FALSE;
     656          42 :       PetscCall(PetscFPTrapPop());
     657             :     }
     658             :   }
     659          16 :   PetscCall(MatDestroy(&H));
     660          16 :   PetscCall(MatDestroy(&K));
     661          16 :   if (!matrix) {
     662           0 :     for (j=0;j<nep->nt;j++) {
     663           0 :       PetscCall(FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j));
     664           0 :       ctx->coeffD[j] *= beta[0];
     665             :     }
     666             :   }
     667          16 :   if (hasmnorm) {
     668             :     norm0 = 0.0;
     669          49 :     for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
     670             :   } else {
     671             :     norm0 = 0.0;
     672           9 :     for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
     673             :   }
     674          16 :   ctx->nmat = ctx->ddmaxit;
     675         184 :   for (k=1;k<ctx->ddmaxit;k++) {
     676         182 :     if (!matrix) {
     677           0 :       PetscCall(NEPNLEIGSEvalNRTFunct(nep,k,s[k],b));
     678           0 :       for (i=0;i<nep->nt;i++) {
     679           0 :         PetscCall(FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i));
     680           0 :         for (j=0;j<k;j++) {
     681           0 :           ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
     682             :         }
     683           0 :         ctx->coeffD[k*nep->nt+i] /= b[k];
     684             :       }
     685             :     }
     686         182 :     if (hasmnorm) {
     687             :       norm = 0.0;
     688         422 :       for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
     689             :     } else {
     690             :       norm = 0.0;
     691         180 :       for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
     692             :     }
     693         182 :     if (k>1 && norm/norm0 < ctx->ddtol) {
     694          14 :       ctx->nmat = k+1;
     695          14 :       break;
     696             :     }
     697             :   }
     698          16 :   if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
     699          16 :   PetscCall(MatIsShellAny(nep->A,nep->nt,&shell));
     700          16 :   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
     701          38 :   for (i=0;i<ctx->nshiftsw;i++) {
     702          22 :     PetscCall(NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs));
     703          22 :     if (!shell) PetscCall(MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T));
     704           3 :     else PetscCall(NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat));
     705          22 :     if (nep->P) { /* user-defined preconditioner */
     706           1 :       PetscCall(MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P));
     707          21 :     } else P=T;
     708          22 :     alpha = 0.0;
     709         340 :     for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
     710          22 :     PetscCall(MatScale(T,alpha));
     711          22 :     if (nep->P) PetscCall(MatScale(P,alpha));
     712          54 :     for (k=1;k<nep->nt;k++) {
     713             :       alpha = 0.0;
     714         416 :       for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
     715          32 :       if (shell) PetscCall(NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat));
     716          32 :       PetscCall(MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr));
     717          32 :       if (nep->P) PetscCall(MatAXPY(P,alpha,nep->P[k],nep->mstrp));
     718          32 :       if (shell) PetscCall(MatDestroy(&Ts));
     719             :     }
     720          22 :     PetscCall(NEP_KSPSetOperators(ctx->ksp[i],T,P));
     721          22 :     PetscCall(KSPSetUp(ctx->ksp[i]));
     722          22 :     PetscCall(MatDestroy(&T));
     723          22 :     if (nep->P) PetscCall(MatDestroy(&P));
     724             :   }
     725          16 :   PetscCall(PetscFree3(b,coeffs,matnorm));
     726          16 :   PetscCall(PetscFree2(pK,pH));
     727          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     728             : }
     729             : 
     730          11 : static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
     731             : {
     732          11 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     733          11 :   PetscInt       k,j,i,maxnmat;
     734          11 :   PetscReal      norm0,norm;
     735          11 :   PetscScalar    *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
     736          11 :   Mat            *D=ctx->D,*DP,T,P;
     737          11 :   PetscBool      shell,has,vec=PETSC_FALSE,precond=(nep->function_pre!=nep->function)?PETSC_TRUE:PETSC_FALSE;
     738          11 :   PetscRandom    rand=NULL;
     739          11 :   Vec            w[2];
     740             : 
     741          11 :   PetscFunctionBegin;
     742          11 :   PetscCall(PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs));
     743          11 :   if (nep->V) PetscCall(BVGetRandomContext(nep->V,&rand));
     744          11 :   T = nep->function;
     745          11 :   P = nep->function_pre;
     746          11 :   PetscCall(NEPComputeFunction(nep,s[0],T,P));
     747          11 :   PetscCall(MatIsShell(T,&shell));
     748          11 :   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
     749          11 :   if (!shell) PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&D[0]));
     750           3 :   else PetscCall(NLEIGSMatToMatShellArray(T,&D[0],maxnmat));
     751          11 :   if (beta[0]!=1.0) PetscCall(MatScale(D[0],1.0/beta[0]));
     752          11 :   PetscCall(MatHasOperation(D[0],MATOP_NORM,&has));
     753          11 :   if (has) PetscCall(MatNorm(D[0],NORM_FROBENIUS,&norm0));
     754             :   else {
     755           3 :     PetscCall(MatCreateVecs(D[0],NULL,&w[0]));
     756           3 :     PetscCall(VecDuplicate(w[0],&w[1]));
     757           3 :     PetscCall(VecDuplicate(w[0],&ctx->vrn));
     758           3 :     PetscCall(VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]));
     759           3 :     PetscCall(VecNormalize(ctx->vrn,NULL));
     760           3 :     vec = PETSC_TRUE;
     761           3 :     PetscCall(MatNormEstimate(D[0],ctx->vrn,w[0],&norm0));
     762             :   }
     763          11 :   if (precond) {
     764           1 :     PetscCall(PetscMalloc1(ctx->ddmaxit,&DP));
     765           1 :     PetscCall(MatDuplicate(P,MAT_COPY_VALUES,&DP[0]));
     766             :   }
     767          11 :   ctx->nmat = ctx->ddmaxit;
     768         122 :   for (k=1;k<ctx->ddmaxit;k++) {
     769         121 :     PetscCall(NEPNLEIGSEvalNRTFunct(nep,k,s[k],b));
     770         121 :     PetscCall(NEPComputeFunction(nep,s[k],T,P));
     771         121 :     if (!shell) PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&D[k]));
     772          63 :     else PetscCall(NLEIGSMatToMatShellArray(T,&D[k],maxnmat));
     773        1271 :     for (j=0;j<k;j++) PetscCall(MatAXPY(D[k],-b[j],D[j],nep->mstr));
     774         121 :     PetscCall(MatScale(D[k],1.0/b[k]));
     775         121 :     PetscCall(MatHasOperation(D[k],MATOP_NORM,&has));
     776         121 :     if (has) PetscCall(MatNorm(D[k],NORM_FROBENIUS,&norm));
     777             :     else {
     778          63 :       if (!vec) {
     779           0 :         PetscCall(MatCreateVecs(D[k],NULL,&w[0]));
     780           0 :         PetscCall(VecDuplicate(w[0],&w[1]));
     781           0 :         PetscCall(VecDuplicate(w[0],&ctx->vrn));
     782           0 :         PetscCall(VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]));
     783           0 :         PetscCall(VecNormalize(ctx->vrn,NULL));
     784             :         vec = PETSC_TRUE;
     785             :       }
     786          63 :       PetscCall(MatNormEstimate(D[k],ctx->vrn,w[0],&norm));
     787             :     }
     788         121 :     if (precond) {
     789           3 :       PetscCall(MatDuplicate(P,MAT_COPY_VALUES,&DP[k]));
     790           9 :       for (j=0;j<k;j++) PetscCall(MatAXPY(DP[k],-b[j],DP[j],nep->mstrp));
     791           3 :       PetscCall(MatScale(DP[k],1.0/b[k]));
     792             :     }
     793         121 :     if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
     794          10 :       ctx->nmat = k+1;
     795          10 :       break;
     796             :     }
     797             :   }
     798          11 :   if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
     799          25 :   for (i=0;i<ctx->nshiftsw;i++) {
     800          14 :     PetscCall(NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs));
     801          14 :     PetscCall(MatDuplicate(D[0],MAT_COPY_VALUES,&T));
     802          14 :     if (coeffs[0]!=1.0) PetscCall(MatScale(T,coeffs[0]));
     803         192 :     for (j=1;j<ctx->nmat;j++) PetscCall(MatAXPY(T,coeffs[j],D[j],nep->mstr));
     804          14 :     if (precond) {
     805           1 :       PetscCall(MatDuplicate(DP[0],MAT_COPY_VALUES,&P));
     806           1 :       if (coeffs[0]!=1.0) PetscCall(MatScale(P,coeffs[0]));
     807           4 :       for (j=1;j<ctx->nmat;j++) PetscCall(MatAXPY(P,coeffs[j],DP[j],nep->mstrp));
     808          13 :     } else P=T;
     809          14 :     PetscCall(NEP_KSPSetOperators(ctx->ksp[i],T,P));
     810          14 :     PetscCall(KSPSetUp(ctx->ksp[i]));
     811          14 :     PetscCall(MatDestroy(&T));
     812             :   }
     813          11 :   PetscCall(PetscFree2(b,coeffs));
     814          11 :   if (vec) {
     815           3 :     PetscCall(VecDestroy(&w[0]));
     816           3 :     PetscCall(VecDestroy(&w[1]));
     817             :   }
     818          11 :   if (precond) {
     819           1 :     PetscCall(MatDestroy(&P));
     820           1 :     PetscCall(MatDestroyMatrices(ctx->nmat,&DP));
     821             :   }
     822          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     823             : }
     824             : 
     825             : /*
     826             :    NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
     827             : */
     828          57 : static PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
     829             : {
     830          57 :   PetscInt       k,newk,marker,inside;
     831          57 :   PetscScalar    re,im;
     832          57 :   PetscReal      resnorm,tt;
     833          57 :   PetscBool      istrivial;
     834          57 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
     835             : 
     836          57 :   PetscFunctionBegin;
     837          57 :   PetscCall(RGIsTrivial(nep->rg,&istrivial));
     838          57 :   marker = -1;
     839          57 :   if (nep->trackall) getall = PETSC_TRUE;
     840         184 :   for (k=kini;k<kini+nits;k++) {
     841             :     /* eigenvalue */
     842         184 :     re = nep->eigr[k];
     843         184 :     im = nep->eigi[k];
     844         184 :     if (!istrivial) {
     845         184 :       if (!ctx->nshifts) PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im));
     846         184 :       PetscCall(RGCheckInside(nep->rg,1,&re,&im,&inside));
     847         184 :       if (marker==-1 && inside<0) marker = k;
     848             :     }
     849         184 :     newk = k;
     850         184 :     PetscCall(DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm));
     851         184 :     tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
     852         184 :     resnorm *=  PetscAbsReal(tt);
     853             :     /* error estimate */
     854         184 :     PetscCall((*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx));
     855         184 :     if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
     856         184 :     if (newk==k+1) {
     857           0 :       nep->errest[k+1] = nep->errest[k];
     858           0 :       k++;
     859             :     }
     860         184 :     if (marker!=-1 && !getall) break;
     861             :   }
     862          57 :   if (marker!=-1) k = marker;
     863          57 :   *kout = k;
     864          57 :   PetscFunctionReturn(PETSC_SUCCESS);
     865             : }
     866             : 
     867          27 : static PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
     868             : {
     869          27 :   PetscInt       k,in;
     870          27 :   PetscScalar    zero=0.0;
     871          27 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     872          27 :   SlepcSC        sc;
     873          27 :   PetscBool      istrivial;
     874             : 
     875          27 :   PetscFunctionBegin;
     876          27 :   PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
     877          27 :   PetscCheck(nep->ncv<=nep->nev+nep->mpd,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
     878          27 :   if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
     879          27 :   if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
     880          27 :   PetscCall(RGIsTrivial(nep->rg,&istrivial));
     881          27 :   PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
     882          27 :   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
     883          27 :   PetscCheck(nep->which==NEP_TARGET_MAGNITUDE || nep->which==NEP_TARGET_REAL || nep->which==NEP_TARGET_IMAGINARY || nep->which==NEP_WHICH_USER,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target selection of eigenvalues");
     884             : 
     885             :   /* Initialize the NLEIGS context structure */
     886          27 :   k = ctx->ddmaxit;
     887          27 :   PetscCall(PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D));
     888          27 :   nep->data = ctx;
     889          27 :   if (nep->tol==(PetscReal)PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
     890          27 :   if (ctx->ddtol==(PetscReal)PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
     891          27 :   if (!ctx->keep) ctx->keep = 0.5;
     892             : 
     893             :   /* Compute Leja-Bagby points and scaling values */
     894          27 :   PetscCall(NEPNLEIGSLejaBagbyPoints(nep));
     895          27 :   if (nep->problem_type!=NEP_RATIONAL) {
     896          24 :     PetscCall(RGCheckInside(nep->rg,1,&nep->target,&zero,&in));
     897          24 :     PetscCheck(in>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
     898             :   }
     899             : 
     900             :   /* Compute the divided difference matrices */
     901          27 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(NEPNLEIGSDividedDifferences_split(nep));
     902          11 :   else PetscCall(NEPNLEIGSDividedDifferences_callback(nep));
     903          27 :   PetscCall(NEPAllocateSolution(nep,ctx->nmat-1));
     904          27 :   PetscCall(NEPSetWorkVecs(nep,4));
     905          27 :   if (!ctx->fullbasis) {
     906          24 :     PetscCheck(!nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
     907             :     /* set-up DS and transfer split operator functions */
     908          45 :     PetscCall(DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP));
     909          24 :     PetscCall(DSAllocate(nep->ds,nep->ncv+1));
     910          24 :     PetscCall(DSGetSlepcSC(nep->ds,&sc));
     911          24 :     if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
     912          24 :     PetscCall(DSSetExtraRow(nep->ds,PETSC_TRUE));
     913          24 :     sc->mapobj        = (PetscObject)nep;
     914          24 :     sc->rg            = nep->rg;
     915          24 :     sc->comparison    = nep->sc->comparison;
     916          24 :     sc->comparisonctx = nep->sc->comparisonctx;
     917          24 :     PetscCall(BVDestroy(&ctx->V));
     918          24 :     PetscCall(BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V));
     919          24 :     nep->ops->solve          = NEPSolve_NLEIGS;
     920          24 :     nep->ops->computevectors = NEPComputeVectors_Schur;
     921             :   } else {
     922           3 :     PetscCall(NEPSetUp_NLEIGS_FullBasis(nep));
     923           3 :     nep->ops->solve          = NEPSolve_NLEIGS_FullBasis;
     924           3 :     nep->ops->computevectors = NULL;
     925             :   }
     926          27 :   PetscFunctionReturn(PETSC_SUCCESS);
     927             : }
     928             : 
     929             : /*
     930             :   Extend the TOAR basis by applying the matrix operator
     931             :   over a vector which is decomposed on the TOAR way
     932             :   Input:
     933             :     - S,V: define the latest Arnoldi vector (nv vectors in V)
     934             :   Output:
     935             :     - t: new vector extending the TOAR basis
     936             :     - r: temporally coefficients to compute the TOAR coefficients
     937             :          for the new Arnoldi vector
     938             :   Workspace: t_ (two vectors)
     939             : */
     940         931 : static PetscErrorCode NEPTOARExtendBasis(NEP nep,PetscInt idxrktg,PetscScalar *S,PetscInt ls,PetscInt nv,BV W,BV V,Vec t,PetscScalar *r,PetscInt lr,Vec *t_)
     941             : {
     942         931 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     943         931 :   PetscInt       deg=ctx->nmat-1,k,j;
     944         931 :   Vec            v=t_[0],q=t_[1],w;
     945         931 :   PetscScalar    *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
     946             : 
     947         931 :   PetscFunctionBegin;
     948         931 :   if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
     949         931 :   sigma = ctx->shifts[idxrktg];
     950         931 :   PetscCall(BVSetActiveColumns(nep->V,0,nv));
     951         931 :   PetscCall(PetscMalloc1(ctx->nmat,&coeffs));
     952         931 :   PetscCheck(PetscAbsScalar(s[deg-2]-sigma)>100*PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
     953             :   /* i-part stored in (i-1) position */
     954       27812 :   for (j=0;j<nv;j++) {
     955       26881 :     r[(deg-2)*lr+j] = (S[(deg-2)*ls+j]+(beta[deg-1]/xi[deg-2])*S[(deg-1)*ls+j])/(s[deg-2]-sigma);
     956             :   }
     957         931 :   PetscCall(BVSetActiveColumns(W,0,deg));
     958         931 :   PetscCall(BVGetColumn(W,deg-1,&w));
     959         931 :   PetscCall(BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls));
     960         931 :   PetscCall(BVRestoreColumn(W,deg-1,&w));
     961         931 :   PetscCall(BVGetColumn(W,deg-2,&w));
     962         931 :   PetscCall(BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr));
     963         931 :   PetscCall(BVRestoreColumn(W,deg-2,&w));
     964        9348 :   for (k=deg-2;k>0;k--) {
     965        8417 :     PetscCheck(PetscAbsScalar(s[k-1]-sigma)>100*PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
     966      272140 :     for (j=0;j<nv;j++) r[(k-1)*lr+j] = (S[(k-1)*ls+j]+(beta[k]/xi[k-1])*S[k*ls+j]-beta[k]*(1.0-sigma/xi[k-1])*r[(k)*lr+j])/(s[k-1]-sigma);
     967        8417 :     PetscCall(BVGetColumn(W,k-1,&w));
     968        8417 :     PetscCall(BVMultVec(V,1.0,0.0,w,r+(k-1)*lr));
     969        8417 :     PetscCall(BVRestoreColumn(W,k-1,&w));
     970             :   }
     971         931 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     972        7391 :     for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
     973         731 :     coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
     974         731 :     PetscCall(BVMultVec(W,1.0,0.0,v,coeffs));
     975         731 :     PetscCall(MatMult(nep->A[0],v,q));
     976        2090 :     for (k=1;k<nep->nt;k++) {
     977       12030 :       for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
     978        1359 :       coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
     979        1359 :       PetscCall(BVMultVec(W,1.0,0,v,coeffs));
     980        1359 :       PetscCall(MatMult(nep->A[k],v,t));
     981        1359 :       PetscCall(VecAXPY(q,1.0,t));
     982             :     }
     983         731 :     PetscCall(KSPSolve(ctx->ksp[idxrktg],q,t));
     984         731 :     PetscCall(VecScale(t,-1.0));
     985             :   } else {
     986        2888 :     for (k=0;k<deg-1;k++) {
     987        2688 :       PetscCall(BVGetColumn(W,k,&w));
     988        2688 :       PetscCall(MatMult(ctx->D[k],w,q));
     989        2688 :       PetscCall(BVRestoreColumn(W,k,&w));
     990        2688 :       PetscCall(BVInsertVec(W,k,q));
     991             :     }
     992         200 :     PetscCall(BVGetColumn(W,deg-1,&w));
     993         200 :     PetscCall(MatMult(ctx->D[deg],w,q));
     994         200 :     PetscCall(BVRestoreColumn(W,k,&w));
     995         200 :     PetscCall(BVInsertVec(W,k,q));
     996        3088 :     for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
     997         200 :     PetscCall(BVMultVec(W,1.0,0.0,q,coeffs));
     998         200 :     PetscCall(KSPSolve(ctx->ksp[idxrktg],q,t));
     999         200 :     PetscCall(VecScale(t,-1.0));
    1000             :   }
    1001         931 :   PetscCall(PetscFree(coeffs));
    1002         931 :   PetscFunctionReturn(PETSC_SUCCESS);
    1003             : }
    1004             : 
    1005             : /*
    1006             :   Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
    1007             : */
    1008         931 : static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
    1009             : {
    1010         931 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
    1011         931 :   PetscInt       k,j,d=ctx->nmat-1;
    1012         931 :   PetscScalar    *t=work;
    1013             : 
    1014         931 :   PetscFunctionBegin;
    1015         931 :   PetscCall(NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t));
    1016       10279 :   for (k=0;k<d-1;k++) {
    1017      308260 :     for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
    1018             :   }
    1019       28691 :   for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
    1020         931 :   PetscFunctionReturn(PETSC_SUCCESS);
    1021             : }
    1022             : 
    1023             : /*
    1024             :   Compute continuation vector coefficients for the Rational-Krylov run.
    1025             :   dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
    1026             : */
    1027         931 : static PetscErrorCode NEPNLEIGS_RKcontinuation(NEP nep,PetscInt ini,PetscInt end,PetscScalar *K,PetscScalar *H,PetscInt ld,PetscScalar sigma,PetscScalar *S,PetscInt lds,PetscScalar *cont,PetscScalar *t,PetscScalar *work)
    1028             : {
    1029         931 :   PetscScalar    *x,*W,*tau,sone=1.0,szero=0.0;
    1030         931 :   PetscInt       i,j,n1,n,nwu=0;
    1031         931 :   PetscBLASInt   info,n_,n1_,one=1,dim,lds_;
    1032         931 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
    1033             : 
    1034         931 :   PetscFunctionBegin;
    1035         931 :   if (!ctx->nshifts || !end) {
    1036         856 :     t[0] = 1;
    1037         856 :     PetscCall(PetscArraycpy(cont,S+end*lds,lds));
    1038             :   } else {
    1039          75 :     n   = end-ini;
    1040          75 :     n1  = n+1;
    1041          75 :     x   = work+nwu;
    1042          75 :     nwu += end+1;
    1043          75 :     tau = work+nwu;
    1044          75 :     nwu += n;
    1045          75 :     W   = work+nwu;
    1046          75 :     nwu += n1*n;
    1047         858 :     for (j=ini;j<end;j++) {
    1048       11421 :       for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
    1049             :     }
    1050          75 :     PetscCall(PetscBLASIntCast(n,&n_));
    1051          75 :     PetscCall(PetscBLASIntCast(n1,&n1_));
    1052          75 :     PetscCall(PetscBLASIntCast(end+1,&dim));
    1053          75 :     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1054          75 :     PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
    1055          75 :     SlepcCheckLapackInfo("geqrf",info);
    1056         858 :     for (i=0;i<end;i++) t[i] = 0.0;
    1057          75 :     t[end] = 1.0;
    1058         858 :     for (j=n-1;j>=0;j--) {
    1059        5319 :       for (i=0;i<ini+j;i++) x[i] = 0.0;
    1060         783 :       x[ini+j] = 1.0;
    1061        6102 :       for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
    1062         783 :       tau[j] = PetscConj(tau[j]);
    1063         783 :       PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
    1064             :     }
    1065          75 :     PetscCall(PetscBLASIntCast(lds,&lds_));
    1066          75 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
    1067          75 :     PetscCall(PetscFPTrapPop());
    1068             :   }
    1069         931 :   PetscFunctionReturn(PETSC_SUCCESS);
    1070             : }
    1071             : 
    1072             : /*
    1073             :   Compute a run of Arnoldi iterations
    1074             : */
    1075          57 : static PetscErrorCode NEPNLEIGSTOARrun(NEP nep,Mat MK,Mat MH,BV W,PetscInt k,PetscInt *M,PetscReal *betah,PetscScalar *betak,PetscBool *breakdown,Vec *t_)
    1076             : {
    1077          57 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
    1078          57 :   PetscInt       i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l,ldh;
    1079          57 :   Vec            t;
    1080          57 :   PetscReal      norm=0.0;
    1081          57 :   PetscScalar    *x,*work,*tt,sigma=1.0,*cont,*S,*K=NULL,*H;
    1082          57 :   PetscBool      lindep;
    1083          57 :   Mat            MS;
    1084             : 
    1085          57 :   PetscFunctionBegin;
    1086          57 :   *betah = 0.0; *betak = 0.0;
    1087          57 :   PetscCall(MatDenseGetArray(MH,&H));
    1088          57 :   if (MK) PetscCall(MatDenseGetArray(MK,&K));
    1089          57 :   PetscCall(MatDenseGetLDA(MH,&ldh));
    1090          57 :   PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1091          57 :   PetscCall(MatDenseGetArray(MS,&S));
    1092          57 :   PetscCall(BVGetSizes(nep->V,NULL,NULL,&ld));
    1093          57 :   lds = ld*deg;
    1094          57 :   PetscCall(BVGetActiveColumns(nep->V,&l,&nqt));
    1095          57 :   lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
    1096          57 :   PetscCall(PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont));
    1097          57 :   PetscCall(BVSetActiveColumns(ctx->V,0,m));
    1098         988 :   for (j=k;j<m;j++) {
    1099         931 :     sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];
    1100             : 
    1101             :     /* Continuation vector */
    1102         931 :     PetscCall(NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work));
    1103             : 
    1104             :     /* apply operator */
    1105         931 :     PetscCall(BVGetColumn(nep->V,nqt,&t));
    1106         931 :     PetscCall(NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_));
    1107         931 :     PetscCall(BVRestoreColumn(nep->V,nqt,&t));
    1108             : 
    1109             :     /* orthogonalize */
    1110         931 :     PetscCall(BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep));
    1111         931 :     if (!lindep) {
    1112         879 :       x[nqt] = norm;
    1113         879 :       PetscCall(BVScaleColumn(nep->V,nqt,1.0/norm));
    1114         879 :       nqt++;
    1115          52 :     } else x[nqt] = 0.0;
    1116             : 
    1117         931 :     PetscCall(NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work));
    1118             : 
    1119             :     /* Level-2 orthogonalization */
    1120         931 :     PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown));
    1121         931 :     H[j+1+ldh*j] = norm;
    1122         931 :     if (ctx->nshifts && MK) {
    1123         939 :       for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
    1124          78 :       K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
    1125             :     }
    1126         931 :     if (*breakdown) {
    1127           0 :       *M = j+1;
    1128           0 :       break;
    1129             :     }
    1130         931 :     PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
    1131         931 :     PetscCall(BVSetActiveColumns(nep->V,l,nqt));
    1132             :   }
    1133          57 :   *betah = norm;
    1134          57 :   if (ctx->nshifts) *betak = norm*sigma;
    1135          57 :   PetscCall(PetscFree4(x,work,tt,cont));
    1136          57 :   PetscCall(MatDenseRestoreArray(MS,&S));
    1137          57 :   PetscCall(MatDenseRestoreArray(MH,&H));
    1138          57 :   if (MK) PetscCall(MatDenseRestoreArray(MK,&K));
    1139          57 :   PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1140          57 :   PetscFunctionReturn(PETSC_SUCCESS);
    1141             : }
    1142             : 
    1143          24 : PetscErrorCode NEPSolve_NLEIGS(NEP nep)
    1144             : {
    1145          24 :   NEP_NLEIGS        *ctx = (NEP_NLEIGS*)nep->data;
    1146          24 :   PetscInt          i,k=0,l,nv=0,ld,lds,nq;
    1147          24 :   PetscInt          deg=ctx->nmat-1,nconv=0,dsn,dsk;
    1148          24 :   PetscScalar       *pU,betak=0,*eigr,*eigi;
    1149          24 :   const PetscScalar *S;
    1150          24 :   PetscReal         betah;
    1151          24 :   PetscBool         falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
    1152          24 :   BV                W;
    1153          24 :   Mat               H,K=NULL,MS,MQ,U;
    1154             : 
    1155          24 :   PetscFunctionBegin;
    1156          24 :   if (ctx->lock) {
    1157             :     /* undocumented option to use a cheaper locking instead of the true locking */
    1158          21 :     PetscCall(PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL));
    1159             :   }
    1160             : 
    1161          24 :   PetscCall(BVGetSizes(nep->V,NULL,NULL,&ld));
    1162          24 :   lds = deg*ld;
    1163          24 :   if (!ctx->nshifts) PetscCall(PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi));
    1164           3 :   else { eigr = nep->eigr; eigi = nep->eigi; }
    1165          24 :   PetscCall(BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W));
    1166             : 
    1167             :   /* clean projected matrix (including the extra-arrow) */
    1168          24 :   PetscCall(DSSetDimensions(nep->ds,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
    1169          24 :   PetscCall(DSGetMat(nep->ds,DS_MAT_A,&H));
    1170          24 :   PetscCall(MatZeroEntries(H));
    1171          24 :   PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&H));
    1172          24 :   if (ctx->nshifts) {
    1173           3 :     PetscCall(DSGetMat(nep->ds,DS_MAT_B,&H));
    1174           3 :     PetscCall(MatZeroEntries(H));
    1175           3 :     PetscCall(DSRestoreMat(nep->ds,DS_MAT_B,&H));
    1176             :   }
    1177             : 
    1178             :   /* Get the starting Arnoldi vector */
    1179          24 :   PetscCall(BVTensorBuildFirstColumn(ctx->V,nep->nini));
    1180             : 
    1181             :   /* Restart loop */
    1182          24 :   l = 0;
    1183          24 :   while (nep->reason == NEP_CONVERGED_ITERATING) {
    1184          57 :     nep->its++;
    1185             : 
    1186             :     /* Compute an nv-step Krylov relation */
    1187          57 :     nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
    1188          57 :     if (ctx->nshifts) PetscCall(DSGetMat(nep->ds,DS_MAT_A,&K));
    1189          57 :     PetscCall(DSGetMat(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H));
    1190          57 :     PetscCall(NEPNLEIGSTOARrun(nep,K,H,W,nep->nconv+l,&nv,&betah,&betak,&breakdown,nep->work));
    1191          57 :     PetscCall(DSRestoreMat(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H));
    1192          57 :     if (ctx->nshifts) PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&K));
    1193          57 :     PetscCall(DSSetDimensions(nep->ds,nv,nep->nconv,nep->nconv+l));
    1194          57 :     if (l==0) PetscCall(DSSetState(nep->ds,DS_STATE_INTERMEDIATE));
    1195          33 :     else PetscCall(DSSetState(nep->ds,DS_STATE_RAW));
    1196             : 
    1197             :     /* Solve projected problem */
    1198          57 :     PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
    1199          57 :     PetscCall(DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL));
    1200          57 :     PetscCall(DSUpdateExtraRow(nep->ds));
    1201          57 :     PetscCall(DSSynchronize(nep->ds,nep->eigr,nep->eigi));
    1202             : 
    1203             :     /* Check convergence */
    1204          57 :     PetscCall(NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work));
    1205          57 :     PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx));
    1206             : 
    1207             :     /* Update l */
    1208          57 :     if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
    1209             :     else {
    1210          33 :       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
    1211          33 :       PetscCall(DSGetTruncateSize(nep->ds,k,nv,&l));
    1212          33 :       if (!breakdown) {
    1213             :         /* Prepare the Rayleigh quotient for restart */
    1214          33 :         PetscCall(DSGetDimensions(nep->ds,&dsn,NULL,&dsk,NULL));
    1215          33 :         PetscCall(DSSetDimensions(nep->ds,dsn,k,dsk));
    1216          33 :         PetscCall(DSTruncate(nep->ds,k+l,PETSC_FALSE));
    1217             :       }
    1218             :     }
    1219          57 :     nconv = k;
    1220          57 :     if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
    1221          57 :     if (l) PetscCall(PetscInfo(nep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
    1222             : 
    1223             :     /* Update S */
    1224         108 :     PetscCall(DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ));
    1225          57 :     PetscCall(BVMultInPlace(ctx->V,MQ,nep->nconv,k+l));
    1226         108 :     PetscCall(DSRestoreMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ));
    1227             : 
    1228             :     /* Copy last column of S */
    1229          57 :     PetscCall(BVCopyColumn(ctx->V,nv,k+l));
    1230             : 
    1231          57 :     if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
    1232             :       /* Stop if breakdown */
    1233           0 :       PetscCall(PetscInfo(nep,"Breakdown (it=%" PetscInt_FMT " norm=%g)\n",nep->its,(double)betah));
    1234           0 :       nep->reason = NEP_DIVERGED_BREAKDOWN;
    1235             :     }
    1236          57 :     if (nep->reason != NEP_CONVERGED_ITERATING) l--;
    1237             :     /* truncate S */
    1238          57 :     PetscCall(BVGetActiveColumns(nep->V,NULL,&nq));
    1239          57 :     if (k+l+deg<=nq) {
    1240          52 :       PetscCall(BVSetActiveColumns(ctx->V,nep->nconv,k+l+1));
    1241          52 :       if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-nep->nconv));
    1242           6 :       else PetscCall(BVTensorCompress(ctx->V,0));
    1243             :     }
    1244          57 :     nep->nconv = k;
    1245          57 :     if (!ctx->nshifts) {
    1246        1440 :       for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
    1247          51 :       PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi));
    1248             :     }
    1249          81 :     PetscCall(NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv));
    1250             :   }
    1251          24 :   nep->nconv = nconv;
    1252          24 :   if (nep->nconv>0) {
    1253          24 :     PetscCall(BVSetActiveColumns(ctx->V,0,nep->nconv));
    1254          24 :     PetscCall(BVGetActiveColumns(nep->V,NULL,&nq));
    1255          24 :     PetscCall(BVSetActiveColumns(nep->V,0,nq));
    1256          24 :     if (nq>nep->nconv) {
    1257           4 :       PetscCall(BVTensorCompress(ctx->V,nep->nconv));
    1258           4 :       PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
    1259           4 :       nq = nep->nconv;
    1260             :     }
    1261          24 :     if (ctx->nshifts) {
    1262           3 :       PetscCall(DSGetMat(nep->ds,DS_MAT_B,&MQ));
    1263           3 :       PetscCall(BVMultInPlace(ctx->V,MQ,0,nep->nconv));
    1264           3 :       PetscCall(DSRestoreMat(nep->ds,DS_MAT_B,&MQ));
    1265             :     }
    1266          24 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1267          24 :     PetscCall(MatDenseGetArrayRead(MS,&S));
    1268          24 :     PetscCall(PetscMalloc1(nq*nep->nconv,&pU));
    1269         145 :     for (i=0;i<nep->nconv;i++) PetscCall(PetscArraycpy(pU+i*nq,S+i*lds,nq));
    1270          24 :     PetscCall(MatDenseRestoreArrayRead(MS,&S));
    1271          24 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1272          24 :     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U));
    1273          24 :     PetscCall(BVSetActiveColumns(nep->V,0,nq));
    1274          24 :     PetscCall(BVMultInPlace(nep->V,U,0,nep->nconv));
    1275          24 :     PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
    1276          24 :     PetscCall(MatDestroy(&U));
    1277          24 :     PetscCall(PetscFree(pU));
    1278          24 :     PetscCall(DSTruncate(nep->ds,nep->nconv,PETSC_TRUE));
    1279             :   }
    1280             : 
    1281             :   /* Map eigenvalues back to the original problem */
    1282          24 :   if (!ctx->nshifts) {
    1283          21 :     PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi));
    1284          21 :     PetscCall(PetscFree2(eigr,eigi));
    1285             :   }
    1286          24 :   PetscCall(BVDestroy(&W));
    1287          24 :   PetscFunctionReturn(PETSC_SUCCESS);
    1288             : }
    1289             : 
    1290          27 : static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
    1291             : {
    1292          27 :   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
    1293             : 
    1294          27 :   PetscFunctionBegin;
    1295          27 :   if (fun) nepctx->computesingularities = fun;
    1296          27 :   if (ctx) nepctx->singularitiesctx     = ctx;
    1297          27 :   nep->state = NEP_STATE_INITIAL;
    1298          27 :   PetscFunctionReturn(PETSC_SUCCESS);
    1299             : }
    1300             : 
    1301             : /*@C
    1302             :    NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
    1303             :    of the singularity set (where T(.) is not analytic).
    1304             : 
    1305             :    Logically Collective
    1306             : 
    1307             :    Input Parameters:
    1308             : +  nep - the NEP context
    1309             : .  fun - user function (if NULL then NEP retains any previously set value)
    1310             : -  ctx - [optional] user-defined context for private data for the function
    1311             :          (may be NULL, in which case NEP retains any previously set value)
    1312             : 
    1313             :    Calling sequence of fun:
    1314             : $  PetscErrorCode fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
    1315             : +   nep   - the NEP context
    1316             : .   maxnp - on input number of requested points in the discretization (can be set)
    1317             : .   xi    - computed values of the discretization
    1318             : -   ctx   - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
    1319             : 
    1320             :    Notes:
    1321             :    The user-defined function can set a smaller value of maxnp if necessary.
    1322             :    It is wrong to return a larger value.
    1323             : 
    1324             :    If the problem type has been set to rational with NEPSetProblemType(),
    1325             :    then it is not necessary to set the singularities explicitly since the
    1326             :    solver will try to determine them automatically.
    1327             : 
    1328             :    Level: intermediate
    1329             : 
    1330             : .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
    1331             : @*/
    1332          27 : PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx),void *ctx)
    1333             : {
    1334          27 :   PetscFunctionBegin;
    1335          27 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1336          27 :   PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
    1337          27 :   PetscFunctionReturn(PETSC_SUCCESS);
    1338             : }
    1339             : 
    1340           1 : static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
    1341             : {
    1342           1 :   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
    1343             : 
    1344           1 :   PetscFunctionBegin;
    1345           1 :   if (fun) *fun = nepctx->computesingularities;
    1346           1 :   if (ctx) *ctx = nepctx->singularitiesctx;
    1347           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1348             : }
    1349             : 
    1350             : /*@C
    1351             :    NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
    1352             :    provided context for computing a discretization of the singularity set.
    1353             : 
    1354             :    Not Collective
    1355             : 
    1356             :    Input Parameter:
    1357             : .  nep - the nonlinear eigensolver context
    1358             : 
    1359             :    Output Parameters:
    1360             : +  fun - location to put the function (or NULL)
    1361             : -  ctx - location to stash the function context (or NULL)
    1362             : 
    1363             :    Level: advanced
    1364             : 
    1365             : .seealso: NEPNLEIGSSetSingularitiesFunction()
    1366             : @*/
    1367           1 : PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
    1368             : {
    1369           1 :   PetscFunctionBegin;
    1370           1 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1371           1 :   PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
    1372           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1373             : }
    1374             : 
    1375           3 : static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
    1376             : {
    1377           3 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1378             : 
    1379           3 :   PetscFunctionBegin;
    1380           3 :   if (keep==(PetscReal)PETSC_DEFAULT) ctx->keep = 0.5;
    1381             :   else {
    1382           3 :     PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
    1383           3 :     ctx->keep = keep;
    1384             :   }
    1385           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1386             : }
    1387             : 
    1388             : /*@
    1389             :    NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
    1390             :    method, in particular the proportion of basis vectors that must be kept
    1391             :    after restart.
    1392             : 
    1393             :    Logically Collective
    1394             : 
    1395             :    Input Parameters:
    1396             : +  nep  - the nonlinear eigensolver context
    1397             : -  keep - the number of vectors to be kept at restart
    1398             : 
    1399             :    Options Database Key:
    1400             : .  -nep_nleigs_restart - Sets the restart parameter
    1401             : 
    1402             :    Notes:
    1403             :    Allowed values are in the range [0.1,0.9]. The default is 0.5.
    1404             : 
    1405             :    Level: advanced
    1406             : 
    1407             : .seealso: NEPNLEIGSGetRestart()
    1408             : @*/
    1409           3 : PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
    1410             : {
    1411           3 :   PetscFunctionBegin;
    1412           3 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1413          12 :   PetscValidLogicalCollectiveReal(nep,keep,2);
    1414           3 :   PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
    1415           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1416             : }
    1417             : 
    1418           6 : static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
    1419             : {
    1420           6 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1421             : 
    1422           6 :   PetscFunctionBegin;
    1423           6 :   *keep = ctx->keep;
    1424           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1425             : }
    1426             : 
    1427             : /*@
    1428             :    NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
    1429             : 
    1430             :    Not Collective
    1431             : 
    1432             :    Input Parameter:
    1433             : .  nep - the nonlinear eigensolver context
    1434             : 
    1435             :    Output Parameter:
    1436             : .  keep - the restart parameter
    1437             : 
    1438             :    Level: advanced
    1439             : 
    1440             : .seealso: NEPNLEIGSSetRestart()
    1441             : @*/
    1442           6 : PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
    1443             : {
    1444           6 :   PetscFunctionBegin;
    1445           6 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1446           6 :   PetscAssertPointer(keep,2);
    1447           6 :   PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
    1448           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1449             : }
    1450             : 
    1451           3 : static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
    1452             : {
    1453           3 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1454             : 
    1455           3 :   PetscFunctionBegin;
    1456           3 :   ctx->lock = lock;
    1457           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1458             : }
    1459             : 
    1460             : /*@
    1461             :    NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
    1462             :    the NLEIGS method.
    1463             : 
    1464             :    Logically Collective
    1465             : 
    1466             :    Input Parameters:
    1467             : +  nep  - the nonlinear eigensolver context
    1468             : -  lock - true if the locking variant must be selected
    1469             : 
    1470             :    Options Database Key:
    1471             : .  -nep_nleigs_locking - Sets the locking flag
    1472             : 
    1473             :    Notes:
    1474             :    The default is to lock converged eigenpairs when the method restarts.
    1475             :    This behaviour can be changed so that all directions are kept in the
    1476             :    working subspace even if already converged to working accuracy (the
    1477             :    non-locking variant).
    1478             : 
    1479             :    Level: advanced
    1480             : 
    1481             : .seealso: NEPNLEIGSGetLocking()
    1482             : @*/
    1483           3 : PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
    1484             : {
    1485           3 :   PetscFunctionBegin;
    1486           3 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1487          12 :   PetscValidLogicalCollectiveBool(nep,lock,2);
    1488           3 :   PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
    1489           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1490             : }
    1491             : 
    1492           6 : static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
    1493             : {
    1494           6 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1495             : 
    1496           6 :   PetscFunctionBegin;
    1497           6 :   *lock = ctx->lock;
    1498           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1499             : }
    1500             : 
    1501             : /*@
    1502             :    NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
    1503             : 
    1504             :    Not Collective
    1505             : 
    1506             :    Input Parameter:
    1507             : .  nep - the nonlinear eigensolver context
    1508             : 
    1509             :    Output Parameter:
    1510             : .  lock - the locking flag
    1511             : 
    1512             :    Level: advanced
    1513             : 
    1514             : .seealso: NEPNLEIGSSetLocking()
    1515             : @*/
    1516           6 : PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
    1517             : {
    1518           6 :   PetscFunctionBegin;
    1519           6 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1520           6 :   PetscAssertPointer(lock,2);
    1521           6 :   PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
    1522           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1523             : }
    1524             : 
    1525           9 : static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
    1526             : {
    1527           9 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
    1528             : 
    1529           9 :   PetscFunctionBegin;
    1530           9 :   if (tol == (PetscReal)PETSC_DEFAULT) {
    1531           6 :     ctx->ddtol = PETSC_DEFAULT;
    1532           6 :     nep->state = NEP_STATE_INITIAL;
    1533             :   } else {
    1534           3 :     PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
    1535           3 :     ctx->ddtol = tol;
    1536             :   }
    1537           9 :   if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
    1538           0 :     ctx->ddmaxit = 0;
    1539           0 :     if (nep->state) PetscCall(NEPReset(nep));
    1540           0 :     nep->state = NEP_STATE_INITIAL;
    1541             :   } else {
    1542           9 :     PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
    1543           9 :     if (ctx->ddmaxit != degree) {
    1544           9 :       ctx->ddmaxit = degree;
    1545           9 :       if (nep->state) PetscCall(NEPReset(nep));
    1546           9 :       nep->state = NEP_STATE_INITIAL;
    1547             :     }
    1548             :   }
    1549           9 :   PetscFunctionReturn(PETSC_SUCCESS);
    1550             : }
    1551             : 
    1552             : /*@
    1553             :    NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
    1554             :    when building the interpolation via divided differences.
    1555             : 
    1556             :    Collective
    1557             : 
    1558             :    Input Parameters:
    1559             : +  nep    - the nonlinear eigensolver context
    1560             : .  tol    - tolerance to stop computing divided differences
    1561             : -  degree - maximum degree of interpolation
    1562             : 
    1563             :    Options Database Key:
    1564             : +  -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
    1565             : -  -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
    1566             : 
    1567             :    Notes:
    1568             :    Use PETSC_DEFAULT for either argument to assign a reasonably good value.
    1569             : 
    1570             :    Level: advanced
    1571             : 
    1572             : .seealso: NEPNLEIGSGetInterpolation()
    1573             : @*/
    1574           9 : PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
    1575             : {
    1576           9 :   PetscFunctionBegin;
    1577           9 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1578          36 :   PetscValidLogicalCollectiveReal(nep,tol,2);
    1579          36 :   PetscValidLogicalCollectiveInt(nep,degree,3);
    1580           9 :   PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
    1581           9 :   PetscFunctionReturn(PETSC_SUCCESS);
    1582             : }
    1583             : 
    1584          31 : static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
    1585             : {
    1586          31 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1587             : 
    1588          31 :   PetscFunctionBegin;
    1589          31 :   if (tol)    *tol    = ctx->ddtol;
    1590          31 :   if (degree) *degree = ctx->ddmaxit;
    1591          31 :   PetscFunctionReturn(PETSC_SUCCESS);
    1592             : }
    1593             : 
    1594             : /*@
    1595             :    NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
    1596             :    when building the interpolation via divided differences.
    1597             : 
    1598             :    Not Collective
    1599             : 
    1600             :    Input Parameter:
    1601             : .  nep - the nonlinear eigensolver context
    1602             : 
    1603             :    Output Parameters:
    1604             : +  tol    - tolerance to stop computing divided differences
    1605             : -  degree - maximum degree of interpolation
    1606             : 
    1607             :    Level: advanced
    1608             : 
    1609             : .seealso: NEPNLEIGSSetInterpolation()
    1610             : @*/
    1611          31 : PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
    1612             : {
    1613          31 :   PetscFunctionBegin;
    1614          31 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1615          31 :   PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
    1616          31 :   PetscFunctionReturn(PETSC_SUCCESS);
    1617             : }
    1618             : 
    1619           3 : static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
    1620             : {
    1621           3 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
    1622           3 :   PetscInt       i;
    1623             : 
    1624           3 :   PetscFunctionBegin;
    1625           3 :   PetscCheck(ns>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
    1626           3 :   if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
    1627           3 :   for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPDestroy(&ctx->ksp[i]));
    1628           3 :   PetscCall(PetscFree(ctx->ksp));
    1629           3 :   ctx->ksp = NULL;
    1630           3 :   if (ns) {
    1631           3 :     PetscCall(PetscMalloc1(ns,&ctx->shifts));
    1632          15 :     for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
    1633             :   }
    1634           3 :   ctx->nshifts = ns;
    1635           3 :   nep->state   = NEP_STATE_INITIAL;
    1636           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1637             : }
    1638             : 
    1639             : /*@
    1640             :    NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
    1641             :    Krylov method.
    1642             : 
    1643             :    Collective
    1644             : 
    1645             :    Input Parameters:
    1646             : +  nep    - the nonlinear eigensolver context
    1647             : .  ns     - number of shifts
    1648             : -  shifts - array of scalar values specifying the shifts
    1649             : 
    1650             :    Options Database Key:
    1651             : .  -nep_nleigs_rk_shifts - Sets the list of shifts
    1652             : 
    1653             :    Notes:
    1654             :    If only one shift is provided, the built subspace built is equivalent to
    1655             :    shift-and-invert Krylov-Schur (provided that the absolute convergence
    1656             :    criterion is used).
    1657             : 
    1658             :    In the case of real scalars, complex shifts are not allowed. In the
    1659             :    command line, a comma-separated list of complex values can be provided with
    1660             :    the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
    1661             :    -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
    1662             : 
    1663             :    Use ns=0 to remove previously set shifts.
    1664             : 
    1665             :    Level: advanced
    1666             : 
    1667             : .seealso: NEPNLEIGSGetRKShifts()
    1668             : @*/
    1669           3 : PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
    1670             : {
    1671           3 :   PetscFunctionBegin;
    1672           3 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1673          12 :   PetscValidLogicalCollectiveInt(nep,ns,2);
    1674           3 :   if (ns) PetscAssertPointer(shifts,3);
    1675           3 :   PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
    1676           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1677             : }
    1678             : 
    1679           1 : static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
    1680             : {
    1681           1 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
    1682           1 :   PetscInt       i;
    1683             : 
    1684           1 :   PetscFunctionBegin;
    1685           1 :   *ns = ctx->nshifts;
    1686           1 :   if (ctx->nshifts) {
    1687           1 :     PetscCall(PetscMalloc1(ctx->nshifts,shifts));
    1688           5 :     for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
    1689             :   }
    1690           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1691             : }
    1692             : 
    1693             : /*@C
    1694             :    NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
    1695             :    Krylov method.
    1696             : 
    1697             :    Not Collective
    1698             : 
    1699             :    Input Parameter:
    1700             : .  nep - the nonlinear eigensolver context
    1701             : 
    1702             :    Output Parameters:
    1703             : +  ns     - number of shifts
    1704             : -  shifts - array of shifts
    1705             : 
    1706             :    Note:
    1707             :    The user is responsible for deallocating the returned array.
    1708             : 
    1709             :    Level: advanced
    1710             : 
    1711             : .seealso: NEPNLEIGSSetRKShifts()
    1712             : @*/
    1713           1 : PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
    1714             : {
    1715           1 :   PetscFunctionBegin;
    1716           1 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1717           1 :   PetscAssertPointer(ns,2);
    1718           1 :   PetscAssertPointer(shifts,3);
    1719           1 :   PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
    1720           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1721             : }
    1722             : 
    1723          25 : static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
    1724             : {
    1725          25 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
    1726          25 :   PetscInt       i;
    1727          25 :   PC             pc;
    1728             : 
    1729          25 :   PetscFunctionBegin;
    1730          25 :   if (!ctx->ksp) {
    1731          25 :     PetscCall(NEPNLEIGSSetShifts(nep,&ctx->nshiftsw));
    1732          25 :     PetscCall(PetscMalloc1(ctx->nshiftsw,&ctx->ksp));
    1733          59 :     for (i=0;i<ctx->nshiftsw;i++) {
    1734          34 :       PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]));
    1735          34 :       PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1));
    1736          34 :       PetscCall(KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix));
    1737          34 :       PetscCall(KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_"));
    1738          34 :       PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options));
    1739          34 :       PetscCall(KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE));
    1740          58 :       PetscCall(KSPSetTolerances(ctx->ksp[i],1e-3*SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
    1741          34 :       PetscCall(KSPGetPC(ctx->ksp[i],&pc));
    1742          34 :       if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
    1743           2 :         PetscCall(KSPSetType(ctx->ksp[i],KSPBCGS));
    1744           2 :         PetscCall(PCSetType(pc,PCBJACOBI));
    1745             :       } else {
    1746          32 :         PetscCall(KSPSetType(ctx->ksp[i],KSPPREONLY));
    1747          34 :         PetscCall(PCSetType(pc,PCLU));
    1748             :       }
    1749             :     }
    1750             :   }
    1751          25 :   if (nsolve) *nsolve = ctx->nshiftsw;
    1752          25 :   if (ksp)    *ksp    = ctx->ksp;
    1753          25 :   PetscFunctionReturn(PETSC_SUCCESS);
    1754             : }
    1755             : 
    1756             : /*@C
    1757             :    NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
    1758             :    the nonlinear eigenvalue solver.
    1759             : 
    1760             :    Collective
    1761             : 
    1762             :    Input Parameter:
    1763             : .  nep - nonlinear eigenvalue solver
    1764             : 
    1765             :    Output Parameters:
    1766             : +  nsolve - number of returned KSP objects
    1767             : -  ksp - array of linear solver object
    1768             : 
    1769             :    Notes:
    1770             :    The number of KSP objects is equal to the number of shifts provided by the user,
    1771             :    or 1 if the user did not provide shifts.
    1772             : 
    1773             :    Level: advanced
    1774             : 
    1775             : .seealso: NEPNLEIGSSetRKShifts()
    1776             : @*/
    1777          25 : PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
    1778             : {
    1779          25 :   PetscFunctionBegin;
    1780          25 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1781          25 :   PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
    1782          25 :   PetscFunctionReturn(PETSC_SUCCESS);
    1783             : }
    1784             : 
    1785           4 : static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
    1786             : {
    1787           4 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1788             : 
    1789           4 :   PetscFunctionBegin;
    1790           4 :   if (fullbasis!=ctx->fullbasis) {
    1791           3 :     ctx->fullbasis = fullbasis;
    1792           3 :     nep->state     = NEP_STATE_INITIAL;
    1793           3 :     nep->useds     = PetscNot(fullbasis);
    1794             :   }
    1795           4 :   PetscFunctionReturn(PETSC_SUCCESS);
    1796             : }
    1797             : 
    1798             : /*@
    1799             :    NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
    1800             :    variants of the NLEIGS method.
    1801             : 
    1802             :    Logically Collective
    1803             : 
    1804             :    Input Parameters:
    1805             : +  nep  - the nonlinear eigensolver context
    1806             : -  fullbasis - true if the full-basis variant must be selected
    1807             : 
    1808             :    Options Database Key:
    1809             : .  -nep_nleigs_full_basis - Sets the full-basis flag
    1810             : 
    1811             :    Notes:
    1812             :    The default is to use a compact representation of the Krylov basis, that is,
    1813             :    V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
    1814             :    the full basis V is explicitly stored and operated with. This variant is more
    1815             :    expensive in terms of memory and computation, but is necessary in some cases,
    1816             :    particularly for two-sided computations, see NEPSetTwoSided().
    1817             : 
    1818             :    In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
    1819             :    solve the linearized eigenproblem, see NEPNLEIGSGetEPS().
    1820             : 
    1821             :    Level: advanced
    1822             : 
    1823             : .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
    1824             : @*/
    1825           4 : PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
    1826             : {
    1827           4 :   PetscFunctionBegin;
    1828           4 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1829          16 :   PetscValidLogicalCollectiveBool(nep,fullbasis,2);
    1830           4 :   PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
    1831           4 :   PetscFunctionReturn(PETSC_SUCCESS);
    1832             : }
    1833             : 
    1834           1 : static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
    1835             : {
    1836           1 :   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
    1837             : 
    1838           1 :   PetscFunctionBegin;
    1839           1 :   *fullbasis = ctx->fullbasis;
    1840           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1841             : }
    1842             : 
    1843             : /*@
    1844             :    NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
    1845             :    full-basis variant.
    1846             : 
    1847             :    Not Collective
    1848             : 
    1849             :    Input Parameter:
    1850             : .  nep - the nonlinear eigensolver context
    1851             : 
    1852             :    Output Parameter:
    1853             : .  fullbasis - the flag
    1854             : 
    1855             :    Level: advanced
    1856             : 
    1857             : .seealso: NEPNLEIGSSetFullBasis()
    1858             : @*/
    1859           1 : PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
    1860             : {
    1861           1 :   PetscFunctionBegin;
    1862           1 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
    1863           1 :   PetscAssertPointer(fullbasis,2);
    1864           1 :   PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
    1865           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1866             : }
    1867             : 
    1868             : #define SHIFTMAX 30
    1869             : 
    1870          25 : static PetscErrorCode NEPSetFromOptions_NLEIGS(NEP nep,PetscOptionItems *PetscOptionsObject)
    1871             : {
    1872          25 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
    1873          25 :   PetscInt       i=0,k;
    1874          25 :   PetscBool      flg1,flg2,b;
    1875          25 :   PetscReal      r;
    1876          25 :   PetscScalar    array[SHIFTMAX];
    1877             : 
    1878          25 :   PetscFunctionBegin;
    1879          25 :   PetscOptionsHeadBegin(PetscOptionsObject,"NEP NLEIGS Options");
    1880             : 
    1881          25 :     PetscCall(PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1));
    1882          25 :     if (flg1) PetscCall(NEPNLEIGSSetRestart(nep,r));
    1883             : 
    1884          25 :     PetscCall(PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1));
    1885          25 :     if (flg1) PetscCall(NEPNLEIGSSetLocking(nep,b));
    1886             : 
    1887          25 :     PetscCall(PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1));
    1888          25 :     if (flg1) PetscCall(NEPNLEIGSSetFullBasis(nep,b));
    1889             : 
    1890          25 :     PetscCall(NEPNLEIGSGetInterpolation(nep,&r,&i));
    1891          25 :     if (!i) i = PETSC_DEFAULT;
    1892          25 :     PetscCall(PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1));
    1893          25 :     PetscCall(PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2));
    1894          25 :     if (flg1 || flg2) PetscCall(NEPNLEIGSSetInterpolation(nep,r,i));
    1895             : 
    1896          25 :     k = SHIFTMAX;
    1897         775 :     for (i=0;i<k;i++) array[i] = 0;
    1898          25 :     PetscCall(PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1));
    1899          25 :     if (flg1) PetscCall(NEPNLEIGSSetRKShifts(nep,k,array));
    1900             : 
    1901          25 :   PetscOptionsHeadEnd();
    1902             : 
    1903          25 :   if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
    1904          59 :   for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPSetFromOptions(ctx->ksp[i]));
    1905             : 
    1906          25 :   if (ctx->fullbasis) {
    1907           3 :     if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
    1908           3 :     PetscCall(EPSSetFromOptions(ctx->eps));
    1909             :   }
    1910          25 :   PetscFunctionReturn(PETSC_SUCCESS);
    1911             : }
    1912             : 
    1913           1 : static PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
    1914             : {
    1915           1 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
    1916           1 :   PetscBool      isascii;
    1917           1 :   PetscInt       i;
    1918           1 :   char           str[50];
    1919             : 
    1920           1 :   PetscFunctionBegin;
    1921           1 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    1922           1 :   if (isascii) {
    1923           1 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
    1924           1 :     if (ctx->fullbasis) PetscCall(PetscViewerASCIIPrintf(viewer,"  using the full-basis variant\n"));
    1925           1 :     else PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
    1926           1 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  divided difference terms: used=%" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->nmat,ctx->ddmaxit));
    1927           1 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  tolerance for divided difference convergence: %g\n",(double)ctx->ddtol));
    1928           1 :     if (ctx->nshifts) {
    1929           1 :       PetscCall(PetscViewerASCIIPrintf(viewer,"  RK shifts: "));
    1930           1 :       PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
    1931           5 :       for (i=0;i<ctx->nshifts;i++) {
    1932           4 :         PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE));
    1933           5 :         PetscCall(PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":""));
    1934             :       }
    1935           1 :       PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
    1936           1 :       PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
    1937             :     }
    1938           1 :     if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
    1939           1 :     PetscCall(PetscViewerASCIIPushTab(viewer));
    1940           1 :     PetscCall(KSPView(ctx->ksp[0],viewer));
    1941           1 :     PetscCall(PetscViewerASCIIPopTab(viewer));
    1942           1 :     if (ctx->fullbasis) {
    1943           0 :       if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
    1944           0 :       PetscCall(PetscViewerASCIIPushTab(viewer));
    1945           0 :       PetscCall(EPSView(ctx->eps,viewer));
    1946           0 :       PetscCall(PetscViewerASCIIPopTab(viewer));
    1947             :     }
    1948             :   }
    1949           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1950             : }
    1951             : 
    1952          27 : static PetscErrorCode NEPReset_NLEIGS(NEP nep)
    1953             : {
    1954          27 :   PetscInt       k;
    1955          27 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
    1956             : 
    1957          27 :   PetscFunctionBegin;
    1958          27 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(PetscFree(ctx->coeffD));
    1959             :   else {
    1960         143 :     for (k=0;k<ctx->nmat;k++) PetscCall(MatDestroy(&ctx->D[k]));
    1961             :   }
    1962          27 :   PetscCall(PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D));
    1963          63 :   for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPReset(ctx->ksp[k]));
    1964          27 :   PetscCall(VecDestroy(&ctx->vrn));
    1965          27 :   if (ctx->fullbasis) {
    1966           3 :     PetscCall(MatDestroy(&ctx->A));
    1967           3 :     PetscCall(EPSReset(ctx->eps));
    1968          15 :     for (k=0;k<4;k++) PetscCall(VecDestroy(&ctx->w[k]));
    1969             :   }
    1970          27 :   PetscFunctionReturn(PETSC_SUCCESS);
    1971             : }
    1972             : 
    1973          40 : static PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
    1974             : {
    1975          40 :   PetscInt       k;
    1976          40 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
    1977             : 
    1978          40 :   PetscFunctionBegin;
    1979          40 :   PetscCall(BVDestroy(&ctx->V));
    1980          74 :   for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPDestroy(&ctx->ksp[k]));
    1981          40 :   PetscCall(PetscFree(ctx->ksp));
    1982          40 :   if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
    1983          40 :   if (ctx->fullbasis) PetscCall(EPSDestroy(&ctx->eps));
    1984          40 :   PetscCall(PetscFree(nep->data));
    1985          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL));
    1986          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL));
    1987          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL));
    1988          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL));
    1989          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL));
    1990          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL));
    1991          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL));
    1992          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL));
    1993          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL));
    1994          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL));
    1995          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL));
    1996          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL));
    1997          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL));
    1998          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL));
    1999          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL));
    2000          40 :   PetscFunctionReturn(PETSC_SUCCESS);
    2001             : }
    2002             : 
    2003          40 : SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
    2004             : {
    2005          40 :   NEP_NLEIGS     *ctx;
    2006             : 
    2007          40 :   PetscFunctionBegin;
    2008          40 :   PetscCall(PetscNew(&ctx));
    2009          40 :   nep->data  = (void*)ctx;
    2010          40 :   ctx->lock  = PETSC_TRUE;
    2011          40 :   ctx->ddtol = PETSC_DEFAULT;
    2012             : 
    2013          40 :   nep->useds = PETSC_TRUE;
    2014             : 
    2015          40 :   nep->ops->setup          = NEPSetUp_NLEIGS;
    2016          40 :   nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
    2017          40 :   nep->ops->view           = NEPView_NLEIGS;
    2018          40 :   nep->ops->destroy        = NEPDestroy_NLEIGS;
    2019          40 :   nep->ops->reset          = NEPReset_NLEIGS;
    2020             : 
    2021          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS));
    2022          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS));
    2023          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS));
    2024          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS));
    2025          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS));
    2026          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS));
    2027          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS));
    2028          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS));
    2029          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS));
    2030          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS));
    2031          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS));
    2032          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS));
    2033          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS));
    2034          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS));
    2035          40 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS));
    2036          40 :   PetscFunctionReturn(PETSC_SUCCESS);
    2037             : }

Generated by: LCOV version 1.14