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             :    Full basis for the linearization of the rational approximation of non-linear eigenproblems
      12             : */
      13             : 
      14             : #include <slepc/private/nepimpl.h>         /*I "slepcnep.h" I*/
      15             : #include "nleigs.h"
      16             : 
      17          84 : static PetscErrorCode MatMult_FullBasis_Sinvert(Mat M,Vec x,Vec y)
      18             : {
      19          84 :   NEP_NLEIGS        *ctx;
      20          84 :   NEP               nep;
      21          84 :   const PetscScalar *px;
      22          84 :   PetscScalar       *beta,*s,*xi,*t,*py,sigma;
      23          84 :   PetscInt          nmat,d,i,k,m;
      24          84 :   Vec               xx,xxx,yy,yyy,w,ww,www;
      25             : 
      26          84 :   PetscFunctionBegin;
      27          84 :   PetscCall(MatShellGetContext(M,&nep));
      28          84 :   ctx = (NEP_NLEIGS*)nep->data;
      29          84 :   beta = ctx->beta; s = ctx->s; xi = ctx->xi;
      30          84 :   sigma = ctx->shifts[0];
      31          84 :   nmat = ctx->nmat;
      32          84 :   d = nmat-1;
      33          84 :   m = nep->nloc;
      34          84 :   PetscCall(PetscMalloc1(ctx->nmat,&t));
      35          84 :   xx = ctx->w[0]; xxx = ctx->w[1]; yy = ctx->w[2]; yyy=ctx->w[3];
      36          84 :   w = nep->work[0]; ww = nep->work[1]; www = nep->work[2];
      37          84 :   PetscCall(VecGetArrayRead(x,&px));
      38          84 :   PetscCall(VecGetArray(y,&py));
      39          84 :   PetscCall(VecPlaceArray(xx,px+(d-1)*m));
      40          84 :   PetscCall(VecPlaceArray(xxx,px+(d-2)*m));
      41          84 :   PetscCall(VecPlaceArray(yy,py+(d-2)*m));
      42          84 :   PetscCall(VecCopy(xxx,yy));
      43          84 :   PetscCall(VecAXPY(yy,beta[d-1]/xi[d-2],xx));
      44          84 :   PetscCall(VecScale(yy,1.0/(s[d-2]-sigma)));
      45          84 :   PetscCall(VecResetArray(xx));
      46          84 :   PetscCall(VecResetArray(xxx));
      47          84 :   PetscCall(VecResetArray(yy));
      48         168 :   for (i=d-3;i>=0;i--) {
      49          84 :     PetscCall(VecPlaceArray(xx,px+(i+1)*m));
      50          84 :     PetscCall(VecPlaceArray(xxx,px+i*m));
      51          84 :     PetscCall(VecPlaceArray(yy,py+i*m));
      52          84 :     PetscCall(VecPlaceArray(yyy,py+(i+1)*m));
      53          84 :     PetscCall(VecCopy(xxx,yy));
      54          84 :     PetscCall(VecAXPY(yy,beta[i+1]/xi[i],xx));
      55          84 :     PetscCall(VecAXPY(yy,-beta[i+1]*(1.0-sigma/xi[i]),yyy));
      56          84 :     PetscCall(VecScale(yy,1.0/(s[i]-sigma)));
      57          84 :     PetscCall(VecResetArray(xx));
      58          84 :     PetscCall(VecResetArray(xxx));
      59          84 :     PetscCall(VecResetArray(yy));
      60          84 :     PetscCall(VecResetArray(yyy));
      61             :   }
      62          84 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
      63          40 :     PetscCall(VecZeroEntries(w));
      64         160 :     for (k=0;k<nep->nt;k++) {
      65         120 :       PetscCall(VecZeroEntries(ww));
      66         120 :       PetscCall(VecPlaceArray(xx,px+(d-1)*m));
      67         120 :       PetscCall(VecAXPY(ww,-ctx->coeffD[k+nep->nt*d]/beta[d],xx));
      68         120 :       PetscCall(VecResetArray(xx));
      69         360 :       for (i=0;i<d-1;i++) {
      70         240 :         PetscCall(VecPlaceArray(yy,py+i*m));
      71         240 :         PetscCall(VecAXPY(ww,-ctx->coeffD[nep->nt*i+k],yy));
      72         240 :         PetscCall(VecResetArray(yy));
      73             :       }
      74         120 :       PetscCall(MatMult(nep->A[k],ww,www));
      75         120 :       PetscCall(VecAXPY(w,1.0,www));
      76             :     }
      77             :   } else {
      78          44 :     PetscCall(VecPlaceArray(xx,px+(d-1)*m));
      79          44 :     PetscCall(MatMult(ctx->D[d],xx,w));
      80          44 :     PetscCall(VecScale(w,-1.0/beta[d]));
      81          44 :     PetscCall(VecResetArray(xx));
      82         132 :     for (i=0;i<d-1;i++) {
      83          88 :       PetscCall(VecPlaceArray(yy,py+i*m));
      84          88 :       PetscCall(MatMult(ctx->D[i],yy,ww));
      85          88 :       PetscCall(VecResetArray(yy));
      86          88 :       PetscCall(VecAXPY(w,-1.0,ww));
      87             :     }
      88             :   }
      89          84 :   PetscCall(VecPlaceArray(yy,py+(d-1)*m));
      90          84 :   PetscCall(KSPSolve(ctx->ksp[0],w,yy));
      91          84 :   PetscCall(NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t));
      92         252 :   for (i=0;i<d-1;i++) {
      93         168 :     PetscCall(VecPlaceArray(yyy,py+i*m));
      94         168 :     PetscCall(VecAXPY(yyy,t[i],yy));
      95         168 :     PetscCall(VecResetArray(yyy));
      96             :   }
      97          84 :   PetscCall(VecScale(yy,t[d-1]));
      98          84 :   PetscCall(VecResetArray(yy));
      99          84 :   PetscCall(VecRestoreArrayRead(x,&px));
     100          84 :   PetscCall(VecRestoreArray(y,&py));
     101          84 :   PetscCall(PetscFree(t));
     102          84 :   PetscFunctionReturn(PETSC_SUCCESS);
     103             : }
     104             : 
     105          66 : static PetscErrorCode MatMultTranspose_FullBasis_Sinvert(Mat M,Vec x,Vec y)
     106             : {
     107          66 :   NEP_NLEIGS        *ctx;
     108          66 :   NEP               nep;
     109          66 :   const PetscScalar *px;
     110          66 :   PetscScalar       *beta,*s,*xi,*t,*py,sigma;
     111          66 :   PetscInt          nmat,d,i,k,m;
     112          66 :   Vec               xx,yy,yyy,w,z0;
     113             : 
     114          66 :   PetscFunctionBegin;
     115          66 :   PetscCall(MatShellGetContext(M,&nep));
     116          66 :   ctx = (NEP_NLEIGS*)nep->data;
     117          66 :   beta = ctx->beta; s = ctx->s; xi = ctx->xi;
     118          66 :   sigma = ctx->shifts[0];
     119          66 :   nmat = ctx->nmat;
     120          66 :   d = nmat-1;
     121          66 :   m = nep->nloc;
     122          66 :   PetscCall(PetscMalloc1(ctx->nmat,&t));
     123          66 :   xx = ctx->w[0]; yy = ctx->w[1]; yyy=ctx->w[2];
     124          66 :   w = nep->work[0]; z0 = nep->work[1];
     125          66 :   PetscCall(VecGetArrayRead(x,&px));
     126          66 :   PetscCall(VecGetArray(y,&py));
     127          66 :   PetscCall(NEPNLEIGSEvalNRTFunct(nep,d,sigma,t));
     128          66 :   PetscCall(VecPlaceArray(xx,px+(d-1)*m));
     129          66 :   PetscCall(VecCopy(xx,w));
     130          66 :   PetscCall(VecScale(w,t[d-1]));
     131          66 :   PetscCall(VecResetArray(xx));
     132         198 :   for (i=0;i<d-1;i++) {
     133         132 :     PetscCall(VecPlaceArray(xx,px+i*m));
     134         132 :     PetscCall(VecAXPY(w,t[i],xx));
     135         132 :     PetscCall(VecResetArray(xx));
     136             :   }
     137          66 :   PetscCall(KSPSolveTranspose(ctx->ksp[0],w,z0));
     138             : 
     139          66 :   PetscCall(VecPlaceArray(yy,py));
     140          66 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     141          40 :     PetscCall(VecZeroEntries(yy));
     142         160 :     for (k=0;k<nep->nt;k++) {
     143         120 :       PetscCall(MatMult(nep->A[k],z0,w));
     144         120 :       PetscCall(VecAXPY(yy,ctx->coeffD[k],w));
     145             :     }
     146          26 :   } else PetscCall(MatMultTranspose(ctx->D[0],z0,yy));
     147          66 :   PetscCall(VecPlaceArray(xx,px));
     148          66 :   PetscCall(VecAXPY(yy,-1.0,xx));
     149          66 :   PetscCall(VecResetArray(xx));
     150          66 :   PetscCall(VecScale(yy,-1.0/(s[0]-sigma)));
     151          66 :   PetscCall(VecResetArray(yy));
     152         132 :   for (i=2;i<d;i++) {
     153          66 :     PetscCall(VecPlaceArray(yy,py+(i-1)*m));
     154          66 :     if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     155          40 :       PetscCall(VecZeroEntries(yy));
     156         160 :       for (k=0;k<nep->nt;k++) {
     157         120 :         PetscCall(MatMult(nep->A[k],z0,w));
     158         120 :         PetscCall(VecAXPY(yy,ctx->coeffD[k+(i-1)*nep->nt],w));
     159             :       }
     160          26 :     } else PetscCall(MatMultTranspose(ctx->D[i-1],z0,yy));
     161          66 :     PetscCall(VecPlaceArray(yyy,py+(i-2)*m));
     162          66 :     PetscCall(VecAXPY(yy,beta[i-1]*(1.0-sigma/xi[i-2]),yyy));
     163          66 :     PetscCall(VecResetArray(yyy));
     164          66 :     PetscCall(VecPlaceArray(xx,px+(i-1)*m));
     165          66 :     PetscCall(VecAXPY(yy,-1.0,xx));
     166          66 :     PetscCall(VecResetArray(xx));
     167          66 :     PetscCall(VecScale(yy,-1.0/(s[i-1]-sigma)));
     168          66 :     PetscCall(VecResetArray(yy));
     169             :   }
     170          66 :   PetscCall(VecPlaceArray(yy,py+(d-1)*m));
     171          66 :   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
     172          40 :     PetscCall(VecZeroEntries(yy));
     173         160 :     for (k=0;k<nep->nt;k++) {
     174         120 :       PetscCall(MatMult(nep->A[k],z0,w));
     175         120 :       PetscCall(VecAXPY(yy,ctx->coeffD[k+d*nep->nt],w));
     176             :     }
     177          26 :   } else PetscCall(MatMultTranspose(ctx->D[d],z0,yy));
     178          66 :   PetscCall(VecScale(yy,-1.0/beta[d]));
     179          66 :   PetscCall(VecPlaceArray(yyy,py+(d-2)*m));
     180          66 :   PetscCall(VecAXPY(yy,beta[d-1]/xi[d-2],yyy));
     181          66 :   PetscCall(VecResetArray(yyy));
     182          66 :   PetscCall(VecResetArray(yy));
     183             : 
     184         132 :   for (i=d-2;i>0;i--) {
     185          66 :     PetscCall(VecPlaceArray(yyy,py+(i-1)*m));
     186          66 :     PetscCall(VecPlaceArray(yy,py+i*m));
     187          66 :     PetscCall(VecAXPY(yy,beta[i]/xi[i-1],yyy));
     188          66 :     PetscCall(VecResetArray(yyy));
     189          66 :     PetscCall(VecResetArray(yy));
     190             :   }
     191             : 
     192          66 :   PetscCall(VecRestoreArrayRead(x,&px));
     193          66 :   PetscCall(VecRestoreArray(y,&py));
     194          66 :   PetscCall(PetscFree(t));
     195          66 :   PetscFunctionReturn(PETSC_SUCCESS);
     196             : }
     197             : 
     198        1746 : static PetscErrorCode BackTransform_FullBasis(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
     199             : {
     200        1746 :   NEP            nep;
     201             : 
     202        1746 :   PetscFunctionBegin;
     203        1746 :   PetscCall(STShellGetContext(st,&nep));
     204        1746 :   PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,n,eigr,eigi));
     205        1746 :   PetscFunctionReturn(PETSC_SUCCESS);
     206             : }
     207             : 
     208          84 : static PetscErrorCode Apply_FullBasis(ST st,Vec x,Vec y)
     209             : {
     210          84 :   NEP            nep;
     211          84 :   NEP_NLEIGS     *ctx;
     212             : 
     213          84 :   PetscFunctionBegin;
     214          84 :   PetscCall(STShellGetContext(st,&nep));
     215          84 :   ctx = (NEP_NLEIGS*)nep->data;
     216          84 :   PetscCall(MatMult(ctx->A,x,y));
     217          84 :   PetscFunctionReturn(PETSC_SUCCESS);
     218             : }
     219             : 
     220          66 : static PetscErrorCode ApplyTranspose_FullBasis(ST st,Vec x,Vec y)
     221             : {
     222          66 :   NEP            nep;
     223          66 :   NEP_NLEIGS     *ctx;
     224             : 
     225          66 :   PetscFunctionBegin;
     226          66 :   PetscCall(STShellGetContext(st,&nep));
     227          66 :   ctx = (NEP_NLEIGS*)nep->data;
     228          66 :   PetscCall(MatMultTranspose(ctx->A,x,y));
     229          66 :   PetscFunctionReturn(PETSC_SUCCESS);
     230             : }
     231             : 
     232           3 : PetscErrorCode NEPSetUp_NLEIGS_FullBasis(NEP nep)
     233             : {
     234           3 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     235           3 :   ST             st;
     236           3 :   Mat            Q;
     237           3 :   PetscInt       i=0,deg=ctx->nmat-1;
     238           3 :   PetscBool      trackall,istrivial,ks;
     239           3 :   PetscScalar    *epsarray,*neparray;
     240           3 :   Vec            veps,w=NULL;
     241           3 :   EPSWhich       which;
     242             : 
     243           3 :   PetscFunctionBegin;
     244           3 :   PetscCheck(ctx->nshifts==0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The full-basis option is not supported with rational Krylov");
     245           3 :   if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
     246           3 :   PetscCall(EPSGetST(ctx->eps,&st));
     247           3 :   PetscCall(EPSSetTarget(ctx->eps,nep->target));
     248           3 :   PetscCall(STSetDefaultShift(st,nep->target));
     249           3 :   if (!((PetscObject)ctx->eps)->type_name) PetscCall(EPSSetType(ctx->eps,EPSKRYLOVSCHUR));
     250             :   else {
     251           3 :     PetscCall(PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks));
     252           3 :     PetscCheck(ks,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Full-basis option only implemented for Krylov-Schur");
     253             :   }
     254           3 :   PetscCall(STSetType(st,STSHELL));
     255           3 :   PetscCall(STShellSetContext(st,nep));
     256           3 :   PetscCall(STShellSetBackTransform(st,BackTransform_FullBasis));
     257           3 :   PetscCall(KSPGetOperators(ctx->ksp[0],&Q,NULL));
     258           3 :   PetscCall(MatCreateVecsEmpty(Q,&ctx->w[0],&ctx->w[1]));
     259           3 :   PetscCall(MatCreateVecsEmpty(Q,&ctx->w[2],&ctx->w[3]));
     260           3 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),deg*nep->nloc,deg*nep->nloc,deg*nep->n,deg*nep->n,nep,&ctx->A));
     261           3 :   PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_FullBasis_Sinvert));
     262           3 :   PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_FullBasis_Sinvert));
     263           3 :   PetscCall(STShellSetApply(st,Apply_FullBasis));
     264           3 :   PetscCall(STShellSetApplyTranspose(st,ApplyTranspose_FullBasis));
     265           3 :   PetscCall(EPSSetOperators(ctx->eps,ctx->A,NULL));
     266           3 :   PetscCall(EPSSetProblemType(ctx->eps,EPS_NHEP));
     267           3 :   switch (nep->which) {
     268             :     case NEP_TARGET_MAGNITUDE:   which = EPS_TARGET_MAGNITUDE; break;
     269           0 :     case NEP_TARGET_REAL:        which = EPS_TARGET_REAL; break;
     270           0 :     case NEP_TARGET_IMAGINARY:   which = EPS_TARGET_IMAGINARY; break;
     271           0 :     case NEP_WHICH_USER:         which = EPS_WHICH_USER;
     272           0 :       PetscCall(EPSSetEigenvalueComparison(ctx->eps,nep->sc->comparison,nep->sc->comparisonctx));
     273             :       break;
     274           0 :     default: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Should set a target selection in NEPSetWhichEigenpairs()");
     275             :   }
     276           3 :   PetscCall(EPSSetWhichEigenpairs(ctx->eps,which));
     277           3 :   PetscCall(RGIsTrivial(nep->rg,&istrivial));
     278           3 :   if (!istrivial) PetscCall(EPSSetRG(ctx->eps,nep->rg));
     279           3 :   PetscCall(EPSSetDimensions(ctx->eps,nep->nev,nep->ncv,nep->mpd));
     280           3 :   PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(nep->tol),nep->max_it));
     281           3 :   PetscCall(EPSSetTwoSided(ctx->eps,nep->twosided));
     282             :   /* Transfer the trackall option from pep to eps */
     283           3 :   PetscCall(NEPGetTrackAll(nep,&trackall));
     284           3 :   PetscCall(EPSSetTrackAll(ctx->eps,trackall));
     285             : 
     286             :   /* process initial vector */
     287           3 :   if (nep->nini<0) {
     288           2 :     PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*nep->nloc,deg*nep->n,&veps));
     289           2 :     PetscCall(VecGetArray(veps,&epsarray));
     290           8 :     for (i=0;i<deg;i++) {
     291           6 :       if (i<-nep->nini) {
     292           2 :         PetscCall(VecGetArray(nep->IS[i],&neparray));
     293           2 :         PetscCall(PetscArraycpy(epsarray+i*nep->nloc,neparray,nep->nloc));
     294           2 :         PetscCall(VecRestoreArray(nep->IS[i],&neparray));
     295             :       } else {
     296           4 :         if (!w) PetscCall(VecDuplicate(nep->IS[0],&w));
     297           4 :         PetscCall(VecSetRandom(w,NULL));
     298           4 :         PetscCall(VecGetArray(w,&neparray));
     299           4 :         PetscCall(PetscArraycpy(epsarray+i*nep->nloc,neparray,nep->nloc));
     300           6 :         PetscCall(VecRestoreArray(w,&neparray));
     301             :       }
     302             :     }
     303           2 :     PetscCall(VecRestoreArray(veps,&epsarray));
     304           2 :     PetscCall(EPSSetInitialSpace(ctx->eps,1,&veps));
     305           2 :     PetscCall(VecDestroy(&veps));
     306           2 :     PetscCall(VecDestroy(&w));
     307           2 :     PetscCall(SlepcBasisDestroy_Private(&nep->nini,&nep->IS));
     308             :   }
     309             : 
     310           3 :   PetscCall(EPSSetUp(ctx->eps));
     311           3 :   PetscCall(EPSGetDimensions(ctx->eps,NULL,&nep->ncv,&nep->mpd));
     312           3 :   PetscCall(EPSGetTolerances(ctx->eps,NULL,&nep->max_it));
     313           3 :   PetscCall(NEPAllocateSolution(nep,0));
     314           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     315             : }
     316             : 
     317             : /*
     318             :    NEPNLEIGSExtract_None - Extracts the first block of the basis
     319             :    and normalizes the columns.
     320             : */
     321           3 : static PetscErrorCode NEPNLEIGSExtract_None(NEP nep,EPS eps)
     322             : {
     323           3 :   PetscInt          i,k,m,d;
     324           3 :   const PetscScalar *px;
     325           3 :   PetscScalar       sigma=nep->target,*b;
     326           3 :   Mat               A;
     327           3 :   Vec               xxr,xxi=NULL,w,t,xx;
     328           3 :   PetscReal         norm;
     329           3 :   NEP_NLEIGS        *ctx=(NEP_NLEIGS*)nep->data;
     330             : 
     331           3 :   PetscFunctionBegin;
     332           3 :   d = ctx->nmat-1;
     333           3 :   PetscCall(EPSGetOperators(eps,&A,NULL));
     334           3 :   PetscCall(MatCreateVecs(A,&xxr,NULL));
     335             : #if !defined(PETSC_USE_COMPLEX)
     336             :   PetscCall(VecDuplicate(xxr,&xxi));
     337             : #endif
     338           3 :   w = nep->work[0];
     339          19 :   for (i=0;i<nep->nconv;i++) {
     340          16 :     PetscCall(EPSGetEigenvector(eps,i,xxr,xxi));
     341          16 :     PetscCall(VecGetArrayRead(xxr,&px));
     342          16 :     PetscCall(VecPlaceArray(w,px));
     343          16 :     PetscCall(BVInsertVec(nep->V,i,w));
     344          16 :     PetscCall(BVNormColumn(nep->V,i,NORM_2,&norm));
     345          16 :     PetscCall(BVScaleColumn(nep->V,i,1.0/norm));
     346          16 :     PetscCall(VecResetArray(w));
     347          16 :     PetscCall(VecRestoreArrayRead(xxr,&px));
     348             :   }
     349           3 :   if (nep->twosided) {
     350           2 :     PetscCall(PetscMalloc1(ctx->nmat,&b));
     351           2 :     PetscCall(NEPNLEIGSEvalNRTFunct(nep,d,sigma,b));
     352           2 :     m = nep->nloc;
     353           2 :     xx = ctx->w[0];
     354           2 :     w = nep->work[0]; t = nep->work[1];
     355          14 :     for (k=0;k<nep->nconv;k++) {
     356          12 :       PetscCall(EPSGetLeftEigenvector(eps,k,xxr,xxi));
     357          12 :       PetscCall(VecGetArrayRead(xxr,&px));
     358          12 :       PetscCall(VecPlaceArray(xx,px+(d-1)*m));
     359          12 :       PetscCall(VecCopy(xx,w));
     360          12 :       PetscCall(VecScale(w,PetscConj(b[d-1])));
     361          12 :       PetscCall(VecResetArray(xx));
     362          36 :       for (i=0;i<d-1;i++) {
     363          24 :         PetscCall(VecPlaceArray(xx,px+i*m));
     364          24 :         PetscCall(VecAXPY(w,PetscConj(b[i]),xx));
     365          24 :         PetscCall(VecResetArray(xx));
     366             :       }
     367          12 :       PetscCall(VecConjugate(w));
     368          12 :       PetscCall(KSPSolveTranspose(ctx->ksp[0],w,t));
     369          12 :       PetscCall(VecConjugate(t));
     370          12 :       PetscCall(BVInsertVec(nep->W,k,t));
     371          12 :       PetscCall(BVNormColumn(nep->W,k,NORM_2,&norm));
     372          12 :       PetscCall(BVScaleColumn(nep->W,k,1.0/norm));
     373          12 :       PetscCall(VecRestoreArrayRead(xxr,&px));
     374             :     }
     375           2 :     PetscCall(PetscFree(b));
     376             :   }
     377           3 :   PetscCall(VecDestroy(&xxr));
     378             : #if !defined(PETSC_USE_COMPLEX)
     379             :   PetscCall(VecDestroy(&xxi));
     380             : #endif
     381           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     382             : }
     383             : 
     384           3 : PetscErrorCode NEPSolve_NLEIGS_FullBasis(NEP nep)
     385             : {
     386           3 :   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
     387           3 :   PetscInt       i;
     388           3 :   PetscScalar    eigi=0.0;
     389             : 
     390           3 :   PetscFunctionBegin;
     391           3 :   PetscCall(EPSSolve(ctx->eps));
     392           3 :   PetscCall(EPSGetConverged(ctx->eps,&nep->nconv));
     393           3 :   PetscCall(EPSGetIterationNumber(ctx->eps,&nep->its));
     394           3 :   PetscCall(EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&nep->reason));
     395             : 
     396             :   /* recover eigenvalues */
     397          19 :   for (i=0;i<nep->nconv;i++) {
     398          16 :     PetscCall(EPSGetEigenpair(ctx->eps,i,&nep->eigr[i],&eigi,NULL,NULL));
     399             : #if !defined(PETSC_USE_COMPLEX)
     400             :     PetscCheck(eigi==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex value requires complex arithmetic");
     401             : #endif
     402             :   }
     403           3 :   PetscCall(NEPNLEIGSExtract_None(nep,ctx->eps));
     404           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     405             : }
     406             : 
     407           0 : PetscErrorCode NEPNLEIGSSetEPS_NLEIGS(NEP nep,EPS eps)
     408             : {
     409           0 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     410             : 
     411           0 :   PetscFunctionBegin;
     412           0 :   PetscCall(PetscObjectReference((PetscObject)eps));
     413           0 :   PetscCall(EPSDestroy(&ctx->eps));
     414           0 :   ctx->eps = eps;
     415           0 :   nep->state = NEP_STATE_INITIAL;
     416           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     417             : }
     418             : 
     419             : /*@
     420             :    NEPNLEIGSSetEPS - Associate an eigensolver object (EPS) to the NLEIGS solver.
     421             : 
     422             :    Collective
     423             : 
     424             :    Input Parameters:
     425             : +  nep - nonlinear eigenvalue solver
     426             : -  eps - the eigensolver object
     427             : 
     428             :    Level: advanced
     429             : 
     430             : .seealso: NEPNLEIGSGetEPS()
     431             : @*/
     432           0 : PetscErrorCode NEPNLEIGSSetEPS(NEP nep,EPS eps)
     433             : {
     434           0 :   PetscFunctionBegin;
     435           0 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     436           0 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
     437           0 :   PetscCheckSameComm(nep,1,eps,2);
     438           0 :   PetscTryMethod(nep,"NEPNLEIGSSetEPS_C",(NEP,EPS),(nep,eps));
     439           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     440             : }
     441             : 
     442           6 : static PetscErrorCode EPSMonitor_NLEIGS(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
     443             : {
     444           6 :   NEP            nep = (NEP)ctx;
     445           6 :   PetscInt       i,nv = PetscMin(nest,nep->ncv);
     446             : 
     447           6 :   PetscFunctionBegin;
     448         129 :   for (i=0;i<nv;i++) {
     449         123 :     nep->eigr[i]   = eigr[i];
     450         123 :     nep->eigi[i]   = eigi[i];
     451         123 :     nep->errest[i] = errest[i];
     452             :   }
     453           6 :   PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nv,nep->eigr,nep->eigi));
     454           6 :   PetscCall(NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest));
     455           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     456             : }
     457             : 
     458           3 : PetscErrorCode NEPNLEIGSGetEPS_NLEIGS(NEP nep,EPS *eps)
     459             : {
     460           3 :   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
     461             : 
     462           3 :   PetscFunctionBegin;
     463           3 :   if (!ctx->eps) {
     464           3 :     PetscCall(EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps));
     465           3 :     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1));
     466           3 :     PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix));
     467           3 :     PetscCall(EPSAppendOptionsPrefix(ctx->eps,"nep_nleigs_"));
     468           3 :     PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options));
     469           3 :     PetscCall(EPSMonitorSet(ctx->eps,EPSMonitor_NLEIGS,nep,NULL));
     470             :   }
     471           3 :   *eps = ctx->eps;
     472           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     473             : }
     474             : 
     475             : /*@
     476             :    NEPNLEIGSGetEPS - Retrieve the eigensolver object (EPS) associated
     477             :    to the nonlinear eigenvalue solver.
     478             : 
     479             :    Collective
     480             : 
     481             :    Input Parameter:
     482             : .  nep - nonlinear eigenvalue solver
     483             : 
     484             :    Output Parameter:
     485             : .  eps - the eigensolver object
     486             : 
     487             :    Level: advanced
     488             : 
     489             : .seealso: NEPNLEIGSSetEPS()
     490             : @*/
     491           3 : PetscErrorCode NEPNLEIGSGetEPS(NEP nep,EPS *eps)
     492             : {
     493           3 :   PetscFunctionBegin;
     494           3 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     495           3 :   PetscAssertPointer(eps,2);
     496           3 :   PetscUseMethod(nep,"NEPNLEIGSGetEPS_C",(NEP,EPS*),(nep,eps));
     497           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     498             : }