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 quadratic eigensolver: "qarnoldi"
      12             : 
      13             :    Method: Q-Arnoldi
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Quadratic Arnoldi with Krylov-Schur type restart.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] K. Meerbergen, "The Quadratic Arnoldi method for the solution
      22             :            of the quadratic eigenvalue problem", SIAM J. Matrix Anal.
      23             :            Appl. 30(4):1462-1482, 2008.
      24             : */
      25             : 
      26             : #include <slepc/private/pepimpl.h>    /*I "slepcpep.h" I*/
      27             : #include <petscblaslapack.h>
      28             : 
      29             : typedef struct {
      30             :   PetscReal keep;         /* restart parameter */
      31             :   PetscBool lock;         /* locking/non-locking variant */
      32             : } PEP_QARNOLDI;
      33             : 
      34          16 : static PetscErrorCode PEPSetUp_QArnoldi(PEP pep)
      35             : {
      36          16 :   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
      37          16 :   PetscBool      flg;
      38             : 
      39          16 :   PetscFunctionBegin;
      40          16 :   PEPCheckQuadratic(pep);
      41          16 :   PEPCheckShiftSinvert(pep);
      42          16 :   PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
      43          16 :   PetscCheck(ctx->lock || pep->mpd>=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
      44          16 :   if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,4*pep->n/pep->ncv);
      45          16 :   if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
      46          16 :   PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
      47             : 
      48          16 :   PetscCall(STGetTransform(pep->st,&flg));
      49          16 :   PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag set, see STSetTransform()");
      50             : 
      51             :   /* set default extraction */
      52          16 :   if (!pep->extract) pep->extract = PEP_EXTRACT_NONE;
      53          16 :   PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_EXTRACT);
      54             : 
      55          16 :   if (!ctx->keep) ctx->keep = 0.5;
      56             : 
      57          16 :   PetscCall(PEPAllocateSolution(pep,0));
      58          16 :   PetscCall(PEPSetWorkVecs(pep,4));
      59             : 
      60          16 :   PetscCall(DSSetType(pep->ds,DSNHEP));
      61          16 :   PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
      62          16 :   PetscCall(DSAllocate(pep->ds,pep->ncv+1));
      63          16 :   PetscFunctionReturn(PETSC_SUCCESS);
      64             : }
      65             : 
      66          16 : static PetscErrorCode PEPExtractVectors_QArnoldi(PEP pep)
      67             : {
      68          16 :   PetscInt       k=pep->nconv;
      69          16 :   Mat            X,X0;
      70             : 
      71          16 :   PetscFunctionBegin;
      72          16 :   if (pep->nconv==0) PetscFunctionReturn(PETSC_SUCCESS);
      73          16 :   PetscCall(DSVectors(pep->ds,DS_MAT_X,NULL,NULL));
      74             : 
      75             :   /* update vectors V = V*X */
      76          16 :   PetscCall(DSGetMat(pep->ds,DS_MAT_X,&X));
      77          16 :   PetscCall(MatDenseGetSubMatrix(X,0,k,0,k,&X0));
      78          16 :   PetscCall(BVMultInPlace(pep->V,X0,0,k));
      79          16 :   PetscCall(MatDenseRestoreSubMatrix(X,&X0));
      80          16 :   PetscCall(DSRestoreMat(pep->ds,DS_MAT_X,&X));
      81          16 :   PetscFunctionReturn(PETSC_SUCCESS);
      82             : }
      83             : 
      84             : /*
      85             :   Compute a step of Classical Gram-Schmidt orthogonalization
      86             : */
      87        2695 : static PetscErrorCode PEPQArnoldiCGS(PEP pep,PetscScalar *H,PetscBLASInt ldh,PetscScalar *h,PetscBLASInt j,BV V,Vec t,Vec v,Vec w,PetscReal *onorm,PetscReal *norm,PetscScalar *work)
      88             : {
      89        2695 :   PetscBLASInt   ione = 1,j_1 = j+1;
      90        2695 :   PetscReal      x,y;
      91        2695 :   PetscScalar    dot,one = 1.0,zero = 0.0;
      92             : 
      93        2695 :   PetscFunctionBegin;
      94             :   /* compute norm of v and w */
      95        2695 :   if (onorm) {
      96        1362 :     PetscCall(VecNorm(v,NORM_2,&x));
      97        1362 :     PetscCall(VecNorm(w,NORM_2,&y));
      98        1362 :     *onorm = SlepcAbs(x,y);
      99             :   }
     100             : 
     101             :   /* orthogonalize: compute h */
     102        2695 :   PetscCall(BVDotVec(V,v,h));
     103        2695 :   PetscCall(BVDotVec(V,w,work));
     104        2695 :   if (j>0) PetscCallBLAS("BLASgemv",BLASgemv_("C",&j_1,&j,&one,H,&ldh,work,&ione,&one,h,&ione));
     105        2695 :   PetscCall(VecDot(w,t,&dot));
     106        2695 :   h[j] += dot;
     107             : 
     108             :   /* orthogonalize: update v and w */
     109        2695 :   PetscCall(BVMultVec(V,-1.0,1.0,v,h));
     110        2695 :   if (j>0) {
     111        2677 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&j_1,&j,&one,H,&ldh,h,&ione,&zero,work,&ione));
     112        2677 :     PetscCall(BVMultVec(V,-1.0,1.0,w,work));
     113             :   }
     114        2695 :   PetscCall(VecAXPY(w,-h[j],t));
     115             : 
     116             :   /* compute norm of v and w */
     117        2695 :   if (norm) {
     118        2566 :     PetscCall(VecNorm(v,NORM_2,&x));
     119        2566 :     PetscCall(VecNorm(w,NORM_2,&y));
     120        2566 :     *norm = SlepcAbs(x,y);
     121             :   }
     122        2695 :   PetscFunctionReturn(PETSC_SUCCESS);
     123             : }
     124             : 
     125             : /*
     126             :   Compute a run of Q-Arnoldi iterations
     127             : */
     128         118 : static PetscErrorCode PEPQArnoldi(PEP pep,Mat A,PetscInt k,PetscInt *M,Vec v,Vec w,PetscReal *beta,PetscBool *breakdown,PetscScalar *work)
     129             : {
     130         118 :   PetscInt           i,j,l,m = *M,ldh;
     131         118 :   PetscBLASInt       jj,ldhh;
     132         118 :   Vec                t = pep->work[2],u = pep->work[3];
     133         118 :   BVOrthogRefineType refinement;
     134         118 :   PetscReal          norm=0.0,onorm,eta;
     135         118 :   PetscScalar        *H,*c = work + m;
     136             : 
     137         118 :   PetscFunctionBegin;
     138         118 :   *beta = 0.0;
     139         118 :   PetscCall(MatDenseGetArray(A,&H));
     140         118 :   PetscCall(MatDenseGetLDA(A,&ldh));
     141         118 :   PetscCall(BVGetOrthogonalization(pep->V,NULL,&refinement,&eta,NULL));
     142         118 :   PetscCall(BVInsertVec(pep->V,k,v));
     143        1539 :   for (j=k;j<m;j++) {
     144             :     /* apply operator */
     145        1421 :     PetscCall(VecCopy(w,t));
     146        1421 :     if (pep->Dr) PetscCall(VecPointwiseMult(v,v,pep->Dr));
     147        1421 :     PetscCall(STMatMult(pep->st,0,v,u));
     148        1421 :     PetscCall(VecCopy(t,v));
     149        1421 :     if (pep->Dr) PetscCall(VecPointwiseMult(t,t,pep->Dr));
     150        1421 :     PetscCall(STMatMult(pep->st,1,t,w));
     151        1421 :     PetscCall(VecAXPY(u,pep->sfactor,w));
     152        1421 :     PetscCall(STMatSolve(pep->st,u,w));
     153        1421 :     PetscCall(VecScale(w,-1.0/(pep->sfactor*pep->sfactor)));
     154        1421 :     if (pep->Dr) PetscCall(VecPointwiseDivide(w,w,pep->Dr));
     155        1421 :     PetscCall(VecCopy(v,t));
     156        1421 :     PetscCall(BVSetActiveColumns(pep->V,0,j+1));
     157             : 
     158             :     /* orthogonalize */
     159        1421 :     PetscCall(PetscBLASIntCast(j,&jj));
     160        1421 :     PetscCall(PetscBLASIntCast(ldh,&ldhh));
     161        1421 :     switch (refinement) {
     162          59 :       case BV_ORTHOG_REFINE_NEVER:
     163          59 :         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,NULL,&norm,work));
     164          59 :         *breakdown = PETSC_FALSE;
     165          59 :         break;
     166         129 :       case BV_ORTHOG_REFINE_ALWAYS:
     167         129 :         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,NULL,NULL,work));
     168         129 :         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,c,jj,pep->V,t,v,w,&onorm,&norm,work));
     169        2033 :         for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
     170         129 :         if (norm < eta * onorm) *breakdown = PETSC_TRUE;
     171         129 :         else *breakdown = PETSC_FALSE;
     172             :         break;
     173        1233 :       case BV_ORTHOG_REFINE_IFNEEDED:
     174        1233 :         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,&onorm,&norm,work));
     175             :         /* ||q|| < eta ||h|| */
     176             :         l = 1;
     177        2378 :         while (l<3 && norm < eta * onorm) {
     178        1145 :           l++;
     179        1145 :           onorm = norm;
     180        1145 :           PetscCall(PEPQArnoldiCGS(pep,H,ldhh,c,jj,pep->V,t,v,w,NULL,&norm,work));
     181       20803 :           for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
     182             :         }
     183        1233 :         if (norm < eta * onorm) *breakdown = PETSC_TRUE;
     184        1233 :         else *breakdown = PETSC_FALSE;
     185             :         break;
     186             :     }
     187        1421 :     PetscCall(VecScale(v,1.0/norm));
     188        1421 :     PetscCall(VecScale(w,1.0/norm));
     189             : 
     190        1421 :     H[j+1+ldh*j] = norm;
     191        1421 :     if (j<m-1) PetscCall(BVInsertVec(pep->V,j+1,v));
     192             :   }
     193         118 :   *beta = norm;
     194         118 :   PetscCall(MatDenseRestoreArray(A,&H));
     195         118 :   PetscFunctionReturn(PETSC_SUCCESS);
     196             : }
     197             : 
     198          16 : static PetscErrorCode PEPSolve_QArnoldi(PEP pep)
     199             : {
     200          16 :   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
     201          16 :   PetscInt       j,k,l,lwork,nv,nconv;
     202          16 :   Vec            v=pep->work[0],w=pep->work[1];
     203          16 :   Mat            Q,S;
     204          16 :   PetscScalar    *work;
     205          16 :   PetscReal      beta,norm,x,y;
     206          16 :   PetscBool      breakdown=PETSC_FALSE,sinv;
     207             : 
     208          16 :   PetscFunctionBegin;
     209          16 :   lwork = 7*pep->ncv;
     210          16 :   PetscCall(PetscMalloc1(lwork,&work));
     211          16 :   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
     212          16 :   PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
     213          16 :   PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
     214             : 
     215             :   /* Get the starting Arnoldi vector */
     216          48 :   for (j=0;j<2;j++) {
     217          32 :     if (j>=pep->nini) PetscCall(BVSetRandomColumn(pep->V,j));
     218             :   }
     219          16 :   PetscCall(BVCopyVec(pep->V,0,v));
     220          16 :   PetscCall(BVCopyVec(pep->V,1,w));
     221          16 :   PetscCall(VecNorm(v,NORM_2,&x));
     222          16 :   PetscCall(VecNorm(w,NORM_2,&y));
     223          16 :   norm = SlepcAbs(x,y);
     224          16 :   PetscCall(VecScale(v,1.0/norm));
     225          16 :   PetscCall(VecScale(w,1.0/norm));
     226             : 
     227             :   /* clean projected matrix (including the extra-arrow) */
     228          16 :   PetscCall(DSSetDimensions(pep->ds,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE));
     229          16 :   PetscCall(DSGetMat(pep->ds,DS_MAT_A,&S));
     230          16 :   PetscCall(MatZeroEntries(S));
     231          16 :   PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&S));
     232             : 
     233             :    /* Restart loop */
     234          16 :   l = 0;
     235          16 :   while (pep->reason == PEP_CONVERGED_ITERATING) {
     236         118 :     pep->its++;
     237             : 
     238             :     /* Compute an nv-step Arnoldi factorization */
     239         118 :     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
     240         118 :     PetscCall(DSGetMat(pep->ds,DS_MAT_A,&S));
     241         118 :     PetscCall(PEPQArnoldi(pep,S,pep->nconv+l,&nv,v,w,&beta,&breakdown,work));
     242         118 :     PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&S));
     243         118 :     PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
     244         118 :     PetscCall(DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
     245         118 :     PetscCall(BVSetActiveColumns(pep->V,pep->nconv,nv));
     246             : 
     247             :     /* Solve projected problem */
     248         118 :     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
     249         118 :     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
     250         118 :     PetscCall(DSUpdateExtraRow(pep->ds));
     251         118 :     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
     252             : 
     253             :     /* Check convergence */
     254         118 :     PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k));
     255         118 :     PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
     256         118 :     nconv = k;
     257             : 
     258             :     /* Update l */
     259         118 :     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
     260             :     else {
     261         102 :       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
     262         102 :       PetscCall(DSGetTruncateSize(pep->ds,k,nv,&l));
     263             :     }
     264         118 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     265         118 :     if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     266             : 
     267         118 :     if (pep->reason == PEP_CONVERGED_ITERATING) {
     268         102 :       if (PetscUnlikely(breakdown)) {
     269             :         /* Stop if breakdown */
     270           0 :         PetscCall(PetscInfo(pep,"Breakdown Quadratic Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
     271           0 :         pep->reason = PEP_DIVERGED_BREAKDOWN;
     272             :       } else {
     273             :         /* Prepare the Rayleigh quotient for restart */
     274         102 :         PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
     275             :       }
     276             :     }
     277             :     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
     278         118 :     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&Q));
     279         118 :     PetscCall(BVMultInPlace(pep->V,Q,pep->nconv,k+l));
     280         118 :     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&Q));
     281             : 
     282         118 :     pep->nconv = k;
     283         134 :     PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
     284             :   }
     285          16 :   PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
     286         107 :   for (j=0;j<pep->nconv;j++) {
     287          91 :     pep->eigr[j] *= pep->sfactor;
     288          91 :     pep->eigi[j] *= pep->sfactor;
     289             :   }
     290             : 
     291          16 :   PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
     292          16 :   PetscCall(RGPopScale(pep->rg));
     293             : 
     294          16 :   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
     295          16 :   PetscCall(PetscFree(work));
     296          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     297             : }
     298             : 
     299           1 : static PetscErrorCode PEPQArnoldiSetRestart_QArnoldi(PEP pep,PetscReal keep)
     300             : {
     301           1 :   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
     302             : 
     303           1 :   PetscFunctionBegin;
     304           1 :   if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
     305             :   else {
     306           1 :     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]");
     307           1 :     ctx->keep = keep;
     308             :   }
     309           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     310             : }
     311             : 
     312             : /*@
     313             :    PEPQArnoldiSetRestart - Sets the restart parameter for the Q-Arnoldi
     314             :    method, in particular the proportion of basis vectors that must be kept
     315             :    after restart.
     316             : 
     317             :    Logically Collective
     318             : 
     319             :    Input Parameters:
     320             : +  pep  - the eigenproblem solver context
     321             : -  keep - the number of vectors to be kept at restart
     322             : 
     323             :    Options Database Key:
     324             : .  -pep_qarnoldi_restart - Sets the restart parameter
     325             : 
     326             :    Notes:
     327             :    Allowed values are in the range [0.1,0.9]. The default is 0.5.
     328             : 
     329             :    Level: advanced
     330             : 
     331             : .seealso: PEPQArnoldiGetRestart()
     332             : @*/
     333           1 : PetscErrorCode PEPQArnoldiSetRestart(PEP pep,PetscReal keep)
     334             : {
     335           1 :   PetscFunctionBegin;
     336           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     337           3 :   PetscValidLogicalCollectiveReal(pep,keep,2);
     338           1 :   PetscTryMethod(pep,"PEPQArnoldiSetRestart_C",(PEP,PetscReal),(pep,keep));
     339           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     340             : }
     341             : 
     342           1 : static PetscErrorCode PEPQArnoldiGetRestart_QArnoldi(PEP pep,PetscReal *keep)
     343             : {
     344           1 :   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
     345             : 
     346           1 :   PetscFunctionBegin;
     347           1 :   *keep = ctx->keep;
     348           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     349             : }
     350             : 
     351             : /*@
     352             :    PEPQArnoldiGetRestart - Gets the restart parameter used in the Q-Arnoldi method.
     353             : 
     354             :    Not Collective
     355             : 
     356             :    Input Parameter:
     357             : .  pep - the eigenproblem solver context
     358             : 
     359             :    Output Parameter:
     360             : .  keep - the restart parameter
     361             : 
     362             :    Level: advanced
     363             : 
     364             : .seealso: PEPQArnoldiSetRestart()
     365             : @*/
     366           1 : PetscErrorCode PEPQArnoldiGetRestart(PEP pep,PetscReal *keep)
     367             : {
     368           1 :   PetscFunctionBegin;
     369           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     370           1 :   PetscAssertPointer(keep,2);
     371           1 :   PetscUseMethod(pep,"PEPQArnoldiGetRestart_C",(PEP,PetscReal*),(pep,keep));
     372           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     373             : }
     374             : 
     375           1 : static PetscErrorCode PEPQArnoldiSetLocking_QArnoldi(PEP pep,PetscBool lock)
     376             : {
     377           1 :   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
     378             : 
     379           1 :   PetscFunctionBegin;
     380           1 :   ctx->lock = lock;
     381           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     382             : }
     383             : 
     384             : /*@
     385             :    PEPQArnoldiSetLocking - Choose between locking and non-locking variants of
     386             :    the Q-Arnoldi method.
     387             : 
     388             :    Logically Collective
     389             : 
     390             :    Input Parameters:
     391             : +  pep  - the eigenproblem solver context
     392             : -  lock - true if the locking variant must be selected
     393             : 
     394             :    Options Database Key:
     395             : .  -pep_qarnoldi_locking - Sets the locking flag
     396             : 
     397             :    Notes:
     398             :    The default is to lock converged eigenpairs when the method restarts.
     399             :    This behaviour can be changed so that all directions are kept in the
     400             :    working subspace even if already converged to working accuracy (the
     401             :    non-locking variant).
     402             : 
     403             :    Level: advanced
     404             : 
     405             : .seealso: PEPQArnoldiGetLocking()
     406             : @*/
     407           1 : PetscErrorCode PEPQArnoldiSetLocking(PEP pep,PetscBool lock)
     408             : {
     409           1 :   PetscFunctionBegin;
     410           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     411           3 :   PetscValidLogicalCollectiveBool(pep,lock,2);
     412           1 :   PetscTryMethod(pep,"PEPQArnoldiSetLocking_C",(PEP,PetscBool),(pep,lock));
     413           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     414             : }
     415             : 
     416           1 : static PetscErrorCode PEPQArnoldiGetLocking_QArnoldi(PEP pep,PetscBool *lock)
     417             : {
     418           1 :   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;
     419             : 
     420           1 :   PetscFunctionBegin;
     421           1 :   *lock = ctx->lock;
     422           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     423             : }
     424             : 
     425             : /*@
     426             :    PEPQArnoldiGetLocking - Gets the locking flag used in the Q-Arnoldi method.
     427             : 
     428             :    Not Collective
     429             : 
     430             :    Input Parameter:
     431             : .  pep - the eigenproblem solver context
     432             : 
     433             :    Output Parameter:
     434             : .  lock - the locking flag
     435             : 
     436             :    Level: advanced
     437             : 
     438             : .seealso: PEPQArnoldiSetLocking()
     439             : @*/
     440           1 : PetscErrorCode PEPQArnoldiGetLocking(PEP pep,PetscBool *lock)
     441             : {
     442           1 :   PetscFunctionBegin;
     443           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     444           1 :   PetscAssertPointer(lock,2);
     445           1 :   PetscUseMethod(pep,"PEPQArnoldiGetLocking_C",(PEP,PetscBool*),(pep,lock));
     446           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     447             : }
     448             : 
     449          14 : static PetscErrorCode PEPSetFromOptions_QArnoldi(PEP pep,PetscOptionItems *PetscOptionsObject)
     450             : {
     451          14 :   PetscBool      flg,lock;
     452          14 :   PetscReal      keep;
     453             : 
     454          14 :   PetscFunctionBegin;
     455          14 :   PetscOptionsHeadBegin(PetscOptionsObject,"PEP Q-Arnoldi Options");
     456             : 
     457          14 :     PetscCall(PetscOptionsReal("-pep_qarnoldi_restart","Proportion of vectors kept after restart","PEPQArnoldiSetRestart",0.5,&keep,&flg));
     458          14 :     if (flg) PetscCall(PEPQArnoldiSetRestart(pep,keep));
     459             : 
     460          14 :     PetscCall(PetscOptionsBool("-pep_qarnoldi_locking","Choose between locking and non-locking variants","PEPQArnoldiSetLocking",PETSC_FALSE,&lock,&flg));
     461          14 :     if (flg) PetscCall(PEPQArnoldiSetLocking(pep,lock));
     462             : 
     463          14 :   PetscOptionsHeadEnd();
     464          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     465             : }
     466             : 
     467           0 : static PetscErrorCode PEPView_QArnoldi(PEP pep,PetscViewer viewer)
     468             : {
     469           0 :   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
     470           0 :   PetscBool      isascii;
     471             : 
     472           0 :   PetscFunctionBegin;
     473           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     474           0 :   if (isascii) {
     475           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
     476           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
     477             :   }
     478           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     479             : }
     480             : 
     481          15 : static PetscErrorCode PEPDestroy_QArnoldi(PEP pep)
     482             : {
     483          15 :   PetscFunctionBegin;
     484          15 :   PetscCall(PetscFree(pep->data));
     485          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",NULL));
     486          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",NULL));
     487          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",NULL));
     488          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",NULL));
     489          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     490             : }
     491             : 
     492          15 : SLEPC_EXTERN PetscErrorCode PEPCreate_QArnoldi(PEP pep)
     493             : {
     494          15 :   PEP_QARNOLDI   *ctx;
     495             : 
     496          15 :   PetscFunctionBegin;
     497          15 :   PetscCall(PetscNew(&ctx));
     498          15 :   pep->data = (void*)ctx;
     499             : 
     500          15 :   pep->lineariz = PETSC_TRUE;
     501          15 :   ctx->lock     = PETSC_TRUE;
     502             : 
     503          15 :   pep->ops->solve          = PEPSolve_QArnoldi;
     504          15 :   pep->ops->setup          = PEPSetUp_QArnoldi;
     505          15 :   pep->ops->setfromoptions = PEPSetFromOptions_QArnoldi;
     506          15 :   pep->ops->destroy        = PEPDestroy_QArnoldi;
     507          15 :   pep->ops->view           = PEPView_QArnoldi;
     508          15 :   pep->ops->backtransform  = PEPBackTransform_Default;
     509          15 :   pep->ops->computevectors = PEPComputeVectors_Default;
     510          15 :   pep->ops->extractvectors = PEPExtractVectors_QArnoldi;
     511          15 :   pep->ops->setdefaultst   = PEPSetDefaultST_Transform;
     512             : 
     513          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",PEPQArnoldiSetRestart_QArnoldi));
     514          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",PEPQArnoldiGetRestart_QArnoldi));
     515          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",PEPQArnoldiSetLocking_QArnoldi));
     516          15 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",PEPQArnoldiGetLocking_QArnoldi));
     517          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     518             : }