LCOV - code coverage report
Current view: top level - pep/impls/jd - pjd.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 1215 1247 97.4 %
Date: 2024-03-28 00:52:16 Functions: 44 45 97.8 %
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 polynomial eigensolver: "jd"
      12             : 
      13             :    Method: Jacobi-Davidson
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Jacobi-Davidson for polynomial eigenvalue problems.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] C. Campos and J.E. Roman, "A polynomial Jacobi-Davidson solver
      22             :            with support for non-monomial bases and deflation", BIT Numer.
      23             :            Math. 60:295-318, 2020.
      24             : 
      25             :        [2] G.L.G. Sleijpen et al., "Jacobi-Davidson type methods for
      26             :            generalized eigenproblems and polynomial eigenproblems", BIT
      27             :            36(3):595-633, 1996.
      28             : 
      29             :        [3] Feng-Nan Hwang, Zih-Hao Wei, Tsung-Ming Huang, Weichung Wang,
      30             :            "A Parallel Additive Schwarz Preconditioned Jacobi-Davidson
      31             :            Algorithm for Polynomial Eigenvalue Problems in Quantum Dot
      32             :            Simulation", J. Comput. Phys. 229(8):2932-2947, 2010.
      33             : */
      34             : 
      35             : #include <slepc/private/pepimpl.h>    /*I "slepcpep.h" I*/
      36             : #include <slepcblaslapack.h>
      37             : 
      38             : static PetscBool  cited = PETSC_FALSE;
      39             : static const char citation[] =
      40             :   "@Article{slepc-slice-qep,\n"
      41             :   "   author = \"C. Campos and J. E. Roman\",\n"
      42             :   "   title = \"A polynomial {Jacobi-Davidson} solver with support for non-monomial bases and deflation\",\n"
      43             :   "   journal = \"{BIT} Numer. Math.\",\n"
      44             :   "   volume = \"60\",\n"
      45             :   "   pages = \"295--318\",\n"
      46             :   "   year = \"2020,\"\n"
      47             :   "   doi = \"https://doi.org/10.1007/s10543-019-00778-z\"\n"
      48             :   "}\n";
      49             : 
      50             : typedef struct {
      51             :   PetscReal   keep;          /* restart parameter */
      52             :   PetscReal   fix;           /* fix parameter */
      53             :   PetscBool   reusepc;       /* flag indicating whether pc is rebuilt or not */
      54             :   BV          V;             /* work basis vectors to store the search space */
      55             :   BV          W;             /* work basis vectors to store the test space */
      56             :   BV          *TV;           /* work basis vectors to store T*V (each TV[i] is the coefficient for \lambda^i of T*V for the extended T) */
      57             :   BV          *AX;           /* work basis vectors to store A_i*X for locked eigenvectors */
      58             :   BV          N[2];          /* auxiliary work BVs */
      59             :   BV          X;             /* locked eigenvectors */
      60             :   PetscScalar *T;            /* matrix of the invariant pair */
      61             :   PetscScalar *Tj;           /* matrix containing the powers of the invariant pair matrix */
      62             :   PetscScalar *XpX;          /* X^H*X */
      63             :   PetscInt    ld;            /* leading dimension for Tj and XpX */
      64             :   PC          pcshell;       /* preconditioner including basic precond+projector */
      65             :   Mat         Pshell;        /* auxiliary shell matrix */
      66             :   PetscInt    nlock;         /* number of locked vectors in the invariant pair */
      67             :   Vec         vtempl;        /* reference nested vector */
      68             :   PetscInt    midx;          /* minimality index */
      69             :   PetscInt    mmidx;         /* maximum allowed minimality index */
      70             :   PEPJDProjection proj;      /* projection type (orthogonal, harmonic) */
      71             : } PEP_JD;
      72             : 
      73             : typedef struct {
      74             :   PEP         pep;
      75             :   PC          pc;            /* basic preconditioner */
      76             :   Vec         Bp[2];         /* preconditioned residual of derivative polynomial, B\p */
      77             :   Vec         u[2];          /* Ritz vector */
      78             :   PetscScalar gamma[2];      /* precomputed scalar u'*B\p */
      79             :   PetscScalar theta;
      80             :   PetscScalar *M;
      81             :   PetscScalar *ps;
      82             :   PetscInt    ld;
      83             :   Vec         *work;
      84             :   Mat         PPr;
      85             :   BV          X;
      86             :   PetscInt    n;
      87             : } PEP_JD_PCSHELL;
      88             : 
      89             : typedef struct {
      90             :   Mat         Pr,Pi;         /* matrix polynomial evaluated at theta */
      91             :   PEP         pep;
      92             :   Vec         *work;
      93             :   PetscScalar theta[2];
      94             : } PEP_JD_MATSHELL;
      95             : 
      96             : /*
      97             :    Duplicate and resize auxiliary basis
      98             : */
      99          55 : static PetscErrorCode PEPJDDuplicateBasis(PEP pep,BV *basis)
     100             : {
     101          55 :   PEP_JD             *pjd = (PEP_JD*)pep->data;
     102          55 :   PetscInt           nloc,m;
     103          55 :   BVType             type;
     104          55 :   BVOrthogType       otype;
     105          55 :   BVOrthogRefineType oref;
     106          55 :   PetscReal          oeta;
     107          55 :   BVOrthogBlockType  oblock;
     108          55 :   VecType            vtype;
     109             : 
     110          55 :   PetscFunctionBegin;
     111          55 :   if (pjd->ld>1) {
     112          55 :     PetscCall(BVCreate(PetscObjectComm((PetscObject)pep),basis));
     113          55 :     PetscCall(BVGetSizes(pep->V,&nloc,NULL,&m));
     114          55 :     nloc += pjd->ld-1;
     115          55 :     PetscCall(BVSetSizes(*basis,nloc,PETSC_DECIDE,m));
     116          55 :     PetscCall(BVGetType(pep->V,&type));
     117          55 :     PetscCall(BVSetType(*basis,type));
     118          55 :     PetscCall(BVGetVecType(pep->V,&vtype));
     119          55 :     PetscCall(BVSetVecType(*basis,vtype));
     120          55 :     PetscCall(BVGetOrthogonalization(pep->V,&otype,&oref,&oeta,&oblock));
     121          55 :     PetscCall(BVSetOrthogonalization(*basis,otype,oref,oeta,oblock));
     122          55 :     PetscCall(PetscObjectStateIncrease((PetscObject)*basis));
     123           0 :   } else PetscCall(BVDuplicate(pep->V,basis));
     124          55 :   PetscFunctionReturn(PETSC_SUCCESS);
     125             : }
     126             : 
     127          11 : static PetscErrorCode PEPSetUp_JD(PEP pep)
     128             : {
     129          11 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
     130          11 :   PetscBool      isprecond,flg;
     131          11 :   PetscRandom    rand;
     132          11 :   PetscInt       i;
     133             : 
     134          11 :   PetscFunctionBegin;
     135          11 :   PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
     136          11 :   if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*pep->n/pep->ncv);
     137          11 :   if (!pep->which) pep->which = PEP_TARGET_MAGNITUDE;
     138          11 :   PetscCheck(pep->which==PEP_TARGET_MAGNITUDE || pep->which==PEP_TARGET_REAL || pep->which==PEP_TARGET_IMAGINARY,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver supports only target which, see PEPSetWhichEigenpairs()");
     139             : 
     140          11 :   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STPRECOND,&isprecond));
     141          11 :   PetscCheck(isprecond,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver only works with PRECOND spectral transformation");
     142             : 
     143          11 :   PetscCall(STGetTransform(pep->st,&flg));
     144          11 :   PetscCheck(!flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The JD solver requires the ST transform flag unset, see STSetTransform()");
     145          11 :   PEPCheckIgnored(pep,PEP_FEATURE_EXTRACT);
     146             : 
     147          11 :   if (!pjd->mmidx) pjd->mmidx = pep->nmat-1;
     148          11 :   pjd->mmidx = PetscMin(pjd->mmidx,pep->nmat-1);
     149          11 :   if (!pjd->keep) pjd->keep = 0.5;
     150          11 :   PetscCall(PEPBasisCoefficients(pep,pep->pbc));
     151          11 :   PetscCall(PEPAllocateSolution(pep,0));
     152          11 :   PetscCall(BVGetRandomContext(pep->V,&rand));  /* make sure the random context is available when duplicating */
     153          11 :   PetscCall(PEPSetWorkVecs(pep,5));
     154          11 :   pjd->ld = pep->nev;
     155             : #if !defined (PETSC_USE_COMPLEX)
     156          11 :   pjd->ld++;
     157             : #endif
     158          11 :   PetscCall(PetscMalloc2(pep->nmat,&pjd->TV,pep->nmat,&pjd->AX));
     159          46 :   for (i=0;i<pep->nmat;i++) PetscCall(PEPJDDuplicateBasis(pep,pjd->TV+i));
     160          11 :   if (pjd->ld>1) {
     161          11 :     PetscCall(PEPJDDuplicateBasis(pep,&pjd->V));
     162          11 :     PetscCall(BVSetFromOptions(pjd->V));
     163          46 :     for (i=0;i<pep->nmat;i++) PetscCall(BVDuplicateResize(pep->V,pjd->ld-1,pjd->AX+i));
     164          11 :     PetscCall(BVDuplicateResize(pep->V,pjd->ld-1,pjd->N));
     165          11 :     PetscCall(BVDuplicateResize(pep->V,pjd->ld-1,pjd->N+1));
     166          11 :     pjd->X = pep->V;
     167          11 :     PetscCall(PetscCalloc3(pjd->ld*pjd->ld,&pjd->XpX,pep->ncv*pep->ncv,&pjd->T,pjd->ld*pjd->ld*pep->nmat,&pjd->Tj));
     168           0 :   } else pjd->V = pep->V;
     169          11 :   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(PEPJDDuplicateBasis(pep,&pjd->W));
     170           2 :   else pjd->W = pjd->V;
     171          11 :   PetscCall(DSSetType(pep->ds,DSPEP));
     172          11 :   PetscCall(DSPEPSetDegree(pep->ds,pep->nmat-1));
     173          11 :   if (pep->basis!=PEP_BASIS_MONOMIAL) PetscCall(DSPEPSetCoefficients(pep->ds,pep->pbc));
     174          11 :   PetscCall(DSAllocate(pep->ds,pep->ncv));
     175          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     176             : }
     177             : 
     178             : /*
     179             :    Updates columns (low to (high-1)) of TV[i]
     180             : */
     181         264 : static PetscErrorCode PEPJDUpdateTV(PEP pep,PetscInt low,PetscInt high,Vec *w)
     182             : {
     183         264 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
     184         264 :   PetscInt       pp,col,i,nloc,nconv;
     185         264 :   Vec            v1,v2,t1,t2;
     186         264 :   PetscScalar    *array1,*array2,*x2,*xx,*N,*Np,*y2=NULL,zero=0.0,sone=1.0,*pT,fact,*psc;
     187         264 :   PetscReal      *cg,*ca,*cb;
     188         264 :   PetscMPIInt    rk,np;
     189         264 :   PetscBLASInt   n_,ld_,one=1;
     190         264 :   Mat            T;
     191         264 :   BV             pbv;
     192             : 
     193         264 :   PetscFunctionBegin;
     194         264 :   ca = pep->pbc; cb = ca+pep->nmat; cg = cb + pep->nmat;
     195         264 :   nconv = pjd->nlock;
     196         264 :   PetscCall(PetscMalloc5(nconv,&x2,nconv,&xx,nconv*nconv,&pT,nconv*nconv,&N,nconv*nconv,&Np));
     197         264 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
     198         264 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
     199         264 :   PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
     200         264 :   t1 = w[0];
     201         264 :   t2 = w[1];
     202         264 :   PetscCall(PetscBLASIntCast(pjd->nlock,&n_));
     203         264 :   PetscCall(PetscBLASIntCast(pjd->ld,&ld_));
     204         264 :   if (nconv) {
     205         189 :     for (i=0;i<nconv;i++) PetscCall(PetscArraycpy(pT+i*nconv,pjd->T+i*pep->ncv,nconv));
     206          68 :     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,pT,&T));
     207             :   }
     208         835 :   for (col=low;col<high;col++) {
     209         571 :     PetscCall(BVGetColumn(pjd->V,col,&v1));
     210         571 :     PetscCall(VecGetArray(v1,&array1));
     211         571 :     if (nconv>0) {
     212         792 :       for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
     213             :     }
     214         571 :     PetscCall(VecPlaceArray(t1,array1));
     215         571 :     if (nconv) {
     216         258 :       PetscCall(BVSetActiveColumns(pjd->N[0],0,nconv));
     217         258 :       PetscCall(BVSetActiveColumns(pjd->N[1],0,nconv));
     218         258 :       PetscCall(BVDotVec(pjd->X,t1,xx));
     219             :     }
     220        2325 :     for (pp=pep->nmat-1;pp>=0;pp--) {
     221        1754 :       PetscCall(BVGetColumn(pjd->TV[pp],col,&v2));
     222        1754 :       PetscCall(VecGetArray(v2,&array2));
     223        1754 :       PetscCall(VecPlaceArray(t2,array2));
     224        1754 :       PetscCall(MatMult(pep->A[pp],t1,t2));
     225        1754 :       if (nconv) {
     226         800 :         if (pp<pep->nmat-3) {
     227          26 :           PetscCall(BVMult(pjd->N[0],1.0,-cg[pp+2],pjd->AX[pp+1],NULL));
     228          26 :           PetscCall(MatShift(T,-cb[pp+1]));
     229          26 :           PetscCall(BVMult(pjd->N[0],1.0/ca[pp],1.0/ca[pp],pjd->N[1],T));
     230          26 :           pbv = pjd->N[0]; pjd->N[0] = pjd->N[1]; pjd->N[1] = pbv;
     231          26 :           PetscCall(BVMultVec(pjd->N[1],1.0,1.0,t2,x2));
     232          26 :           PetscCall(MatShift(T,cb[pp+1]));
     233         774 :         } else if (pp==pep->nmat-3) {
     234         258 :           PetscCall(BVCopy(pjd->AX[pp+2],pjd->N[0]));
     235         258 :           PetscCall(BVScale(pjd->N[0],1/ca[pp+1]));
     236         258 :           PetscCall(BVCopy(pjd->AX[pp+1],pjd->N[1]));
     237         258 :           PetscCall(MatShift(T,-cb[pp+1]));
     238         258 :           PetscCall(BVMult(pjd->N[1],1.0/ca[pp],1.0/ca[pp],pjd->N[0],T));
     239         258 :           PetscCall(BVMultVec(pjd->N[1],1.0,1.0,t2,x2));
     240         258 :           PetscCall(MatShift(T,cb[pp+1]));
     241         516 :         } else if (pp==pep->nmat-2) PetscCall(BVMultVec(pjd->AX[pp+1],1.0/ca[pp],1.0,t2,x2));
     242         800 :         if (pp<pjd->midx) {
     243         397 :           y2 = array2+nloc;
     244         397 :           PetscCallBLAS("BLASgemv",BLASgemv_("C",&n_,&n_,&sone,pjd->Tj+pjd->ld*pjd->ld*pp,&ld_,xx,&one,&zero,y2,&one));
     245         397 :           if (pp<pjd->midx-2) {
     246           0 :             fact = -cg[pp+2];
     247           0 :             PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+(pp+1)*pjd->ld*pjd->ld,&ld_,pjd->XpX,&ld_,&fact,Np,&n_));
     248           0 :             fact = 1/ca[pp];
     249           0 :             PetscCall(MatShift(T,-cb[pp+1]));
     250           0 :             PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&fact,N,&n_,pT,&n_,&fact,Np,&n_));
     251           0 :             PetscCall(MatShift(T,cb[pp+1]));
     252           0 :             psc = Np; Np = N; N = psc;
     253           0 :             PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,N,&n_,x2,&one,&sone,y2,&one));
     254         397 :           } else if (pp==pjd->midx-2) {
     255         139 :             fact = 1/ca[pp];
     256         139 :             PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&fact,pjd->Tj+(pp+1)*pjd->ld*pjd->ld,&ld_,pjd->XpX,&ld_,&zero,N,&n_));
     257         139 :             PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,N,&n_,x2,&one,&sone,y2,&one));
     258         258 :           } else if (pp==pjd->midx-1) PetscCall(PetscArrayzero(Np,nconv*nconv));
     259             :         }
     260        2428 :         for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
     261             :       }
     262        1754 :       PetscCall(VecResetArray(t2));
     263        1754 :       PetscCall(VecRestoreArray(v2,&array2));
     264        1754 :       PetscCall(BVRestoreColumn(pjd->TV[pp],col,&v2));
     265             :     }
     266         571 :     PetscCall(VecResetArray(t1));
     267         571 :     PetscCall(VecRestoreArray(v1,&array1));
     268         571 :     PetscCall(BVRestoreColumn(pjd->V,col,&v1));
     269             :   }
     270         264 :   if (nconv) PetscCall(MatDestroy(&T));
     271         264 :   PetscCall(PetscFree5(x2,xx,pT,N,Np));
     272         264 :   PetscFunctionReturn(PETSC_SUCCESS);
     273             : }
     274             : 
     275             : /*
     276             :    RRQR of X. Xin*P=Xou*R. Rank of R is rk
     277             : */
     278          47 : static PetscErrorCode PEPJDOrthogonalize(PetscInt row,PetscInt col,PetscScalar *X,PetscInt ldx,PetscInt *rk,PetscInt *P,PetscScalar *R,PetscInt ldr)
     279             : {
     280          47 :   PetscInt       i,j,n,r;
     281          47 :   PetscBLASInt   row_,col_,ldx_,*p,lwork,info,n_;
     282          47 :   PetscScalar    *tau,*work;
     283          47 :   PetscReal      tol,*rwork;
     284             : 
     285          47 :   PetscFunctionBegin;
     286          47 :   PetscCall(PetscBLASIntCast(row,&row_));
     287          47 :   PetscCall(PetscBLASIntCast(col,&col_));
     288          47 :   PetscCall(PetscBLASIntCast(ldx,&ldx_));
     289          47 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     290          47 :   n = PetscMin(row,col);
     291          47 :   PetscCall(PetscBLASIntCast(n,&n_));
     292          47 :   lwork = 3*col_+1;
     293          47 :   PetscCall(PetscMalloc4(col,&p,n,&tau,lwork,&work,2*col,&rwork));
     294         442 :   for (i=1;i<col;i++) p[i] = 0;
     295          47 :   p[0] = 1;
     296             : 
     297             :   /* rank revealing QR */
     298             : #if defined(PETSC_USE_COMPLEX)
     299             :   PetscCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,rwork,&info));
     300             : #else
     301          47 :   PetscCallBLAS("LAPACKgeqp3",LAPACKgeqp3_(&row_,&col_,X,&ldx_,p,tau,work,&lwork,&info));
     302             : #endif
     303          47 :   SlepcCheckLapackInfo("geqp3",info);
     304         241 :   if (P) for (i=0;i<col;i++) P[i] = p[i]-1;
     305             : 
     306             :   /* rank computation */
     307          47 :   tol = PetscMax(row,col)*PETSC_MACHINE_EPSILON*PetscAbsScalar(X[0]);
     308          47 :   r = 1;
     309         442 :   for (i=1;i<n;i++) {
     310         395 :     if (PetscAbsScalar(X[i+ldx*i])>tol) r++;
     311             :     else break;
     312             :   }
     313          47 :   if (rk) *rk=r;
     314             : 
     315             :   /* copy upper triangular matrix if requested */
     316          47 :   if (R) {
     317         209 :      for (i=0;i<r;i++) {
     318         194 :        PetscCall(PetscArrayzero(R+i*ldr,r));
     319        1703 :        for (j=0;j<=i;j++) R[i*ldr+j] = X[i*ldx+j];
     320             :      }
     321             :   }
     322          47 :   PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&row_,&n_,&n_,X,&ldx_,tau,work,&lwork,&info));
     323          47 :   SlepcCheckLapackInfo("orgqr",info);
     324          47 :   PetscCall(PetscFPTrapPop());
     325          47 :   PetscCall(PetscFree4(p,tau,work,rwork));
     326          47 :   PetscFunctionReturn(PETSC_SUCCESS);
     327             : }
     328             : 
     329             : /*
     330             :    Application of extended preconditioner
     331             : */
     332        4925 : static PetscErrorCode PEPJDExtendedPCApply(PC pc,Vec x,Vec y)
     333             : {
     334        4925 :   PetscInt          i,j,nloc,n,ld=0;
     335        4925 :   PetscMPIInt       np;
     336        4925 :   Vec               tx,ty;
     337        4925 :   PEP_JD_PCSHELL    *ctx;
     338        4925 :   const PetscScalar *array1;
     339        4925 :   PetscScalar       *x2=NULL,*t=NULL,*ps=NULL,*array2,zero=0.0,sone=1.0;
     340        4925 :   PetscBLASInt      one=1,ld_,n_,ncv_;
     341        4925 :   PEP_JD            *pjd=NULL;
     342             : 
     343        4925 :   PetscFunctionBegin;
     344        4925 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&np));
     345        4925 :   PetscCall(PCShellGetContext(pc,&ctx));
     346        4925 :   n  = ctx->n;
     347        4925 :   if (n) {
     348        2196 :     pjd = (PEP_JD*)ctx->pep->data;
     349        2196 :     ps = ctx->ps;
     350        2196 :     ld = pjd->ld;
     351        2196 :     PetscCall(PetscMalloc2(n,&x2,n,&t));
     352        2196 :     PetscCall(VecGetLocalSize(ctx->work[0],&nloc));
     353        2196 :     PetscCall(VecGetArrayRead(x,&array1));
     354        4982 :     for (i=0;i<n;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
     355        2196 :     PetscCall(VecRestoreArrayRead(x,&array1));
     356             :   }
     357             : 
     358             :   /* y = B\x apply PC */
     359        4925 :   tx = ctx->work[0];
     360        4925 :   ty = ctx->work[1];
     361        4925 :   PetscCall(VecGetArrayRead(x,&array1));
     362        4925 :   PetscCall(VecPlaceArray(tx,array1));
     363        4925 :   PetscCall(VecGetArray(y,&array2));
     364        4925 :   PetscCall(VecPlaceArray(ty,array2));
     365        4925 :   PetscCall(PCApply(ctx->pc,tx,ty));
     366        4925 :   if (n) {
     367        2196 :     PetscCall(PetscBLASIntCast(ld,&ld_));
     368        2196 :     PetscCall(PetscBLASIntCast(n,&n_));
     369        4982 :     for (i=0;i<n;i++) {
     370        2786 :       t[i] = 0.0;
     371        7096 :       for (j=0;j<n;j++) t[i] += ctx->M[i+j*ld]*x2[j];
     372             :     }
     373        2196 :     if (pjd->midx==1) {
     374        1808 :       PetscCall(PetscBLASIntCast(ctx->pep->ncv,&ncv_));
     375        3664 :       for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] -= ctx->theta;
     376        1808 :       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,pjd->T,&ncv_,t,&one,&zero,x2,&one));
     377        3664 :       for (i=0;i<n;i++) pjd->T[i*(1+ctx->pep->ncv)] += ctx->theta;
     378        3664 :       for (i=0;i<n;i++) array2[nloc+i] = x2[i];
     379        3664 :       for (i=0;i<n;i++) x2[i] = -t[i];
     380             :     } else {
     381        1318 :       for (i=0;i<n;i++) array2[nloc+i] = t[i];
     382         388 :       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,ps,&ld_,t,&one,&zero,x2,&one));
     383             :     }
     384        4982 :     for (i=0;i<n;i++) array2[nloc+i] /= PetscSqrtReal(np);
     385        2196 :     PetscCall(BVSetActiveColumns(pjd->X,0,n));
     386        2196 :     PetscCall(BVMultVec(pjd->X,-1.0,1.0,ty,x2));
     387        2196 :     PetscCall(PetscFree2(x2,t));
     388             :   }
     389        4925 :   PetscCall(VecResetArray(tx));
     390        4925 :   PetscCall(VecResetArray(ty));
     391        4925 :   PetscCall(VecRestoreArrayRead(x,&array1));
     392        4925 :   PetscCall(VecRestoreArray(y,&array2));
     393        4925 :   PetscFunctionReturn(PETSC_SUCCESS);
     394             : }
     395             : 
     396             : /*
     397             :    Application of shell preconditioner:
     398             :       y = B\x - eta*B\p,  with eta = (u'*B\x)/(u'*B\p)
     399             : */
     400        4437 : static PetscErrorCode PCShellApply_PEPJD(PC pc,Vec x,Vec y)
     401             : {
     402        4437 :   PetscScalar    rr,eta;
     403        4437 :   PEP_JD_PCSHELL *ctx;
     404        4437 :   PetscInt       sz;
     405        4437 :   const Vec      *xs,*ys;
     406             : #if !defined(PETSC_USE_COMPLEX)
     407        4437 :   PetscScalar    rx,xr,xx;
     408             : #endif
     409             : 
     410        4437 :   PetscFunctionBegin;
     411        4437 :   PetscCall(PCShellGetContext(pc,&ctx));
     412        4437 :   PetscCall(VecCompGetSubVecs(x,&sz,&xs));
     413        4437 :   PetscCall(VecCompGetSubVecs(y,NULL,&ys));
     414             :   /* y = B\x apply extended PC */
     415        4437 :   PetscCall(PEPJDExtendedPCApply(pc,xs[0],ys[0]));
     416             : #if !defined(PETSC_USE_COMPLEX)
     417        4437 :   if (sz==2) PetscCall(PEPJDExtendedPCApply(pc,xs[1],ys[1]));
     418             : #endif
     419             : 
     420             :   /* Compute eta = u'*y / u'*Bp */
     421        4437 :   PetscCall(VecDot(ys[0],ctx->u[0],&rr));
     422        4437 :   eta  = -rr*ctx->gamma[0];
     423             : #if !defined(PETSC_USE_COMPLEX)
     424        4437 :   if (sz==2) {
     425         227 :     PetscCall(VecDot(ys[0],ctx->u[1],&xr));
     426         227 :     PetscCall(VecDot(ys[1],ctx->u[0],&rx));
     427         227 :     PetscCall(VecDot(ys[1],ctx->u[1],&xx));
     428         227 :     eta += -ctx->gamma[0]*xx-ctx->gamma[1]*(-xr+rx);
     429             :   }
     430             : #endif
     431        4437 :   eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
     432             : 
     433             :   /* y = y - eta*Bp */
     434        4437 :   PetscCall(VecAXPY(ys[0],eta,ctx->Bp[0]));
     435             : #if !defined(PETSC_USE_COMPLEX)
     436        4437 :   if (sz==2) {
     437         227 :     PetscCall(VecAXPY(ys[1],eta,ctx->Bp[1]));
     438         227 :     eta = -ctx->gamma[1]*(rr+xx)+ctx->gamma[0]*(-xr+rx);
     439         227 :     eta /= ctx->gamma[0]*ctx->gamma[0]+ctx->gamma[1]*ctx->gamma[1];
     440         227 :     PetscCall(VecAXPY(ys[0],eta,ctx->Bp[1]));
     441         227 :     PetscCall(VecAXPY(ys[1],-eta,ctx->Bp[0]));
     442             :   }
     443             : #endif
     444        4437 :   PetscFunctionReturn(PETSC_SUCCESS);
     445             : }
     446             : 
     447         219 : static PetscErrorCode PEPJDCopyToExtendedVec(PEP pep,Vec v,PetscScalar *a,PetscInt na,PetscInt off,Vec vex,PetscBool back)
     448             : {
     449         219 :   PetscMPIInt    np,rk,count;
     450         219 :   PetscScalar    *array1,*array2;
     451         219 :   PetscInt       nloc;
     452             : 
     453         219 :   PetscFunctionBegin;
     454         219 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
     455         219 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
     456         219 :   PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
     457         219 :   if (v) {
     458          33 :     PetscCall(VecGetArray(v,&array1));
     459          33 :     PetscCall(VecGetArray(vex,&array2));
     460          33 :     if (back) PetscCall(PetscArraycpy(array1,array2,nloc));
     461           0 :     else PetscCall(PetscArraycpy(array2,array1,nloc));
     462          33 :     PetscCall(VecRestoreArray(v,&array1));
     463          33 :     PetscCall(VecRestoreArray(vex,&array2));
     464             :   }
     465         219 :   if (a) {
     466         219 :     PetscCall(VecGetArray(vex,&array2));
     467         219 :     if (back) {
     468          33 :       PetscCall(PetscArraycpy(a,array2+nloc+off,na));
     469          33 :       PetscCall(PetscMPIIntCast(na,&count));
     470          66 :       PetscCallMPI(MPI_Bcast(a,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
     471             :     } else {
     472         186 :       PetscCall(PetscArraycpy(array2+nloc+off,a,na));
     473         186 :       PetscCall(PetscMPIIntCast(na,&count));
     474         372 :       PetscCallMPI(MPI_Bcast(array2+nloc+off,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
     475             :     }
     476         219 :     PetscCall(VecRestoreArray(vex,&array2));
     477             :   }
     478         219 :   PetscFunctionReturn(PETSC_SUCCESS);
     479             : }
     480             : 
     481             : /* Computes Phi^hat(lambda) times a vector or its derivative (depends on beval)
     482             :      if no vector is provided returns a matrix
     483             :  */
     484        9293 : static PetscErrorCode PEPJDEvaluateHatBasis(PEP pep,PetscInt n,PetscScalar *H,PetscInt ldh,PetscScalar *beval,PetscScalar *t,PetscInt idx,PetscScalar *qpp,PetscScalar *qp,PetscScalar *q)
     485             : {
     486        9293 :   PetscInt       j,i;
     487        9293 :   PetscBLASInt   n_,ldh_,one=1;
     488        9293 :   PetscReal      *a,*b,*g;
     489        9293 :   PetscScalar    sone=1.0,zero=0.0;
     490             : 
     491        9293 :   PetscFunctionBegin;
     492        9293 :   a = pep->pbc; b=a+pep->nmat; g=b+pep->nmat;
     493        9293 :   PetscCall(PetscBLASIntCast(n,&n_));
     494        9293 :   PetscCall(PetscBLASIntCast(ldh,&ldh_));
     495        9293 :   if (idx<1) PetscCall(PetscArrayzero(q,t?n:n*n));
     496        6730 :   else if (idx==1) {
     497        6107 :     if (t) {for (j=0;j<n;j++) q[j] = t[j]*beval[idx-1]/a[0];}
     498             :     else {
     499           0 :       PetscCall(PetscArrayzero(q,n*n));
     500           0 :       for (j=0;j<n;j++) q[(j+1)*n] = beval[idx-1]/a[0];
     501             :     }
     502             :   } else {
     503        4167 :     if (t) {
     504        4144 :       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&n_,&sone,H,&ldh_,qp,&one,&zero,q,&one));
     505        9269 :       for (j=0;j<n;j++) {
     506        5125 :         q[j] += beval[idx-1]*t[j]-b[idx-1]*qp[j]-g[idx-1]*qpp[j];
     507        5125 :         q[j] /= a[idx-1];
     508             :       }
     509             :     } else {
     510          23 :       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,H,&ldh_,qp,&n_,&zero,q,&n_));
     511          81 :       for (j=0;j<n;j++) {
     512          58 :         q[j+n*j] += beval[idx-1];
     513         210 :         for (i=0;i<n;i++) {
     514         152 :           q[i+n*j] += -b[idx-1]*qp[j*n+i]-g[idx-1]*qpp[j*n+i];
     515         152 :           q[i+n*j] /= a[idx-1];
     516             :         }
     517             :       }
     518             :     }
     519             :   }
     520        9293 :   PetscFunctionReturn(PETSC_SUCCESS);
     521             : }
     522             : 
     523         705 : static PetscErrorCode PEPJDComputeResidual(PEP pep,PetscBool derivative,PetscInt sz,Vec *u,PetscScalar *theta,Vec *p,Vec *work)
     524             : {
     525         705 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
     526         705 :   PetscMPIInt    rk,np,count;
     527         705 :   Vec            tu,tp,w;
     528         705 :   PetscScalar    *dval,*dvali,*array1,*array2,*x2=NULL,*y2,*qj=NULL,*tt=NULL,*xx=NULL,*xxi=NULL,sone=1.0;
     529         705 :   PetscInt       i,j,nconv,nloc;
     530         705 :   PetscBLASInt   n,ld,one=1;
     531             : #if !defined(PETSC_USE_COMPLEX)
     532         705 :   Vec            tui=NULL,tpi=NULL;
     533         705 :   PetscScalar    *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
     534             : #endif
     535             : 
     536         705 :   PetscFunctionBegin;
     537         705 :   nconv = pjd->nlock;
     538         705 :   if (!nconv) PetscCall(PetscMalloc1(2*sz*pep->nmat,&dval));
     539             :   else {
     540         623 :     PetscCall(PetscMalloc5(2*pep->nmat,&dval,2*nconv,&xx,nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*pep->nmat,&qj));
     541         316 :     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pep),&rk));
     542         316 :     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
     543         316 :     PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
     544         316 :     PetscCall(VecGetArray(u[0],&array1));
     545         942 :     for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]*PetscSqrtReal(np);
     546         316 :     PetscCall(VecRestoreArray(u[0],&array1));
     547             : #if !defined(PETSC_USE_COMPLEX)
     548         316 :     if (sz==2) {
     549           9 :       x2i = x2+nconv;
     550           9 :       PetscCall(VecGetArray(u[1],&arrayi1));
     551          27 :       for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]*PetscSqrtReal(np);
     552           9 :       PetscCall(VecRestoreArray(u[1],&arrayi1));
     553             :     }
     554             : #endif
     555             :   }
     556         705 :   dvali = dval+pep->nmat;
     557         705 :   tu = work[0];
     558         705 :   tp = work[1];
     559         705 :   w  = work[2];
     560         705 :   PetscCall(VecGetArray(u[0],&array1));
     561         705 :   PetscCall(VecPlaceArray(tu,array1));
     562         705 :   PetscCall(VecGetArray(p[0],&array2));
     563         705 :   PetscCall(VecPlaceArray(tp,array2));
     564         705 :   PetscCall(VecSet(tp,0.0));
     565             : #if !defined(PETSC_USE_COMPLEX)
     566         705 :   if (sz==2) {
     567          57 :     tui = work[3];
     568          57 :     tpi = work[4];
     569          57 :     PetscCall(VecGetArray(u[1],&arrayi1));
     570          57 :     PetscCall(VecPlaceArray(tui,arrayi1));
     571          57 :     PetscCall(VecGetArray(p[1],&arrayi2));
     572          57 :     PetscCall(VecPlaceArray(tpi,arrayi2));
     573          57 :     PetscCall(VecSet(tpi,0.0));
     574             :   }
     575             : #endif
     576         705 :   if (derivative) PetscCall(PEPEvaluateBasisDerivative(pep,theta[0],theta[1],dval,dvali));
     577         459 :   else PetscCall(PEPEvaluateBasis(pep,theta[0],theta[1],dval,dvali));
     578        2647 :   for (i=derivative?1:0;i<pep->nmat;i++) {
     579        1942 :     PetscCall(MatMult(pep->A[i],tu,w));
     580        1942 :     PetscCall(VecAXPY(tp,dval[i],w));
     581             : #if !defined(PETSC_USE_COMPLEX)
     582        1942 :     if (sz==2) {
     583         147 :       PetscCall(VecAXPY(tpi,dvali[i],w));
     584         147 :       PetscCall(MatMult(pep->A[i],tui,w));
     585         147 :       PetscCall(VecAXPY(tpi,dval[i],w));
     586        1942 :       PetscCall(VecAXPY(tp,-dvali[i],w));
     587             :     }
     588             : #endif
     589             :   }
     590         705 :   if (nconv) {
     591        1305 :     for (i=0;i<pep->nmat;i++) PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dval,x2,i,i>1?qj+(i-2)*nconv:NULL,i>0?qj+(i-1)*nconv:NULL,qj+i*nconv));
     592             : #if !defined(PETSC_USE_COMPLEX)
     593         316 :     if (sz==2) {
     594           9 :       qji = qj+nconv*pep->nmat;
     595           9 :       qq = qji+nconv*pep->nmat;
     596          36 :       for (i=0;i<pep->nmat;i++) PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dvali,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
     597          63 :       for (i=0;i<nconv*pep->nmat;i++) qj[i] -= qji[i];
     598          36 :       for (i=0;i<pep->nmat;i++) {
     599          27 :         PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dval,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
     600          27 :         PetscCall(PEPJDEvaluateHatBasis(pep,nconv,pjd->T,pep->ncv,dvali,x2,i,i>1?qq+(i-2)*nconv:NULL,i>0?qq+(i-1)*nconv:NULL,qq+i*nconv));
     601             :       }
     602          63 :       for (i=0;i<nconv*pep->nmat;i++) qji[i] += qq[i];
     603          32 :       for (i=derivative?2:1;i<pep->nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,tpi,qji+i*nconv));
     604             :     }
     605             : #endif
     606        1241 :     for (i=derivative?2:1;i<pep->nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,tp,qj+i*nconv));
     607             : 
     608             :     /* extended vector part */
     609         316 :     PetscCall(BVSetActiveColumns(pjd->X,0,nconv));
     610         316 :     PetscCall(BVDotVec(pjd->X,tu,xx));
     611         316 :     xxi = xx+nconv;
     612             : #if !defined(PETSC_USE_COMPLEX)
     613         316 :     if (sz==2) PetscCall(BVDotVec(pjd->X,tui,xxi));
     614             : #endif
     615         316 :     if (sz==1) PetscCall(PetscArrayzero(xxi,nconv));
     616         316 :     if (rk==np-1) {
     617         277 :       PetscCall(PetscBLASIntCast(nconv,&n));
     618         277 :       PetscCall(PetscBLASIntCast(pjd->ld,&ld));
     619         277 :       y2  = array2+nloc;
     620         277 :       PetscCall(PetscArrayzero(y2,nconv));
     621         619 :       for (j=derivative?1:0;j<pjd->midx;j++) {
     622        1081 :         for (i=0;i<nconv;i++) tt[i] = dval[j]*xx[i]-dvali[j]*xxi[i];
     623         342 :         PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
     624         342 :         PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
     625             :       }
     626         805 :       for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
     627             : #if !defined(PETSC_USE_COMPLEX)
     628         277 :       if (sz==2) {
     629           9 :         y2i = arrayi2+nloc;
     630           9 :         PetscCall(PetscArrayzero(y2i,nconv));
     631          23 :         for (j=derivative?1:0;j<pjd->midx;j++) {
     632          42 :           for (i=0;i<nconv;i++) tt[i] = dval[j]*xxi[i]+dvali[j]*xx[i];
     633          14 :           PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
     634          14 :           PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
     635             :         }
     636          27 :         for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
     637             :       }
     638             : #endif
     639             :     }
     640         316 :     PetscCall(PetscMPIIntCast(nconv,&count));
     641         632 :     PetscCallMPI(MPI_Bcast(array2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
     642             : #if !defined(PETSC_USE_COMPLEX)
     643         325 :     if (sz==2) PetscCallMPI(MPI_Bcast(arrayi2+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
     644             : #endif
     645             :   }
     646         705 :   if (nconv) PetscCall(PetscFree5(dval,xx,tt,x2,qj));
     647         389 :   else PetscCall(PetscFree(dval));
     648         705 :   PetscCall(VecResetArray(tu));
     649         705 :   PetscCall(VecRestoreArray(u[0],&array1));
     650         705 :   PetscCall(VecResetArray(tp));
     651         705 :   PetscCall(VecRestoreArray(p[0],&array2));
     652             : #if !defined(PETSC_USE_COMPLEX)
     653         705 :   if (sz==2) {
     654          57 :     PetscCall(VecResetArray(tui));
     655          57 :     PetscCall(VecRestoreArray(u[1],&arrayi1));
     656          57 :     PetscCall(VecResetArray(tpi));
     657          57 :     PetscCall(VecRestoreArray(p[1],&arrayi2));
     658             :   }
     659             : #endif
     660         705 :   PetscFunctionReturn(PETSC_SUCCESS);
     661             : }
     662             : 
     663          11 : static PetscErrorCode PEPJDProcessInitialSpace(PEP pep,Vec *w)
     664             : {
     665          11 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
     666          11 :   PetscScalar    *tt,target[2];
     667          11 :   Vec            vg,wg;
     668          11 :   PetscInt       i;
     669          11 :   PetscReal      norm;
     670             : 
     671          11 :   PetscFunctionBegin;
     672          11 :   PetscCall(PetscMalloc1(pjd->ld-1,&tt));
     673          11 :   PetscCheck(pep->nini==0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Support for initial vectors not implemented yet");
     674          11 :   PetscCall(BVSetRandomColumn(pjd->V,0));
     675          44 :   for (i=0;i<pjd->ld-1;i++) tt[i] = 0.0;
     676          11 :   PetscCall(BVGetColumn(pjd->V,0,&vg));
     677          11 :   PetscCall(PEPJDCopyToExtendedVec(pep,NULL,tt,pjd->ld-1,0,vg,PETSC_FALSE));
     678          11 :   PetscCall(BVRestoreColumn(pjd->V,0,&vg));
     679          11 :   PetscCall(BVNormColumn(pjd->V,0,NORM_2,&norm));
     680          11 :   PetscCall(BVScaleColumn(pjd->V,0,1.0/norm));
     681          11 :   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
     682           9 :     PetscCall(BVGetColumn(pjd->V,0,&vg));
     683           9 :     PetscCall(BVGetColumn(pjd->W,0,&wg));
     684           9 :     PetscCall(VecSet(wg,0.0));
     685           9 :     target[0] = pep->target; target[1] = 0.0;
     686           9 :     PetscCall(PEPJDComputeResidual(pep,PETSC_TRUE,1,&vg,target,&wg,w));
     687           9 :     PetscCall(BVRestoreColumn(pjd->W,0,&wg));
     688           9 :     PetscCall(BVRestoreColumn(pjd->V,0,&vg));
     689           9 :     PetscCall(BVNormColumn(pjd->W,0,NORM_2,&norm));
     690           9 :     PetscCall(BVScaleColumn(pjd->W,0,1.0/norm));
     691             :   }
     692          11 :   PetscCall(PetscFree(tt));
     693          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     694             : }
     695             : 
     696        4200 : static PetscErrorCode MatMult_PEPJD(Mat P,Vec x,Vec y)
     697             : {
     698        4200 :   PEP_JD_MATSHELL *matctx;
     699        4200 :   PEP_JD          *pjd;
     700        4200 :   PetscInt        i,j,nconv,nloc,nmat,ldt,ncv,sz;
     701        4200 :   Vec             tx,ty;
     702        4200 :   const Vec       *xs,*ys;
     703        4200 :   PetscScalar     *array1,*array2,*x2=NULL,*y2,*tt=NULL,*xx=NULL,*xxi,theta[2],sone=1.0,*qj,*val,*vali=NULL;
     704        4200 :   PetscBLASInt    n,ld,one=1;
     705        4200 :   PetscMPIInt     np;
     706             : #if !defined(PETSC_USE_COMPLEX)
     707        4200 :   Vec             txi=NULL,tyi=NULL;
     708        4200 :   PetscScalar     *x2i=NULL,*qji=NULL,*qq,*y2i,*arrayi1,*arrayi2;
     709             : #endif
     710             : 
     711        4200 :   PetscFunctionBegin;
     712        4200 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)P),&np));
     713        4200 :   PetscCall(MatShellGetContext(P,&matctx));
     714        4200 :   pjd   = (PEP_JD*)matctx->pep->data;
     715        4200 :   nconv = pjd->nlock;
     716        4200 :   nmat  = matctx->pep->nmat;
     717        4200 :   ncv   = matctx->pep->ncv;
     718        4200 :   ldt   = pjd->ld;
     719        4200 :   PetscCall(VecCompGetSubVecs(x,&sz,&xs));
     720        4200 :   PetscCall(VecCompGetSubVecs(y,NULL,&ys));
     721        4200 :   theta[0] = matctx->theta[0];
     722        4200 :   theta[1] = (sz==2)?matctx->theta[1]:0.0;
     723        4200 :   if (nconv>0) {
     724        3880 :     PetscCall(PetscMalloc5(nconv,&tt,sz*nconv,&x2,(sz==2?3:1)*nconv*nmat,&qj,2*nconv,&xx,2*nmat,&val));
     725        1980 :     PetscCall(BVGetSizes(matctx->pep->V,&nloc,NULL,NULL));
     726        1980 :     PetscCall(VecGetArray(xs[0],&array1));
     727        4364 :     for (i=0;i<nconv;i++) x2[i] = array1[nloc+i]* PetscSqrtReal(np);
     728        1980 :     PetscCall(VecRestoreArray(xs[0],&array1));
     729             : #if !defined(PETSC_USE_COMPLEX)
     730        1980 :     if (sz==2) {
     731          80 :       x2i = x2+nconv;
     732          80 :       PetscCall(VecGetArray(xs[1],&arrayi1));
     733         240 :       for (i=0;i<nconv;i++) x2i[i] = arrayi1[nloc+i]* PetscSqrtReal(np);
     734          80 :       PetscCall(VecRestoreArray(xs[1],&arrayi1));
     735             :     }
     736             : #endif
     737        1980 :     vali = val+nmat;
     738             :   }
     739        4200 :   tx = matctx->work[0];
     740        4200 :   ty = matctx->work[1];
     741        4200 :   PetscCall(VecGetArray(xs[0],&array1));
     742        4200 :   PetscCall(VecPlaceArray(tx,array1));
     743        4200 :   PetscCall(VecGetArray(ys[0],&array2));
     744        4200 :   PetscCall(VecPlaceArray(ty,array2));
     745        4200 :   PetscCall(MatMult(matctx->Pr,tx,ty));
     746             : #if !defined(PETSC_USE_COMPLEX)
     747        4200 :   if (sz==2) {
     748         203 :     txi = matctx->work[2];
     749         203 :     tyi = matctx->work[3];
     750         203 :     PetscCall(VecGetArray(xs[1],&arrayi1));
     751         203 :     PetscCall(VecPlaceArray(txi,arrayi1));
     752         203 :     PetscCall(VecGetArray(ys[1],&arrayi2));
     753         203 :     PetscCall(VecPlaceArray(tyi,arrayi2));
     754         203 :     PetscCall(MatMult(matctx->Pr,txi,tyi));
     755         203 :     if (theta[1]!=0.0) {
     756         144 :       PetscCall(MatMult(matctx->Pi,txi,matctx->work[4]));
     757         144 :       PetscCall(VecAXPY(ty,-1.0,matctx->work[4]));
     758         144 :       PetscCall(MatMult(matctx->Pi,tx,matctx->work[4]));
     759         144 :       PetscCall(VecAXPY(tyi,1.0,matctx->work[4]));
     760             :     }
     761             :   }
     762             : #endif
     763        4200 :   if (nconv>0) {
     764        1980 :     PetscCall(PEPEvaluateBasis(matctx->pep,theta[0],theta[1],val,vali));
     765        9460 :     for (i=0;i<nmat;i++) PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,ncv,val,x2,i,i>1?qj+(i-2)*nconv:NULL,i>0?qj+(i-1)*nconv:NULL,qj+i*nconv));
     766             : #if !defined(PETSC_USE_COMPLEX)
     767        1980 :     if (sz==2) {
     768          80 :       qji = qj+nconv*nmat;
     769          80 :       qq = qji+nconv*nmat;
     770         320 :       for (i=0;i<nmat;i++) PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,matctx->pep->ncv,vali,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
     771         560 :       for (i=0;i<nconv*nmat;i++) qj[i] -= qji[i];
     772         320 :       for (i=0;i<nmat;i++) {
     773         240 :         PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,matctx->pep->ncv,val,x2i,i,i>1?qji+(i-2)*nconv:NULL,i>0?qji+(i-1)*nconv:NULL,qji+i*nconv));
     774         240 :         PetscCall(PEPJDEvaluateHatBasis(matctx->pep,nconv,pjd->T,matctx->pep->ncv,vali,x2,i,i>1?qq+(i-2)*nconv:NULL,i>0?qq+(i-1)*nconv:NULL,qq+i*nconv));
     775             :       }
     776         560 :       for (i=0;i<nconv*nmat;i++) qji[i] += qq[i];
     777         240 :       for (i=1;i<matctx->pep->nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,tyi,qji+i*nconv));
     778             :     }
     779             : #endif
     780        7480 :     for (i=1;i<nmat;i++) PetscCall(BVMultVec(pjd->AX[i],1.0,1.0,ty,qj+i*nconv));
     781             : 
     782             :     /* extended vector part */
     783        1980 :     PetscCall(BVSetActiveColumns(pjd->X,0,nconv));
     784        1980 :     PetscCall(BVDotVec(pjd->X,tx,xx));
     785        1980 :     xxi = xx+nconv;
     786             : #if !defined(PETSC_USE_COMPLEX)
     787        1980 :     if (sz==2) PetscCall(BVDotVec(pjd->X,txi,xxi));
     788             : #endif
     789        1980 :     if (sz==1) PetscCall(PetscArrayzero(xxi,nconv));
     790        1980 :     PetscCall(PetscBLASIntCast(pjd->nlock,&n));
     791        1980 :     PetscCall(PetscBLASIntCast(ldt,&ld));
     792        1980 :     y2 = array2+nloc;
     793        1980 :     PetscCall(PetscArrayzero(y2,nconv));
     794        4208 :     for (j=0;j<pjd->midx;j++) {
     795        5232 :       for (i=0;i<nconv;i++) tt[i] = val[j]*xx[i]-vali[j]*xxi[i];
     796        2228 :       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qj+j*nconv,&one,&sone,tt,&one));
     797        2228 :       PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2,&one));
     798             :     }
     799             : #if !defined(PETSC_USE_COMPLEX)
     800        1980 :     if (sz==2) {
     801          80 :       y2i = arrayi2+nloc;
     802          80 :       PetscCall(PetscArrayzero(y2i,nconv));
     803         240 :       for (j=0;j<pjd->midx;j++) {
     804         480 :         for (i=0;i<nconv;i++) tt[i] = val[j]*xxi[i]+vali[j]*xx[i];
     805         160 :         PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,pjd->XpX,&ld,qji+j*nconv,&one,&sone,tt,&one));
     806         160 :         PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&sone,pjd->Tj+j*ld*ld,&ld,tt,&one,&sone,y2i,&one));
     807             :       }
     808         240 :       for (i=0;i<nconv;i++) arrayi2[nloc+i] /= PetscSqrtReal(np);
     809             :     }
     810             : #endif
     811        4364 :     for (i=0;i<nconv;i++) array2[nloc+i] /= PetscSqrtReal(np);
     812        1980 :     PetscCall(PetscFree5(tt,x2,qj,xx,val));
     813             :   }
     814        4200 :   PetscCall(VecResetArray(tx));
     815        4200 :   PetscCall(VecRestoreArray(xs[0],&array1));
     816        4200 :   PetscCall(VecResetArray(ty));
     817        4200 :   PetscCall(VecRestoreArray(ys[0],&array2));
     818             : #if !defined(PETSC_USE_COMPLEX)
     819        4200 :   if (sz==2) {
     820         203 :     PetscCall(VecResetArray(txi));
     821         203 :     PetscCall(VecRestoreArray(xs[1],&arrayi1));
     822         203 :     PetscCall(VecResetArray(tyi));
     823         203 :     PetscCall(VecRestoreArray(ys[1],&arrayi2));
     824             :   }
     825             : #endif
     826        4200 :   PetscFunctionReturn(PETSC_SUCCESS);
     827             : }
     828             : 
     829          11 : static PetscErrorCode MatCreateVecs_PEPJD(Mat A,Vec *right,Vec *left)
     830             : {
     831          11 :   PEP_JD_MATSHELL *matctx;
     832          11 :   PEP_JD          *pjd;
     833          11 :   PetscInt        kspsf=1,i;
     834          11 :   Vec             v[2];
     835             : 
     836          11 :   PetscFunctionBegin;
     837          11 :   PetscCall(MatShellGetContext(A,&matctx));
     838          11 :   pjd   = (PEP_JD*)matctx->pep->data;
     839             : #if !defined (PETSC_USE_COMPLEX)
     840          11 :   kspsf = 2;
     841             : #endif
     842          33 :   for (i=0;i<kspsf;i++) PetscCall(BVCreateVec(pjd->V,v+i));
     843          11 :   if (right) PetscCall(VecCreateCompWithVecs(v,kspsf,pjd->vtempl,right));
     844          11 :   if (left) PetscCall(VecCreateCompWithVecs(v,kspsf,pjd->vtempl,left));
     845          33 :   for (i=0;i<kspsf;i++) PetscCall(VecDestroy(&v[i]));
     846          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     847             : }
     848             : 
     849         119 : static PetscErrorCode PEPJDUpdateExtendedPC(PEP pep,PetscScalar theta)
     850             : {
     851         119 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
     852         119 :   PEP_JD_PCSHELL *pcctx;
     853         119 :   PetscInt       i,j,k,n=pjd->nlock,ld=pjd->ld,deg=pep->nmat-1;
     854         119 :   PetscScalar    *M,*ps,*work,*U,*V,*S,*Sp,*Spp,snone=-1.0,sone=1.0,zero=0.0,*val;
     855         119 :   PetscReal      tol,maxeig=0.0,*sg,*rwork;
     856         119 :   PetscBLASInt   n_,info,ld_,*p,lw_,rk=0;
     857             : 
     858         119 :   PetscFunctionBegin;
     859         119 :   if (n) {
     860          48 :     PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
     861          48 :     pcctx->theta = theta;
     862          48 :     pcctx->n = n;
     863          48 :     M  = pcctx->M;
     864          48 :     PetscCall(PetscBLASIntCast(n,&n_));
     865          48 :     PetscCall(PetscBLASIntCast(ld,&ld_));
     866          48 :     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     867          48 :     if (pjd->midx==1) {
     868          25 :       PetscCall(PetscArraycpy(M,pjd->XpX,ld*ld));
     869          25 :       PetscCall(PetscCalloc2(10*n,&work,n,&p));
     870             :     } else {
     871          23 :       ps = pcctx->ps;
     872          23 :       PetscCall(PetscCalloc7(2*n*n,&U,3*n*n,&S,n,&sg,10*n,&work,5*n,&rwork,n,&p,deg+1,&val));
     873          23 :       V = U+n*n;
     874             :       /* pseudo-inverse */
     875          81 :       for (j=0;j<n;j++) {
     876         210 :         for (i=0;i<n;i++) S[n*j+i] = -pjd->T[pep->ncv*j+i];
     877          58 :         S[n*j+j] += theta;
     878             :       }
     879          23 :       lw_ = 10*n_;
     880             : #if !defined (PETSC_USE_COMPLEX)
     881          23 :       PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,&info));
     882             : #else
     883             :       PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&n_,&n_,S,&n_,sg,U,&n_,V,&n_,work,&lw_,rwork,&info));
     884             : #endif
     885          23 :       SlepcCheckLapackInfo("gesvd",info);
     886          81 :       for (i=0;i<n;i++) maxeig = PetscMax(maxeig,sg[i]);
     887          23 :       tol = 10*PETSC_MACHINE_EPSILON*n*maxeig;
     888          81 :       for (j=0;j<n;j++) {
     889          58 :         if (sg[j]>tol) {
     890         210 :           for (i=0;i<n;i++) U[j*n+i] /= sg[j];
     891          58 :           rk++;
     892             :         } else break;
     893             :       }
     894          23 :       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&rk,&sone,U,&n_,V,&n_,&zero,ps,&ld_));
     895             : 
     896             :       /* compute M */
     897          23 :       PetscCall(PEPEvaluateBasis(pep,theta,0.0,val,NULL));
     898          23 :       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&snone,pjd->XpX,&ld_,ps,&ld_,&zero,M,&ld_));
     899          23 :       PetscCall(PetscArrayzero(S,2*n*n));
     900          23 :       Sp = S+n*n;
     901          81 :       for (j=0;j<n;j++) S[j*(n+1)] = 1.0;
     902          46 :       for (k=1;k<pjd->midx;k++) {
     903         233 :         for (j=0;j<n;j++) for (i=0;i<n;i++) V[j*n+i] = S[j*n+i] - ps[j*ld+i]*val[k];
     904          23 :         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pjd->XpX,&ld_,V,&n_,&zero,U,&n_));
     905          23 :         PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,pjd->Tj+k*ld*ld,&ld_,U,&n_,&sone,M,&ld_));
     906          23 :         Spp = Sp; Sp = S;
     907          23 :         PetscCall(PEPJDEvaluateHatBasis(pep,n,pjd->T,pep->ncv,val,NULL,k+1,Spp,Sp,S));
     908             :       }
     909             :     }
     910             :     /* inverse */
     911          48 :     PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,M,&ld_,p,&info));
     912          48 :     SlepcCheckLapackInfo("getrf",info);
     913          48 :     PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&n_,M,&ld_,p,work,&n_,&info));
     914          48 :     SlepcCheckLapackInfo("getri",info);
     915          48 :     PetscCall(PetscFPTrapPop());
     916          48 :     if (pjd->midx==1) PetscCall(PetscFree2(work,p));
     917          23 :     else PetscCall(PetscFree7(U,S,sg,work,rwork,p,val));
     918             :   }
     919         119 :   PetscFunctionReturn(PETSC_SUCCESS);
     920             : }
     921             : 
     922         248 : static PetscErrorCode PEPJDMatSetUp(PEP pep,PetscInt sz,PetscScalar *theta)
     923             : {
     924         248 :   PEP_JD          *pjd = (PEP_JD*)pep->data;
     925         248 :   PEP_JD_MATSHELL *matctx;
     926         248 :   PEP_JD_PCSHELL  *pcctx;
     927         248 :   MatStructure    str;
     928         248 :   PetscScalar     *vals,*valsi;
     929         248 :   PetscBool       skipmat=PETSC_FALSE;
     930         248 :   PetscInt        i;
     931         248 :   Mat             Pr=NULL;
     932             : 
     933         248 :   PetscFunctionBegin;
     934         248 :   if (sz==2 && theta[1]==0.0) sz = 1;
     935         248 :   PetscCall(MatShellGetContext(pjd->Pshell,&matctx));
     936         248 :   PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
     937         248 :   if (matctx->Pr && matctx->theta[0]==theta[0] && ((!matctx->Pi && sz==1) || (sz==2 && matctx->theta[1]==theta[1]))) {
     938         127 :     if (pcctx->n == pjd->nlock) PetscFunctionReturn(PETSC_SUCCESS);
     939             :     skipmat = PETSC_TRUE;
     940             :   }
     941             :   if (!skipmat) {
     942         121 :     PetscCall(PetscMalloc2(pep->nmat,&vals,pep->nmat,&valsi));
     943         121 :     PetscCall(STGetMatStructure(pep->st,&str));
     944         121 :     PetscCall(PEPEvaluateBasis(pep,theta[0],theta[1],vals,valsi));
     945         121 :     if (!matctx->Pr) PetscCall(MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pr));
     946         110 :     else PetscCall(MatCopy(pep->A[0],matctx->Pr,str));
     947         374 :     for (i=1;i<pep->nmat;i++) PetscCall(MatAXPY(matctx->Pr,vals[i],pep->A[i],str));
     948         121 :     if (!pjd->reusepc) {
     949         118 :       if (pcctx->PPr && sz==2) {
     950           6 :         PetscCall(MatCopy(matctx->Pr,pcctx->PPr,str));
     951           6 :         Pr = pcctx->PPr;
     952         112 :       } else Pr = matctx->Pr;
     953             :     }
     954         121 :     matctx->theta[0] = theta[0];
     955             : #if !defined(PETSC_USE_COMPLEX)
     956         121 :     if (sz==2) {
     957           8 :       if (!matctx->Pi) PetscCall(MatDuplicate(pep->A[0],MAT_COPY_VALUES,&matctx->Pi));
     958           5 :       else PetscCall(MatCopy(pep->A[1],matctx->Pi,str));
     959           8 :       PetscCall(MatScale(matctx->Pi,valsi[1]));
     960          16 :       for (i=2;i<pep->nmat;i++) PetscCall(MatAXPY(matctx->Pi,valsi[i],pep->A[i],str));
     961           8 :       matctx->theta[1] = theta[1];
     962             :     }
     963             : #endif
     964         121 :     PetscCall(PetscFree2(vals,valsi));
     965             :   }
     966         122 :   if (!pjd->reusepc) {
     967         119 :     if (!skipmat) {
     968         118 :       PetscCall(PCSetOperators(pcctx->pc,Pr,Pr));
     969         118 :       PetscCall(PCSetUp(pcctx->pc));
     970             :     }
     971         119 :     PetscCall(PEPJDUpdateExtendedPC(pep,theta[0]));
     972             :   }
     973         122 :   PetscFunctionReturn(PETSC_SUCCESS);
     974             : }
     975             : 
     976          11 : static PetscErrorCode PEPJDCreateShellPC(PEP pep,Vec *ww)
     977             : {
     978          11 :   PEP_JD          *pjd = (PEP_JD*)pep->data;
     979          11 :   PEP_JD_PCSHELL  *pcctx;
     980          11 :   PEP_JD_MATSHELL *matctx;
     981          11 :   KSP             ksp;
     982          11 :   PetscInt        nloc,mloc,kspsf=1;
     983          11 :   Vec             v[2];
     984          11 :   PetscScalar     target[2];
     985          11 :   Mat             Pr;
     986             : 
     987          11 :   PetscFunctionBegin;
     988             :   /* Create the reference vector */
     989          11 :   PetscCall(BVGetColumn(pjd->V,0,&v[0]));
     990          11 :   v[1] = v[0];
     991             : #if !defined (PETSC_USE_COMPLEX)
     992          11 :   kspsf = 2;
     993             : #endif
     994          11 :   PetscCall(VecCreateCompWithVecs(v,kspsf,NULL,&pjd->vtempl));
     995          11 :   PetscCall(BVRestoreColumn(pjd->V,0,&v[0]));
     996             : 
     997             :   /* Replace preconditioner with one containing projectors */
     998          11 :   PetscCall(PCCreate(PetscObjectComm((PetscObject)pep),&pjd->pcshell));
     999          11 :   PetscCall(PCSetType(pjd->pcshell,PCSHELL));
    1000          11 :   PetscCall(PCShellSetName(pjd->pcshell,"PCPEPJD"));
    1001          11 :   PetscCall(PCShellSetApply(pjd->pcshell,PCShellApply_PEPJD));
    1002          11 :   PetscCall(PetscNew(&pcctx));
    1003          11 :   PetscCall(PCShellSetContext(pjd->pcshell,pcctx));
    1004          11 :   PetscCall(STGetKSP(pep->st,&ksp));
    1005          11 :   PetscCall(BVCreateVec(pjd->V,&pcctx->Bp[0]));
    1006          11 :   PetscCall(VecDuplicate(pcctx->Bp[0],&pcctx->Bp[1]));
    1007          11 :   PetscCall(KSPGetPC(ksp,&pcctx->pc));
    1008          11 :   PetscCall(PetscObjectReference((PetscObject)pcctx->pc));
    1009          11 :   PetscCall(MatGetLocalSize(pep->A[0],&mloc,&nloc));
    1010          11 :   if (pjd->ld>1) {
    1011          11 :     nloc += pjd->ld-1; mloc += pjd->ld-1;
    1012             :   }
    1013          11 :   PetscCall(PetscNew(&matctx));
    1014          11 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),kspsf*nloc,kspsf*mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,&pjd->Pshell));
    1015          11 :   PetscCall(MatShellSetOperation(pjd->Pshell,MATOP_MULT,(void(*)(void))MatMult_PEPJD));
    1016          11 :   PetscCall(MatShellSetOperation(pjd->Pshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_PEPJD));
    1017          11 :   matctx->pep = pep;
    1018          11 :   target[0] = pep->target; target[1] = 0.0;
    1019          11 :   PetscCall(PEPJDMatSetUp(pep,1,target));
    1020          11 :   Pr = matctx->Pr;
    1021          11 :   pcctx->PPr = NULL;
    1022             : #if !defined(PETSC_USE_COMPLEX)
    1023          11 :   if (!pjd->reusepc) {
    1024          10 :     PetscCall(MatDuplicate(matctx->Pr,MAT_COPY_VALUES,&pcctx->PPr));
    1025          10 :     Pr = pcctx->PPr;
    1026             :   }
    1027             : #endif
    1028          11 :   PetscCall(PCSetOperators(pcctx->pc,Pr,Pr));
    1029          11 :   PetscCall(PCSetErrorIfFailure(pcctx->pc,PETSC_TRUE));
    1030          11 :   PetscCall(KSPSetPC(ksp,pjd->pcshell));
    1031          11 :   if (pjd->reusepc) {
    1032           1 :     PetscCall(PCSetReusePreconditioner(pcctx->pc,PETSC_TRUE));
    1033           1 :     PetscCall(KSPSetReusePreconditioner(ksp,PETSC_TRUE));
    1034             :   }
    1035          11 :   PetscCall(PEP_KSPSetOperators(ksp,pjd->Pshell,pjd->Pshell));
    1036          11 :   PetscCall(KSPSetUp(ksp));
    1037          11 :   if (pjd->ld>1) {
    1038          11 :     PetscCall(PetscMalloc2(pjd->ld*pjd->ld,&pcctx->M,pjd->ld*pjd->ld,&pcctx->ps));
    1039          11 :     pcctx->pep = pep;
    1040             :   }
    1041          11 :   matctx->work = ww;
    1042          11 :   pcctx->work  = ww;
    1043          11 :   PetscFunctionReturn(PETSC_SUCCESS);
    1044             : }
    1045             : 
    1046          11 : static PetscErrorCode PEPJDEigenvectors(PEP pep)
    1047             : {
    1048          11 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
    1049          11 :   PetscBLASInt   ld,nconv,info,nc;
    1050          11 :   PetscScalar    *Z;
    1051          11 :   PetscReal      *wr;
    1052          11 :   Mat            U;
    1053             : #if defined(PETSC_USE_COMPLEX)
    1054             :   PetscScalar    *w;
    1055             : #endif
    1056             : 
    1057          11 :   PetscFunctionBegin;
    1058          11 :   PetscCall(PetscBLASIntCast(pep->ncv,&ld));
    1059          11 :   PetscCall(PetscBLASIntCast(pep->nconv,&nconv));
    1060          11 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1061             : #if !defined(PETSC_USE_COMPLEX)
    1062          11 :   PetscCall(PetscMalloc2(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr));
    1063          11 :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,wr,&info));
    1064             : #else
    1065             :   PetscCall(PetscMalloc3(pep->nconv*pep->nconv,&Z,3*pep->ncv,&wr,2*pep->ncv,&w));
    1066             :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","A",NULL,&nconv,pjd->T,&ld,NULL,&nconv,Z,&nconv,&nconv,&nc,w,wr,&info));
    1067             : #endif
    1068          11 :   PetscCall(PetscFPTrapPop());
    1069          11 :   SlepcCheckLapackInfo("trevc",info);
    1070          11 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nconv,nconv,Z,&U));
    1071          11 :   PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv));
    1072          11 :   PetscCall(BVMultInPlace(pjd->X,U,0,pep->nconv));
    1073          11 :   PetscCall(BVNormalize(pjd->X,pep->eigi));
    1074          11 :   PetscCall(MatDestroy(&U));
    1075             : #if !defined(PETSC_USE_COMPLEX)
    1076          11 :   PetscCall(PetscFree2(Z,wr));
    1077             : #else
    1078             :   PetscCall(PetscFree3(Z,wr,w));
    1079             : #endif
    1080          11 :   PetscFunctionReturn(PETSC_SUCCESS);
    1081             : }
    1082             : 
    1083          15 : static PetscErrorCode PEPJDLockConverged(PEP pep,PetscInt *nv,PetscInt sz)
    1084             : {
    1085          15 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
    1086          15 :   PetscInt       j,i,*P,ldds,rk=0,nvv=*nv;
    1087          15 :   Vec            v,x,w;
    1088          15 :   PetscScalar    *R,*r,*pX,target[2];
    1089          15 :   Mat            X;
    1090          15 :   PetscBLASInt   sz_,rk_,nv_,info;
    1091          15 :   PetscMPIInt    np;
    1092             : 
    1093          15 :   PetscFunctionBegin;
    1094             :   /* update AX and XpX */
    1095          34 :   for (i=sz;i>0;i--) {
    1096          19 :     PetscCall(BVGetColumn(pjd->X,pjd->nlock-i,&x));
    1097          78 :     for (j=0;j<pep->nmat;j++) {
    1098          59 :       PetscCall(BVGetColumn(pjd->AX[j],pjd->nlock-i,&v));
    1099          59 :       PetscCall(MatMult(pep->A[j],x,v));
    1100          59 :       PetscCall(BVRestoreColumn(pjd->AX[j],pjd->nlock-i,&v));
    1101          59 :       PetscCall(BVSetActiveColumns(pjd->AX[j],0,pjd->nlock-i+1));
    1102             :     }
    1103          19 :     PetscCall(BVRestoreColumn(pjd->X,pjd->nlock-i,&x));
    1104          19 :     PetscCall(BVDotColumn(pjd->X,(pjd->nlock-i),pjd->XpX+(pjd->nlock-i)*pjd->ld));
    1105          19 :     pjd->XpX[(pjd->nlock-i)*(1+pjd->ld)] = 1.0;
    1106          35 :     for (j=0;j<pjd->nlock-i;j++) pjd->XpX[j*pjd->ld+pjd->nlock-i] = PetscConj(pjd->XpX[(pjd->nlock-i)*pjd->ld+j]);
    1107             :   }
    1108             : 
    1109             :   /* minimality index */
    1110          15 :   pjd->midx = PetscMin(pjd->mmidx,pjd->nlock);
    1111             : 
    1112             :   /* evaluate the polynomial basis in T */
    1113          15 :   PetscCall(PetscArrayzero(pjd->Tj,pjd->ld*pjd->ld*pep->nmat));
    1114          62 :   for (j=0;j<pep->nmat;j++) PetscCall(PEPEvaluateBasisMat(pep,pjd->nlock,pjd->T,pep->ncv,j,(j>1)?pjd->Tj+(j-2)*pjd->ld*pjd->ld:NULL,pjd->ld,j?pjd->Tj+(j-1)*pjd->ld*pjd->ld:NULL,pjd->ld,pjd->Tj+j*pjd->ld*pjd->ld,pjd->ld));
    1115             : 
    1116             :   /* Extend search space */
    1117          15 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
    1118          15 :   PetscCall(PetscCalloc3(nvv,&P,nvv*nvv,&R,nvv*sz,&r));
    1119          15 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
    1120          15 :   PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
    1121          15 :   PetscCall(PEPJDOrthogonalize(nvv,nvv,pX,ldds,&rk,P,R,nvv));
    1122          34 :   for (j=0;j<sz;j++) {
    1123         262 :     for (i=0;i<rk;i++) r[i*sz+j] = PetscConj(R[nvv*i+j]*pep->eigr[P[i]]); /* first row scaled with permuted diagonal */
    1124             :   }
    1125          15 :   PetscCall(PetscBLASIntCast(rk,&rk_));
    1126          15 :   PetscCall(PetscBLASIntCast(sz,&sz_));
    1127          15 :   PetscCall(PetscBLASIntCast(nvv,&nv_));
    1128          15 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1129          15 :   PetscCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&rk_,R,&nv_,&info));
    1130          15 :   PetscCall(PetscFPTrapPop());
    1131          15 :   SlepcCheckLapackInfo("trtri",info);
    1132          34 :   for (i=0;i<sz;i++) PetscCallBLAS("BLAStrmv",BLAStrmv_("U","C","N",&rk_,R,&nv_,r+i,&sz_));
    1133         258 :   for (i=0;i<sz*rk;i++) r[i] = PetscConj(r[i])/PetscSqrtReal(np); /* revert */
    1134          15 :   PetscCall(BVSetActiveColumns(pjd->V,0,nvv));
    1135          15 :   rk -= sz;
    1136         190 :   for (j=0;j<rk;j++) PetscCall(PetscArraycpy(R+j*nvv,pX+(j+sz)*ldds,nvv));
    1137          15 :   PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
    1138          15 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nvv,rk,R,&X));
    1139          15 :   PetscCall(BVMultInPlace(pjd->V,X,0,rk));
    1140          15 :   PetscCall(MatDestroy(&X));
    1141          15 :   PetscCall(BVSetActiveColumns(pjd->V,0,rk));
    1142         190 :   for (j=0;j<rk;j++) {
    1143         175 :     PetscCall(BVGetColumn(pjd->V,j,&v));
    1144         175 :     PetscCall(PEPJDCopyToExtendedVec(pep,NULL,r+sz*(j+sz),sz,pjd->nlock-sz,v,PETSC_FALSE));
    1145         175 :     PetscCall(BVRestoreColumn(pjd->V,j,&v));
    1146             :   }
    1147          15 :   PetscCall(BVOrthogonalize(pjd->V,NULL));
    1148             : 
    1149          15 :   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
    1150         181 :     for (j=0;j<rk;j++) {
    1151             :       /* W = P(target)*V */
    1152         168 :       PetscCall(BVGetColumn(pjd->W,j,&w));
    1153         168 :       PetscCall(BVGetColumn(pjd->V,j,&v));
    1154         168 :       target[0] = pep->target; target[1] = 0.0;
    1155         168 :       PetscCall(PEPJDComputeResidual(pep,PETSC_FALSE,1,&v,target,&w,pep->work));
    1156         168 :       PetscCall(BVRestoreColumn(pjd->V,j,&v));
    1157         168 :       PetscCall(BVRestoreColumn(pjd->W,j,&w));
    1158             :     }
    1159          13 :     PetscCall(BVSetActiveColumns(pjd->W,0,rk));
    1160          13 :     PetscCall(BVOrthogonalize(pjd->W,NULL));
    1161             :   }
    1162          15 :   *nv = rk;
    1163          15 :   PetscCall(PetscFree3(P,R,r));
    1164          15 :   PetscFunctionReturn(PETSC_SUCCESS);
    1165             : }
    1166             : 
    1167         237 : static PetscErrorCode PEPJDSystemSetUp(PEP pep,PetscInt sz,PetscScalar *theta,Vec *u,Vec *p,Vec *ww)
    1168             : {
    1169         237 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
    1170         237 :   PEP_JD_PCSHELL *pcctx;
    1171             : #if !defined(PETSC_USE_COMPLEX)
    1172         237 :   PetscScalar    s[2];
    1173             : #endif
    1174             : 
    1175         237 :   PetscFunctionBegin;
    1176         237 :   PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
    1177         237 :   PetscCall(PEPJDMatSetUp(pep,sz,theta));
    1178         237 :   pcctx->u[0] = u[0]; pcctx->u[1] = u[1];
    1179             :   /* Compute r'. p is a work space vector */
    1180         237 :   PetscCall(PEPJDComputeResidual(pep,PETSC_TRUE,sz,u,theta,p,ww));
    1181         237 :   PetscCall(PEPJDExtendedPCApply(pjd->pcshell,p[0],pcctx->Bp[0]));
    1182         237 :   PetscCall(VecDot(pcctx->Bp[0],u[0],pcctx->gamma));
    1183             : #if !defined(PETSC_USE_COMPLEX)
    1184         237 :   if (sz==2) {
    1185          24 :     PetscCall(PEPJDExtendedPCApply(pjd->pcshell,p[1],pcctx->Bp[1]));
    1186          24 :     PetscCall(VecDot(pcctx->Bp[0],u[1],pcctx->gamma+1));
    1187          24 :     PetscCall(VecMDot(pcctx->Bp[1],2,u,s));
    1188          24 :     pcctx->gamma[0] += s[1];
    1189          24 :     pcctx->gamma[1] = -pcctx->gamma[1]+s[0];
    1190             :   }
    1191             : #endif
    1192         237 :   if (sz==1) {
    1193         213 :     PetscCall(VecZeroEntries(pcctx->Bp[1]));
    1194         213 :     pcctx->gamma[1] = 0.0;
    1195             :   }
    1196         237 :   PetscFunctionReturn(PETSC_SUCCESS);
    1197             : }
    1198             : 
    1199          11 : static PetscErrorCode PEPSolve_JD(PEP pep)
    1200             : {
    1201          11 :   PEP_JD          *pjd = (PEP_JD*)pep->data;
    1202          11 :   PetscInt        k,nv,nvc,ld,minv,dim,bupdated=0,sz=1,kspsf=1,idx,off,maxits,nloc;
    1203          11 :   PetscMPIInt     np,count;
    1204          11 :   PetscScalar     theta[2]={0.0,0.0},ritz[2]={0.0,0.0},*pX,*eig,*eigi,*array;
    1205          11 :   PetscReal       norm,*res,tol=0.0,rtol,abstol, dtol;
    1206          11 :   PetscBool       lindep,ini=PETSC_TRUE;
    1207          11 :   Vec             tc,t[2]={NULL,NULL},u[2]={NULL,NULL},p[2]={NULL,NULL};
    1208          11 :   Vec             rc,rr[2],r[2]={NULL,NULL},*ww=pep->work,v[2];
    1209          11 :   Mat             G,X,Y;
    1210          11 :   KSP             ksp;
    1211          11 :   PEP_JD_PCSHELL  *pcctx;
    1212          11 :   PEP_JD_MATSHELL *matctx;
    1213             : #if !defined(PETSC_USE_COMPLEX)
    1214          11 :   PetscReal       norm1;
    1215             : #endif
    1216             : 
    1217          11 :   PetscFunctionBegin;
    1218          11 :   PetscCall(PetscCitationsRegister(citation,&cited));
    1219          11 :   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
    1220          11 :   PetscCall(BVGetSizes(pep->V,&nloc,NULL,NULL));
    1221          11 :   PetscCall(DSGetLeadingDimension(pep->ds,&ld));
    1222          11 :   PetscCall(PetscCalloc3(pep->ncv+pep->nev,&eig,pep->ncv+pep->nev,&eigi,pep->ncv+pep->nev,&res));
    1223          11 :   pjd->nlock = 0;
    1224          11 :   PetscCall(STGetKSP(pep->st,&ksp));
    1225          11 :   PetscCall(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
    1226             : #if !defined (PETSC_USE_COMPLEX)
    1227          11 :   kspsf = 2;
    1228             : #endif
    1229          11 :   PetscCall(PEPJDProcessInitialSpace(pep,ww));
    1230          11 :   nv = (pep->nini)?pep->nini:1;
    1231             : 
    1232             :   /* Replace preconditioner with one containing projectors */
    1233          11 :   PetscCall(PEPJDCreateShellPC(pep,ww));
    1234          11 :   PetscCall(PCShellGetContext(pjd->pcshell,&pcctx));
    1235             : 
    1236             :   /* Create auxiliary vectors */
    1237          11 :   PetscCall(BVCreateVec(pjd->V,&u[0]));
    1238          11 :   PetscCall(VecDuplicate(u[0],&p[0]));
    1239          11 :   PetscCall(VecDuplicate(u[0],&r[0]));
    1240             : #if !defined (PETSC_USE_COMPLEX)
    1241          11 :   PetscCall(VecDuplicate(u[0],&u[1]));
    1242          11 :   PetscCall(VecDuplicate(u[0],&p[1]));
    1243          11 :   PetscCall(VecDuplicate(u[0],&r[1]));
    1244             : #endif
    1245             : 
    1246             :   /* Restart loop */
    1247         275 :   while (pep->reason == PEP_CONVERGED_ITERATING) {
    1248         264 :     pep->its++;
    1249         264 :     PetscCall(DSSetDimensions(pep->ds,nv,0,0));
    1250         264 :     PetscCall(BVSetActiveColumns(pjd->V,bupdated,nv));
    1251         264 :     PetscCall(PEPJDUpdateTV(pep,bupdated,nv,ww));
    1252         264 :     if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(BVSetActiveColumns(pjd->W,bupdated,nv));
    1253        1084 :     for (k=0;k<pep->nmat;k++) {
    1254         820 :       PetscCall(BVSetActiveColumns(pjd->TV[k],bupdated,nv));
    1255         820 :       PetscCall(DSGetMat(pep->ds,DSMatExtra[k],&G));
    1256         820 :       PetscCall(BVMatProject(pjd->TV[k],NULL,pjd->W,G));
    1257         820 :       PetscCall(DSRestoreMat(pep->ds,DSMatExtra[k],&G));
    1258             :     }
    1259         264 :     PetscCall(BVSetActiveColumns(pjd->V,0,nv));
    1260         264 :     PetscCall(BVSetActiveColumns(pjd->W,0,nv));
    1261             : 
    1262             :     /* Solve projected problem */
    1263         264 :     PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
    1264         264 :     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
    1265         264 :     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
    1266         264 :     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
    1267             :     idx = 0;
    1268         291 :     do {
    1269         291 :       ritz[0] = pep->eigr[idx];
    1270             : #if !defined(PETSC_USE_COMPLEX)
    1271         291 :       ritz[1] = pep->eigi[idx];
    1272         291 :       sz = (ritz[1]==0.0)?1:2;
    1273             : #endif
    1274             :       /* Compute Ritz vector u=V*X(:,1) */
    1275         291 :       PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
    1276         291 :       PetscCall(BVSetActiveColumns(pjd->V,0,nv));
    1277         291 :       PetscCall(BVMultVec(pjd->V,1.0,0.0,u[0],pX+idx*ld));
    1278             : #if !defined(PETSC_USE_COMPLEX)
    1279         291 :       if (sz==2) PetscCall(BVMultVec(pjd->V,1.0,0.0,u[1],pX+(idx+1)*ld));
    1280             : #endif
    1281         291 :       PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
    1282         291 :       PetscCall(PEPJDComputeResidual(pep,PETSC_FALSE,sz,u,ritz,r,ww));
    1283             :       /* Check convergence */
    1284         291 :       PetscCall(VecNorm(r[0],NORM_2,&norm));
    1285             : #if !defined(PETSC_USE_COMPLEX)
    1286         291 :       if (sz==2) {
    1287          33 :         PetscCall(VecNorm(r[1],NORM_2,&norm1));
    1288          33 :         norm = SlepcAbs(norm,norm1);
    1289             :       }
    1290             : #endif
    1291         291 :       PetscCall((*pep->converged)(pep,ritz[0],ritz[1],norm,&pep->errest[pep->nconv],pep->convergedctx));
    1292         291 :       if (sz==2) pep->errest[pep->nconv+1] = pep->errest[pep->nconv];
    1293         291 :       if (ini) {
    1294          40 :         tol = PetscMin(.1,pep->errest[pep->nconv]); ini = PETSC_FALSE;
    1295         251 :       } else tol = PetscMin(pep->errest[pep->nconv],tol/2);
    1296         291 :       PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,(pep->errest[pep->nconv]<pep->tol)?pep->nconv+sz:pep->nconv,pep->nev,&pep->reason,pep->stoppingctx));
    1297         291 :       if (pep->errest[pep->nconv]<pep->tol) {
    1298             :         /* Ritz pair converged */
    1299          29 :         ini = PETSC_TRUE;
    1300          29 :         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
    1301          29 :         if (pjd->ld>1) {
    1302          29 :           PetscCall(BVGetColumn(pjd->X,pep->nconv,&v[0]));
    1303          29 :           PetscCall(PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*pep->nconv,pjd->ld-1,0,u[0],PETSC_TRUE));
    1304          29 :           PetscCall(BVRestoreColumn(pjd->X,pep->nconv,&v[0]));
    1305          29 :           PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv+1));
    1306          29 :           PetscCall(BVNormColumn(pjd->X,pep->nconv,NORM_2,&norm));
    1307          29 :           PetscCall(BVScaleColumn(pjd->X,pep->nconv,1.0/norm));
    1308          63 :           for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*pep->nconv+k] *= PetscSqrtReal(np)/norm;
    1309          29 :           pjd->T[(pep->ncv+1)*pep->nconv] = ritz[0];
    1310          29 :           eig[pep->nconv] = ritz[0];
    1311          29 :           idx++;
    1312             : #if !defined(PETSC_USE_COMPLEX)
    1313          29 :           if (sz==2) {
    1314           4 :             PetscCall(BVGetColumn(pjd->X,pep->nconv+1,&v[0]));
    1315           4 :             PetscCall(PEPJDCopyToExtendedVec(pep,v[0],pjd->T+pep->ncv*(pep->nconv+1),pjd->ld-1,0,u[1],PETSC_TRUE));
    1316           4 :             PetscCall(BVRestoreColumn(pjd->X,pep->nconv+1,&v[0]));
    1317           4 :             PetscCall(BVSetActiveColumns(pjd->X,0,pep->nconv+2));
    1318           4 :             PetscCall(BVNormColumn(pjd->X,pep->nconv+1,NORM_2,&norm1));
    1319           4 :             PetscCall(BVScaleColumn(pjd->X,pep->nconv+1,1.0/norm1));
    1320           6 :             for (k=0;k<pep->nconv;k++) pjd->T[pep->ncv*(pep->nconv+1)+k] *= PetscSqrtReal(np)/norm1;
    1321           4 :             pjd->T[(pep->ncv+1)*(pep->nconv+1)] = ritz[0];
    1322           4 :             pjd->T[(pep->ncv+1)*pep->nconv+1] = -ritz[1]*norm1/norm;
    1323           4 :             pjd->T[(pep->ncv+1)*(pep->nconv+1)-1] = ritz[1]*norm/norm1;
    1324           4 :             eig[pep->nconv+1] = ritz[0];
    1325           4 :             eigi[pep->nconv] = ritz[1]; eigi[pep->nconv+1] = -ritz[1];
    1326           4 :             idx++;
    1327             :           }
    1328             : #endif
    1329           0 :         } else PetscCall(BVInsertVec(pep->V,pep->nconv,u[0]));
    1330          29 :         pep->nconv += sz;
    1331             :       }
    1332         291 :     } while (pep->errest[pep->nconv]<pep->tol && pep->nconv<nv);
    1333             : 
    1334         264 :     if (pep->reason==PEP_CONVERGED_ITERATING) {
    1335         253 :       nvc = nv;
    1336         253 :       if (idx) {
    1337          15 :         pjd->nlock +=idx;
    1338          15 :         PetscCall(PEPJDLockConverged(pep,&nv,idx));
    1339             :       }
    1340         253 :       if (nv+sz>=pep->ncv-1) {
    1341             :         /* Basis full, force restart */
    1342          16 :         minv = PetscMin(nv,(PetscInt)(pjd->keep*pep->ncv));
    1343          16 :         PetscCall(DSGetDimensions(pep->ds,&dim,NULL,NULL,NULL));
    1344          16 :         PetscCall(DSGetArray(pep->ds,DS_MAT_X,&pX));
    1345          16 :         PetscCall(PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld));
    1346          16 :         PetscCall(DSRestoreArray(pep->ds,DS_MAT_X,&pX));
    1347          16 :         PetscCall(DSGetArray(pep->ds,DS_MAT_Y,&pX));
    1348          16 :         PetscCall(PEPJDOrthogonalize(dim,minv,pX,ld,&minv,NULL,NULL,ld));
    1349          16 :         PetscCall(DSRestoreArray(pep->ds,DS_MAT_Y,&pX));
    1350          16 :         PetscCall(DSGetMat(pep->ds,DS_MAT_X,&X));
    1351          16 :         PetscCall(BVMultInPlace(pjd->V,X,0,minv));
    1352          16 :         PetscCall(DSRestoreMat(pep->ds,DS_MAT_X,&X));
    1353          16 :         if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
    1354          12 :          PetscCall(DSGetMat(pep->ds,DS_MAT_Y,&Y));
    1355          12 :          PetscCall(BVMultInPlace(pjd->W,Y,pep->nconv,minv));
    1356          12 :          PetscCall(DSRestoreMat(pep->ds,DS_MAT_Y,&Y));
    1357             :         }
    1358          16 :         nv = minv;
    1359          16 :         bupdated = 0;
    1360             :       } else {
    1361         237 :         if (!idx && pep->errest[pep->nconv]<pjd->fix) {theta[0] = ritz[0]; theta[1] = ritz[1];}
    1362         153 :         else {theta[0] = pep->target; theta[1] = 0.0;}
    1363             :         /* Update system mat */
    1364         237 :         PetscCall(PEPJDSystemSetUp(pep,sz,theta,u,p,ww));
    1365             :         /* Solve correction equation to expand basis */
    1366         237 :         PetscCall(BVGetColumn(pjd->V,nv,&t[0]));
    1367         237 :         rr[0] = r[0];
    1368         237 :         if (sz==2) {
    1369          24 :           PetscCall(BVGetColumn(pjd->V,nv+1,&t[1]));
    1370          24 :           rr[1] = r[1];
    1371             :         } else {
    1372         213 :           t[1] = NULL;
    1373         213 :           rr[1] = NULL;
    1374             :         }
    1375         237 :         PetscCall(VecCreateCompWithVecs(t,kspsf,pjd->vtempl,&tc));
    1376         237 :         PetscCall(VecCreateCompWithVecs(rr,kspsf,pjd->vtempl,&rc));
    1377         237 :         PetscCall(VecCompSetSubVecs(pjd->vtempl,sz,NULL));
    1378         237 :         tol  = PetscMax(rtol,tol/2);
    1379         237 :         PetscCall(KSPSetTolerances(ksp,tol,abstol,dtol,maxits));
    1380         237 :         PetscCall(KSPSolve(ksp,rc,tc));
    1381         237 :         PetscCall(VecDestroy(&tc));
    1382         237 :         PetscCall(VecDestroy(&rc));
    1383         237 :         PetscCall(VecGetArray(t[0],&array));
    1384         237 :         PetscCall(PetscMPIIntCast(pep->nconv,&count));
    1385         474 :         PetscCallMPI(MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
    1386         237 :         PetscCall(VecRestoreArray(t[0],&array));
    1387         237 :         PetscCall(BVRestoreColumn(pjd->V,nv,&t[0]));
    1388         237 :         PetscCall(BVOrthogonalizeColumn(pjd->V,nv,NULL,&norm,&lindep));
    1389         237 :         if (lindep || norm==0.0) {
    1390           0 :           PetscCheck(sz!=1,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
    1391             :           off = 1;
    1392             :         } else {
    1393         237 :           off = 0;
    1394         237 :           PetscCall(BVScaleColumn(pjd->V,nv,1.0/norm));
    1395             :         }
    1396             : #if !defined(PETSC_USE_COMPLEX)
    1397         237 :         if (sz==2) {
    1398          24 :           PetscCall(VecGetArray(t[1],&array));
    1399          48 :           PetscCallMPI(MPI_Bcast(array+nloc,count,MPIU_SCALAR,np-1,PetscObjectComm((PetscObject)pep)));
    1400          24 :           PetscCall(VecRestoreArray(t[1],&array));
    1401          24 :           PetscCall(BVRestoreColumn(pjd->V,nv+1,&t[1]));
    1402          24 :           if (off) PetscCall(BVCopyColumn(pjd->V,nv+1,nv));
    1403          24 :           PetscCall(BVOrthogonalizeColumn(pjd->V,nv+1-off,NULL,&norm,&lindep));
    1404          24 :           if (lindep || norm==0.0) {
    1405           0 :             PetscCheck(off==0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
    1406             :             off = 1;
    1407          24 :           } else PetscCall(BVScaleColumn(pjd->V,nv+1-off,1.0/norm));
    1408             :         }
    1409             : #endif
    1410         237 :         if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) {
    1411         204 :           PetscCall(BVInsertVec(pjd->W,nv,r[0]));
    1412         204 :           if (sz==2 && !off) PetscCall(BVInsertVec(pjd->W,nv+1,r[1]));
    1413         204 :           PetscCall(BVOrthogonalizeColumn(pjd->W,nv,NULL,&norm,&lindep));
    1414         204 :           PetscCheck(!lindep && norm>0.0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
    1415         204 :           PetscCall(BVScaleColumn(pjd->W,nv,1.0/norm));
    1416         204 :           if (sz==2 && !off) {
    1417          18 :             PetscCall(BVOrthogonalizeColumn(pjd->W,nv+1,NULL,&norm,&lindep));
    1418          18 :             PetscCheck(!lindep && norm>0.0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linearly dependent continuation vector");
    1419          18 :             PetscCall(BVScaleColumn(pjd->W,nv+1,1.0/norm));
    1420             :           }
    1421             :         }
    1422         237 :         bupdated = idx?0:nv;
    1423         237 :         nv += sz-off;
    1424             :       }
    1425        2745 :       for (k=0;k<nvc;k++) {
    1426        2492 :         eig[pep->nconv-idx+k] = pep->eigr[k];
    1427             : #if !defined(PETSC_USE_COMPLEX)
    1428        2492 :         eigi[pep->nconv-idx+k] = pep->eigi[k];
    1429             : #endif
    1430             :       }
    1431         528 :       PetscCall(PEPMonitor(pep,pep->its,pep->nconv,eig,eigi,pep->errest,pep->nconv+1));
    1432             :     }
    1433             :   }
    1434          11 :   if (pjd->ld>1) {
    1435          44 :     for (k=0;k<pep->nconv;k++) {
    1436          33 :       pep->eigr[k] = eig[k];
    1437          33 :       pep->eigi[k] = eigi[k];
    1438             :     }
    1439          11 :     if (pep->nconv>0) PetscCall(PEPJDEigenvectors(pep));
    1440          11 :     PetscCall(PetscFree2(pcctx->M,pcctx->ps));
    1441             :   }
    1442          11 :   PetscCall(VecDestroy(&u[0]));
    1443          11 :   PetscCall(VecDestroy(&r[0]));
    1444          11 :   PetscCall(VecDestroy(&p[0]));
    1445             : #if !defined (PETSC_USE_COMPLEX)
    1446          11 :   PetscCall(VecDestroy(&u[1]));
    1447          11 :   PetscCall(VecDestroy(&r[1]));
    1448          11 :   PetscCall(VecDestroy(&p[1]));
    1449             : #endif
    1450          11 :   PetscCall(KSPSetTolerances(ksp,rtol,abstol,dtol,maxits));
    1451          11 :   PetscCall(KSPSetPC(ksp,pcctx->pc));
    1452          11 :   PetscCall(VecDestroy(&pcctx->Bp[0]));
    1453          11 :   PetscCall(VecDestroy(&pcctx->Bp[1]));
    1454          11 :   PetscCall(MatShellGetContext(pjd->Pshell,&matctx));
    1455          11 :   PetscCall(MatDestroy(&matctx->Pr));
    1456          11 :   PetscCall(MatDestroy(&matctx->Pi));
    1457          11 :   PetscCall(MatDestroy(&pjd->Pshell));
    1458          11 :   PetscCall(MatDestroy(&pcctx->PPr));
    1459          11 :   PetscCall(PCDestroy(&pcctx->pc));
    1460          11 :   PetscCall(PetscFree(pcctx));
    1461          11 :   PetscCall(PetscFree(matctx));
    1462          11 :   PetscCall(PCDestroy(&pjd->pcshell));
    1463          11 :   PetscCall(PetscFree3(eig,eigi,res));
    1464          11 :   PetscCall(VecDestroy(&pjd->vtempl));
    1465          11 :   PetscFunctionReturn(PETSC_SUCCESS);
    1466             : }
    1467             : 
    1468           3 : static PetscErrorCode PEPJDSetRestart_JD(PEP pep,PetscReal keep)
    1469             : {
    1470           3 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1471             : 
    1472           3 :   PetscFunctionBegin;
    1473           3 :   if (keep==(PetscReal)PETSC_DEFAULT) pjd->keep = 0.5;
    1474             :   else {
    1475           3 :     PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
    1476           3 :     pjd->keep = keep;
    1477             :   }
    1478           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1479             : }
    1480             : 
    1481             : /*@
    1482             :    PEPJDSetRestart - Sets the restart parameter for the Jacobi-Davidson
    1483             :    method, in particular the proportion of basis vectors that must be kept
    1484             :    after restart.
    1485             : 
    1486             :    Logically Collective
    1487             : 
    1488             :    Input Parameters:
    1489             : +  pep  - the eigenproblem solver context
    1490             : -  keep - the number of vectors to be kept at restart
    1491             : 
    1492             :    Options Database Key:
    1493             : .  -pep_jd_restart - Sets the restart parameter
    1494             : 
    1495             :    Notes:
    1496             :    Allowed values are in the range [0.1,0.9]. The default is 0.5.
    1497             : 
    1498             :    Level: advanced
    1499             : 
    1500             : .seealso: PEPJDGetRestart()
    1501             : @*/
    1502           3 : PetscErrorCode PEPJDSetRestart(PEP pep,PetscReal keep)
    1503             : {
    1504           3 :   PetscFunctionBegin;
    1505           3 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1506          12 :   PetscValidLogicalCollectiveReal(pep,keep,2);
    1507           3 :   PetscTryMethod(pep,"PEPJDSetRestart_C",(PEP,PetscReal),(pep,keep));
    1508           3 :   PetscFunctionReturn(PETSC_SUCCESS);
    1509             : }
    1510             : 
    1511           2 : static PetscErrorCode PEPJDGetRestart_JD(PEP pep,PetscReal *keep)
    1512             : {
    1513           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1514             : 
    1515           2 :   PetscFunctionBegin;
    1516           2 :   *keep = pjd->keep;
    1517           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1518             : }
    1519             : 
    1520             : /*@
    1521             :    PEPJDGetRestart - Gets the restart parameter used in the Jacobi-Davidson method.
    1522             : 
    1523             :    Not Collective
    1524             : 
    1525             :    Input Parameter:
    1526             : .  pep - the eigenproblem solver context
    1527             : 
    1528             :    Output Parameter:
    1529             : .  keep - the restart parameter
    1530             : 
    1531             :    Level: advanced
    1532             : 
    1533             : .seealso: PEPJDSetRestart()
    1534             : @*/
    1535           2 : PetscErrorCode PEPJDGetRestart(PEP pep,PetscReal *keep)
    1536             : {
    1537           2 :   PetscFunctionBegin;
    1538           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1539           2 :   PetscAssertPointer(keep,2);
    1540           2 :   PetscUseMethod(pep,"PEPJDGetRestart_C",(PEP,PetscReal*),(pep,keep));
    1541           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1542             : }
    1543             : 
    1544           2 : static PetscErrorCode PEPJDSetFix_JD(PEP pep,PetscReal fix)
    1545             : {
    1546           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1547             : 
    1548           2 :   PetscFunctionBegin;
    1549           2 :   if (fix == (PetscReal)PETSC_DEFAULT || fix == (PetscReal)PETSC_DECIDE) pjd->fix = 0.01;
    1550             :   else {
    1551           2 :     PetscCheck(fix>=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid fix value, must be >0");
    1552           2 :     pjd->fix = fix;
    1553             :   }
    1554           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1555             : }
    1556             : 
    1557             : /*@
    1558             :    PEPJDSetFix - Sets the threshold for changing the target in the correction
    1559             :    equation.
    1560             : 
    1561             :    Logically Collective
    1562             : 
    1563             :    Input Parameters:
    1564             : +  pep - the eigenproblem solver context
    1565             : -  fix - threshold for changing the target
    1566             : 
    1567             :    Options Database Key:
    1568             : .  -pep_jd_fix - the fix value
    1569             : 
    1570             :    Note:
    1571             :    The target in the correction equation is fixed at the first iterations.
    1572             :    When the norm of the residual vector is lower than the fix value,
    1573             :    the target is set to the corresponding eigenvalue.
    1574             : 
    1575             :    Level: advanced
    1576             : 
    1577             : .seealso: PEPJDGetFix()
    1578             : @*/
    1579           2 : PetscErrorCode PEPJDSetFix(PEP pep,PetscReal fix)
    1580             : {
    1581           2 :   PetscFunctionBegin;
    1582           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1583           8 :   PetscValidLogicalCollectiveReal(pep,fix,2);
    1584           2 :   PetscTryMethod(pep,"PEPJDSetFix_C",(PEP,PetscReal),(pep,fix));
    1585           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1586             : }
    1587             : 
    1588           2 : static PetscErrorCode PEPJDGetFix_JD(PEP pep,PetscReal *fix)
    1589             : {
    1590           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1591             : 
    1592           2 :   PetscFunctionBegin;
    1593           2 :   *fix = pjd->fix;
    1594           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1595             : }
    1596             : 
    1597             : /*@
    1598             :    PEPJDGetFix - Returns the threshold for changing the target in the correction
    1599             :    equation.
    1600             : 
    1601             :    Not Collective
    1602             : 
    1603             :    Input Parameter:
    1604             : .  pep - the eigenproblem solver context
    1605             : 
    1606             :    Output Parameter:
    1607             : .  fix - threshold for changing the target
    1608             : 
    1609             :    Note:
    1610             :    The target in the correction equation is fixed at the first iterations.
    1611             :    When the norm of the residual vector is lower than the fix value,
    1612             :    the target is set to the corresponding eigenvalue.
    1613             : 
    1614             :    Level: advanced
    1615             : 
    1616             : .seealso: PEPJDSetFix()
    1617             : @*/
    1618           2 : PetscErrorCode PEPJDGetFix(PEP pep,PetscReal *fix)
    1619             : {
    1620           2 :   PetscFunctionBegin;
    1621           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1622           2 :   PetscAssertPointer(fix,2);
    1623           2 :   PetscUseMethod(pep,"PEPJDGetFix_C",(PEP,PetscReal*),(pep,fix));
    1624           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1625             : }
    1626             : 
    1627           1 : static PetscErrorCode PEPJDSetReusePreconditioner_JD(PEP pep,PetscBool reusepc)
    1628             : {
    1629           1 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1630             : 
    1631           1 :   PetscFunctionBegin;
    1632           1 :   pjd->reusepc = reusepc;
    1633           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1634             : }
    1635             : 
    1636             : /*@
    1637             :    PEPJDSetReusePreconditioner - Sets a flag indicating whether the preconditioner
    1638             :    must be reused or not.
    1639             : 
    1640             :    Logically Collective
    1641             : 
    1642             :    Input Parameters:
    1643             : +  pep     - the eigenproblem solver context
    1644             : -  reusepc - the reuse flag
    1645             : 
    1646             :    Options Database Key:
    1647             : .  -pep_jd_reuse_preconditioner - the reuse flag
    1648             : 
    1649             :    Note:
    1650             :    The default value is False. If set to True, the preconditioner is built
    1651             :    only at the beginning, using the target value. Otherwise, it may be rebuilt
    1652             :    (depending on the fix parameter) at each iteration from the Ritz value.
    1653             : 
    1654             :    Level: advanced
    1655             : 
    1656             : .seealso: PEPJDGetReusePreconditioner(), PEPJDSetFix()
    1657             : @*/
    1658           1 : PetscErrorCode PEPJDSetReusePreconditioner(PEP pep,PetscBool reusepc)
    1659             : {
    1660           1 :   PetscFunctionBegin;
    1661           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1662           4 :   PetscValidLogicalCollectiveBool(pep,reusepc,2);
    1663           1 :   PetscTryMethod(pep,"PEPJDSetReusePreconditioner_C",(PEP,PetscBool),(pep,reusepc));
    1664           1 :   PetscFunctionReturn(PETSC_SUCCESS);
    1665             : }
    1666             : 
    1667           2 : static PetscErrorCode PEPJDGetReusePreconditioner_JD(PEP pep,PetscBool *reusepc)
    1668             : {
    1669           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1670             : 
    1671           2 :   PetscFunctionBegin;
    1672           2 :   *reusepc = pjd->reusepc;
    1673           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1674             : }
    1675             : 
    1676             : /*@
    1677             :    PEPJDGetReusePreconditioner - Returns the flag for reusing the preconditioner.
    1678             : 
    1679             :    Not Collective
    1680             : 
    1681             :    Input Parameter:
    1682             : .  pep - the eigenproblem solver context
    1683             : 
    1684             :    Output Parameter:
    1685             : .  reusepc - the reuse flag
    1686             : 
    1687             :    Level: advanced
    1688             : 
    1689             : .seealso: PEPJDSetReusePreconditioner()
    1690             : @*/
    1691           2 : PetscErrorCode PEPJDGetReusePreconditioner(PEP pep,PetscBool *reusepc)
    1692             : {
    1693           2 :   PetscFunctionBegin;
    1694           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1695           2 :   PetscAssertPointer(reusepc,2);
    1696           2 :   PetscUseMethod(pep,"PEPJDGetReusePreconditioner_C",(PEP,PetscBool*),(pep,reusepc));
    1697           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1698             : }
    1699             : 
    1700           2 : static PetscErrorCode PEPJDSetMinimalityIndex_JD(PEP pep,PetscInt mmidx)
    1701             : {
    1702           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1703             : 
    1704           2 :   PetscFunctionBegin;
    1705           2 :   if (mmidx == PETSC_DEFAULT || mmidx == PETSC_DECIDE) {
    1706           0 :     if (pjd->mmidx != 1) pep->state = PEP_STATE_INITIAL;
    1707           0 :     pjd->mmidx = 1;
    1708             :   } else {
    1709           2 :     PetscCheck(mmidx>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mmidx value, should be >0");
    1710           2 :     if (pjd->mmidx != mmidx) pep->state = PEP_STATE_INITIAL;
    1711           2 :     pjd->mmidx = mmidx;
    1712             :   }
    1713           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1714             : }
    1715             : 
    1716             : /*@
    1717             :    PEPJDSetMinimalityIndex - Sets the maximum allowed value for the minimality index.
    1718             : 
    1719             :    Logically Collective
    1720             : 
    1721             :    Input Parameters:
    1722             : +  pep   - the eigenproblem solver context
    1723             : -  mmidx - maximum minimality index
    1724             : 
    1725             :    Options Database Key:
    1726             : .  -pep_jd_minimality_index - the minimality index value
    1727             : 
    1728             :    Note:
    1729             :    The default value is equal to the degree of the polynomial. A smaller value
    1730             :    can be used if the wanted eigenvectors are known to be linearly independent.
    1731             : 
    1732             :    Level: advanced
    1733             : 
    1734             : .seealso: PEPJDGetMinimalityIndex()
    1735             : @*/
    1736           2 : PetscErrorCode PEPJDSetMinimalityIndex(PEP pep,PetscInt mmidx)
    1737             : {
    1738           2 :   PetscFunctionBegin;
    1739           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1740           8 :   PetscValidLogicalCollectiveInt(pep,mmidx,2);
    1741           2 :   PetscTryMethod(pep,"PEPJDSetMinimalityIndex_C",(PEP,PetscInt),(pep,mmidx));
    1742           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1743             : }
    1744             : 
    1745           2 : static PetscErrorCode PEPJDGetMinimalityIndex_JD(PEP pep,PetscInt *mmidx)
    1746             : {
    1747           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1748             : 
    1749           2 :   PetscFunctionBegin;
    1750           2 :   *mmidx = pjd->mmidx;
    1751           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1752             : }
    1753             : 
    1754             : /*@
    1755             :    PEPJDGetMinimalityIndex - Returns the maximum allowed value of the minimality
    1756             :    index.
    1757             : 
    1758             :    Not Collective
    1759             : 
    1760             :    Input Parameter:
    1761             : .  pep - the eigenproblem solver context
    1762             : 
    1763             :    Output Parameter:
    1764             : .  mmidx - minimality index
    1765             : 
    1766             :    Level: advanced
    1767             : 
    1768             : .seealso: PEPJDSetMinimalityIndex()
    1769             : @*/
    1770           2 : PetscErrorCode PEPJDGetMinimalityIndex(PEP pep,PetscInt *mmidx)
    1771             : {
    1772           2 :   PetscFunctionBegin;
    1773           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1774           2 :   PetscAssertPointer(mmidx,2);
    1775           2 :   PetscUseMethod(pep,"PEPJDGetMinimalityIndex_C",(PEP,PetscInt*),(pep,mmidx));
    1776           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1777             : }
    1778             : 
    1779           2 : static PetscErrorCode PEPJDSetProjection_JD(PEP pep,PEPJDProjection proj)
    1780             : {
    1781           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1782             : 
    1783           2 :   PetscFunctionBegin;
    1784           2 :   switch (proj) {
    1785           2 :     case PEP_JD_PROJECTION_HARMONIC:
    1786             :     case PEP_JD_PROJECTION_ORTHOGONAL:
    1787           2 :       if (pjd->proj != proj) {
    1788           2 :         pep->state = PEP_STATE_INITIAL;
    1789           2 :         pjd->proj = proj;
    1790             :       }
    1791           2 :       break;
    1792           0 :     default:
    1793           0 :       SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid 'proj' value");
    1794             :   }
    1795           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1796             : }
    1797             : 
    1798             : /*@
    1799             :    PEPJDSetProjection - Sets the type of projection to be used in the Jacobi-Davidson solver.
    1800             : 
    1801             :    Logically Collective
    1802             : 
    1803             :    Input Parameters:
    1804             : +  pep  - the eigenproblem solver context
    1805             : -  proj - the type of projection
    1806             : 
    1807             :    Options Database Key:
    1808             : .  -pep_jd_projection - the projection type, either orthogonal or harmonic
    1809             : 
    1810             :    Level: advanced
    1811             : 
    1812             : .seealso: PEPJDGetProjection()
    1813             : @*/
    1814           2 : PetscErrorCode PEPJDSetProjection(PEP pep,PEPJDProjection proj)
    1815             : {
    1816           2 :   PetscFunctionBegin;
    1817           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1818           8 :   PetscValidLogicalCollectiveEnum(pep,proj,2);
    1819           2 :   PetscTryMethod(pep,"PEPJDSetProjection_C",(PEP,PEPJDProjection),(pep,proj));
    1820           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1821             : }
    1822             : 
    1823           2 : static PetscErrorCode PEPJDGetProjection_JD(PEP pep,PEPJDProjection *proj)
    1824             : {
    1825           2 :   PEP_JD *pjd = (PEP_JD*)pep->data;
    1826             : 
    1827           2 :   PetscFunctionBegin;
    1828           2 :   *proj = pjd->proj;
    1829           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1830             : }
    1831             : 
    1832             : /*@
    1833             :    PEPJDGetProjection - Returns the type of projection used by the Jacobi-Davidson solver.
    1834             : 
    1835             :    Not Collective
    1836             : 
    1837             :    Input Parameter:
    1838             : .  pep - the eigenproblem solver context
    1839             : 
    1840             :    Output Parameter:
    1841             : .  proj - the type of projection
    1842             : 
    1843             :    Level: advanced
    1844             : 
    1845             : .seealso: PEPJDSetProjection()
    1846             : @*/
    1847           2 : PetscErrorCode PEPJDGetProjection(PEP pep,PEPJDProjection *proj)
    1848             : {
    1849           2 :   PetscFunctionBegin;
    1850           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1851           2 :   PetscAssertPointer(proj,2);
    1852           2 :   PetscUseMethod(pep,"PEPJDGetProjection_C",(PEP,PEPJDProjection*),(pep,proj));
    1853           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1854             : }
    1855             : 
    1856          10 : static PetscErrorCode PEPSetFromOptions_JD(PEP pep,PetscOptionItems *PetscOptionsObject)
    1857             : {
    1858          10 :   PetscBool       flg,b1;
    1859          10 :   PetscReal       r1;
    1860          10 :   PetscInt        i1;
    1861          10 :   PEPJDProjection proj;
    1862             : 
    1863          10 :   PetscFunctionBegin;
    1864          10 :   PetscOptionsHeadBegin(PetscOptionsObject,"PEP JD Options");
    1865             : 
    1866          10 :     PetscCall(PetscOptionsReal("-pep_jd_restart","Proportion of vectors kept after restart","PEPJDSetRestart",0.5,&r1,&flg));
    1867          10 :     if (flg) PetscCall(PEPJDSetRestart(pep,r1));
    1868             : 
    1869          10 :     PetscCall(PetscOptionsReal("-pep_jd_fix","Tolerance for changing the target in the correction equation","PEPJDSetFix",0.01,&r1,&flg));
    1870          10 :     if (flg) PetscCall(PEPJDSetFix(pep,r1));
    1871             : 
    1872          10 :     PetscCall(PetscOptionsBool("-pep_jd_reuse_preconditioner","Whether to reuse the preconditioner","PEPJDSetReusePreconditoiner",PETSC_FALSE,&b1,&flg));
    1873          10 :     if (flg) PetscCall(PEPJDSetReusePreconditioner(pep,b1));
    1874             : 
    1875          10 :     PetscCall(PetscOptionsInt("-pep_jd_minimality_index","Maximum allowed minimality index","PEPJDSetMinimalityIndex",1,&i1,&flg));
    1876          10 :     if (flg) PetscCall(PEPJDSetMinimalityIndex(pep,i1));
    1877             : 
    1878          10 :     PetscCall(PetscOptionsEnum("-pep_jd_projection","Type of projection","PEPJDSetProjection",PEPJDProjectionTypes,(PetscEnum)PEP_JD_PROJECTION_HARMONIC,(PetscEnum*)&proj,&flg));
    1879          10 :     if (flg) PetscCall(PEPJDSetProjection(pep,proj));
    1880             : 
    1881          10 :   PetscOptionsHeadEnd();
    1882          10 :   PetscFunctionReturn(PETSC_SUCCESS);
    1883             : }
    1884             : 
    1885           0 : static PetscErrorCode PEPView_JD(PEP pep,PetscViewer viewer)
    1886             : {
    1887           0 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
    1888           0 :   PetscBool      isascii;
    1889             : 
    1890           0 :   PetscFunctionBegin;
    1891           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    1892           0 :   if (isascii) {
    1893           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*pjd->keep)));
    1894           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  threshold for changing the target in the correction equation (fix): %g\n",(double)pjd->fix));
    1895           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  projection type: %s\n",PEPJDProjectionTypes[pjd->proj]));
    1896           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum allowed minimality index: %" PetscInt_FMT "\n",pjd->mmidx));
    1897           0 :     if (pjd->reusepc) PetscCall(PetscViewerASCIIPrintf(viewer,"  reusing the preconditioner\n"));
    1898             :   }
    1899           0 :   PetscFunctionReturn(PETSC_SUCCESS);
    1900             : }
    1901             : 
    1902          21 : static PetscErrorCode PEPSetDefaultST_JD(PEP pep)
    1903             : {
    1904          21 :   KSP            ksp;
    1905             : 
    1906          21 :   PetscFunctionBegin;
    1907          21 :   if (!((PetscObject)pep->st)->type_name) {
    1908          10 :     PetscCall(STSetType(pep->st,STPRECOND));
    1909          10 :     PetscCall(STPrecondSetKSPHasMat(pep->st,PETSC_TRUE));
    1910             :   }
    1911          21 :   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
    1912          21 :   PetscCall(STGetKSP(pep->st,&ksp));
    1913          21 :   if (!((PetscObject)ksp)->type_name) {
    1914          10 :     PetscCall(KSPSetType(ksp,KSPBCGSL));
    1915          10 :     PetscCall(KSPSetTolerances(ksp,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,100));
    1916             :   }
    1917          21 :   PetscFunctionReturn(PETSC_SUCCESS);
    1918             : }
    1919             : 
    1920          12 : static PetscErrorCode PEPReset_JD(PEP pep)
    1921             : {
    1922          12 :   PEP_JD         *pjd = (PEP_JD*)pep->data;
    1923          12 :   PetscInt       i;
    1924             : 
    1925          12 :   PetscFunctionBegin;
    1926          47 :   for (i=0;i<pep->nmat;i++) PetscCall(BVDestroy(pjd->TV+i));
    1927          12 :   if (pjd->proj==PEP_JD_PROJECTION_HARMONIC) PetscCall(BVDestroy(&pjd->W));
    1928          12 :   if (pjd->ld>1) {
    1929          12 :     PetscCall(BVDestroy(&pjd->V));
    1930          47 :     for (i=0;i<pep->nmat;i++) PetscCall(BVDestroy(pjd->AX+i));
    1931          12 :     PetscCall(BVDestroy(&pjd->N[0]));
    1932          12 :     PetscCall(BVDestroy(&pjd->N[1]));
    1933          12 :     PetscCall(PetscFree3(pjd->XpX,pjd->T,pjd->Tj));
    1934             :   }
    1935          12 :   PetscCall(PetscFree2(pjd->TV,pjd->AX));
    1936          12 :   PetscFunctionReturn(PETSC_SUCCESS);
    1937             : }
    1938             : 
    1939          10 : static PetscErrorCode PEPDestroy_JD(PEP pep)
    1940             : {
    1941          10 :   PetscFunctionBegin;
    1942          10 :   PetscCall(PetscFree(pep->data));
    1943          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",NULL));
    1944          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",NULL));
    1945          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",NULL));
    1946          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",NULL));
    1947          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",NULL));
    1948          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",NULL));
    1949          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",NULL));
    1950          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",NULL));
    1951          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",NULL));
    1952          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",NULL));
    1953          10 :   PetscFunctionReturn(PETSC_SUCCESS);
    1954             : }
    1955             : 
    1956          10 : SLEPC_EXTERN PetscErrorCode PEPCreate_JD(PEP pep)
    1957             : {
    1958          10 :   PEP_JD         *pjd;
    1959             : 
    1960          10 :   PetscFunctionBegin;
    1961          10 :   PetscCall(PetscNew(&pjd));
    1962          10 :   pep->data = (void*)pjd;
    1963             : 
    1964          10 :   pep->lineariz = PETSC_FALSE;
    1965          10 :   pjd->fix      = 0.01;
    1966          10 :   pjd->mmidx    = 0;
    1967             : 
    1968          10 :   pep->ops->solve          = PEPSolve_JD;
    1969          10 :   pep->ops->setup          = PEPSetUp_JD;
    1970          10 :   pep->ops->setfromoptions = PEPSetFromOptions_JD;
    1971          10 :   pep->ops->destroy        = PEPDestroy_JD;
    1972          10 :   pep->ops->reset          = PEPReset_JD;
    1973          10 :   pep->ops->view           = PEPView_JD;
    1974          10 :   pep->ops->setdefaultst   = PEPSetDefaultST_JD;
    1975             : 
    1976          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetRestart_C",PEPJDSetRestart_JD));
    1977          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetRestart_C",PEPJDGetRestart_JD));
    1978          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetFix_C",PEPJDSetFix_JD));
    1979          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetFix_C",PEPJDGetFix_JD));
    1980          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetReusePreconditioner_C",PEPJDSetReusePreconditioner_JD));
    1981          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetReusePreconditioner_C",PEPJDGetReusePreconditioner_JD));
    1982          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetMinimalityIndex_C",PEPJDSetMinimalityIndex_JD));
    1983          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetMinimalityIndex_C",PEPJDGetMinimalityIndex_JD));
    1984          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDSetProjection_C",PEPJDSetProjection_JD));
    1985          10 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPJDGetProjection_C",PEPJDGetProjection_JD));
    1986          10 :   PetscFunctionReturn(PETSC_SUCCESS);
    1987             : }

Generated by: LCOV version 1.14