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             :    Newton refinement for NEP, simple version
      12             : */
      13             : 
      14             : #include <slepc/private/nepimpl.h>
      15             : #include <slepcblaslapack.h>
      16             : 
      17             : #define NREF_MAXIT 10
      18             : 
      19             : typedef struct {
      20             :   VecScatter    *scatter_id,nst;
      21             :   Mat           *A;
      22             :   Vec           nv,vg,v,w;
      23             :   FN            *fn;
      24             : } NEPSimpNRefctx;
      25             : 
      26             : typedef struct {
      27             :   Mat          M1;
      28             :   Vec          M2,M3;
      29             :   PetscScalar  M4,m3;
      30             : } NEP_REFINE_MATSHELL;
      31             : 
      32          72 : static PetscErrorCode MatMult_FS(Mat M ,Vec x,Vec y)
      33             : {
      34          72 :   NEP_REFINE_MATSHELL *ctx;
      35          72 :   PetscScalar         t;
      36             : 
      37          72 :   PetscFunctionBegin;
      38          72 :   PetscCall(MatShellGetContext(M,&ctx));
      39          72 :   PetscCall(VecDot(x,ctx->M3,&t));
      40          72 :   t *= ctx->m3/ctx->M4;
      41          72 :   PetscCall(MatMult(ctx->M1,x,y));
      42          72 :   PetscCall(VecAXPY(y,-t,ctx->M2));
      43          72 :   PetscFunctionReturn(PETSC_SUCCESS);
      44             : }
      45             : 
      46          13 : static PetscErrorCode NEPSimpleNRefSetUp(NEP nep,NEPSimpNRefctx **ctx_)
      47             : {
      48          13 :   PetscInt       i,si,j,n0,m0,nloc,*idx1,*idx2,ne;
      49          13 :   IS             is1,is2;
      50          13 :   NEPSimpNRefctx *ctx;
      51          13 :   Vec            v;
      52          13 :   PetscMPIInt    rank,size;
      53          13 :   MPI_Comm       child;
      54             : 
      55          13 :   PetscFunctionBegin;
      56          13 :   PetscCall(PetscCalloc1(1,ctx_));
      57          13 :   ctx = *ctx_;
      58          13 :   if (nep->npart==1) {
      59           5 :     ctx->A  = nep->A;
      60           5 :     ctx->fn = nep->f;
      61           5 :     nep->refinesubc = NULL;
      62           5 :     ctx->scatter_id = NULL;
      63             :   } else {
      64           8 :     PetscCall(PetscSubcommGetChild(nep->refinesubc,&child));
      65           8 :     PetscCall(PetscMalloc2(nep->nt,&ctx->A,nep->npart,&ctx->scatter_id));
      66             : 
      67             :     /* Duplicate matrices */
      68          32 :     for (i=0;i<nep->nt;i++) PetscCall(MatCreateRedundantMatrix(nep->A[i],0,child,MAT_INITIAL_MATRIX,&ctx->A[i]));
      69           8 :     PetscCall(MatCreateVecs(ctx->A[0],&ctx->v,NULL));
      70             : 
      71             :     /* Duplicate FNs */
      72           8 :     PetscCall(PetscMalloc1(nep->nt,&ctx->fn));
      73          32 :     for (i=0;i<nep->nt;i++) PetscCall(FNDuplicate(nep->f[i],child,&ctx->fn[i]));
      74             : 
      75             :     /* Create scatters for sending vectors to each subcommucator */
      76           8 :     PetscCall(BVGetColumn(nep->V,0,&v));
      77           8 :     PetscCall(VecGetOwnershipRange(v,&n0,&m0));
      78           8 :     PetscCall(BVRestoreColumn(nep->V,0,&v));
      79           8 :     PetscCall(VecGetLocalSize(ctx->v,&nloc));
      80           8 :     PetscCall(PetscMalloc2(m0-n0,&idx1,m0-n0,&idx2));
      81           8 :     PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)nep),nloc,PETSC_DECIDE,&ctx->vg));
      82          24 :     for (si=0;si<nep->npart;si++) {
      83          16 :       j = 0;
      84         616 :       for (i=n0;i<m0;i++) {
      85         600 :         idx1[j]   = i;
      86         600 :         idx2[j++] = i+nep->n*si;
      87             :       }
      88          16 :       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)nep),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
      89          16 :       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)nep),(m0-n0),idx2,PETSC_COPY_VALUES,&is2));
      90          16 :       PetscCall(BVGetColumn(nep->V,0,&v));
      91          16 :       PetscCall(VecScatterCreate(v,is1,ctx->vg,is2,&ctx->scatter_id[si]));
      92          16 :       PetscCall(BVRestoreColumn(nep->V,0,&v));
      93          16 :       PetscCall(ISDestroy(&is1));
      94          16 :       PetscCall(ISDestroy(&is2));
      95             :     }
      96           8 :     PetscCall(PetscFree2(idx1,idx2));
      97             :   }
      98          13 :   if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT) {
      99           8 :     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ctx->A[0]),&rank));
     100           8 :     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ctx->A[0]),&size));
     101           8 :     if (size>1) {
     102           4 :       if (nep->npart==1) PetscCall(BVGetColumn(nep->V,0,&v));
     103           4 :       else v = ctx->v;
     104           4 :       PetscCall(VecGetOwnershipRange(v,&n0,&m0));
     105           4 :       ne = (rank == size-1)?nep->n:0;
     106           4 :       PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->A[0]),ne,PETSC_DECIDE,&ctx->nv));
     107           4 :       PetscCall(PetscMalloc1(m0-n0,&idx1));
     108         204 :       for (i=n0;i<m0;i++) {
     109         200 :         idx1[i-n0] = i;
     110             :       }
     111           4 :       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ctx->A[0]),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
     112           4 :       PetscCall(VecScatterCreate(v,is1,ctx->nv,is1,&ctx->nst));
     113           4 :       if (nep->npart==1) PetscCall(BVRestoreColumn(nep->V,0,&v));
     114           4 :       PetscCall(PetscFree(idx1));
     115           4 :       PetscCall(ISDestroy(&is1));
     116             :     }
     117          13 :   }  PetscFunctionReturn(PETSC_SUCCESS);
     118             : }
     119             : 
     120             : /*
     121             :   Gather Eigenpair idx from subcommunicator with color sc
     122             : */
     123         113 : static PetscErrorCode NEPSimpleNRefGatherEigenpair(NEP nep,NEPSimpNRefctx *ctx,PetscInt sc,PetscInt idx,PetscInt *fail)
     124             : {
     125         113 :   PetscMPIInt    nproc,p;
     126         113 :   MPI_Comm       comm=((PetscObject)nep)->comm;
     127         113 :   Vec            v;
     128         113 :   PetscScalar    *array;
     129             : 
     130         113 :   PetscFunctionBegin;
     131         113 :   PetscCallMPI(MPI_Comm_size(comm,&nproc));
     132         113 :   PetscCall(PetscMPIIntCast((nproc/nep->npart)*(sc+1)+PetscMin(sc+1,nproc%nep->npart)-1,&p));
     133         113 :   if (nep->npart>1) {
     134             :     /* Communicate convergence successful */
     135         160 :     PetscCallMPI(MPI_Bcast(fail,1,MPIU_INT,p,comm));
     136          80 :     if (!(*fail)) {
     137             :       /* Process 0 of subcommunicator sc broadcasts the eigenvalue */
     138         160 :       PetscCallMPI(MPI_Bcast(&nep->eigr[idx],1,MPIU_SCALAR,p,comm));
     139             :       /* Gather nep->V[idx] from the subcommuniator sc */
     140          80 :       PetscCall(BVGetColumn(nep->V,idx,&v));
     141          80 :       if (nep->refinesubc->color==sc) {
     142          40 :         PetscCall(VecGetArray(ctx->v,&array));
     143          40 :         PetscCall(VecPlaceArray(ctx->vg,array));
     144             :       }
     145          80 :       PetscCall(VecScatterBegin(ctx->scatter_id[sc],ctx->vg,v,INSERT_VALUES,SCATTER_REVERSE));
     146          80 :       PetscCall(VecScatterEnd(ctx->scatter_id[sc],ctx->vg,v,INSERT_VALUES,SCATTER_REVERSE));
     147          80 :       if (nep->refinesubc->color==sc) {
     148          40 :         PetscCall(VecResetArray(ctx->vg));
     149          40 :         PetscCall(VecRestoreArray(ctx->v,&array));
     150             :       }
     151          80 :       PetscCall(BVRestoreColumn(nep->V,idx,&v));
     152             :     }
     153             :   } else {
     154          43 :     if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT && !(*fail)) PetscCallMPI(MPI_Bcast(&nep->eigr[idx],1,MPIU_SCALAR,p,comm));
     155             :   }
     156         113 :   PetscFunctionReturn(PETSC_SUCCESS);
     157             : }
     158             : 
     159          83 : static PetscErrorCode NEPSimpleNRefScatterEigenvector(NEP nep,NEPSimpNRefctx *ctx,PetscInt sc,PetscInt idx)
     160             : {
     161          83 :   Vec            v;
     162          83 :   PetscScalar    *array;
     163             : 
     164          83 :   PetscFunctionBegin;
     165          83 :   if (nep->npart>1) {
     166          56 :     PetscCall(BVGetColumn(nep->V,idx,&v));
     167          56 :     if (nep->refinesubc->color==sc) {
     168          28 :       PetscCall(VecGetArray(ctx->v,&array));
     169          28 :       PetscCall(VecPlaceArray(ctx->vg,array));
     170             :     }
     171          56 :     PetscCall(VecScatterBegin(ctx->scatter_id[sc],v,ctx->vg,INSERT_VALUES,SCATTER_FORWARD));
     172          56 :     PetscCall(VecScatterEnd(ctx->scatter_id[sc],v,ctx->vg,INSERT_VALUES,SCATTER_FORWARD));
     173          56 :     if (nep->refinesubc->color==sc) {
     174          28 :       PetscCall(VecResetArray(ctx->vg));
     175          28 :       PetscCall(VecRestoreArray(ctx->v,&array));
     176             :     }
     177          56 :     PetscCall(BVRestoreColumn(nep->V,idx,&v));
     178             :   }
     179          83 :   PetscFunctionReturn(PETSC_SUCCESS);
     180             : }
     181             : 
     182          73 : static PetscErrorCode NEPSimpleNRefSetUpSystem(NEP nep,NEPSimpNRefctx *ctx,Mat *A,PetscInt idx,Mat *Mt,Mat *T,Mat *P,PetscBool ini,Vec t,Vec v)
     183             : {
     184          73 :   PetscInt            i,st,ml,m0,n0,m1,mg;
     185          73 :   PetscInt            ncols,*cols2=NULL,nt=nep->nt;
     186          73 :   PetscScalar         zero=0.0,*coeffs,*coeffs2;
     187          73 :   PetscMPIInt         rank,size;
     188          73 :   MPI_Comm            comm;
     189          73 :   const PetscInt      *cols;
     190          73 :   const PetscScalar   *vals,*array;
     191          73 :   NEP_REFINE_MATSHELL *fctx;
     192          73 :   Vec                 w=ctx->w;
     193          73 :   Mat                 M;
     194             : 
     195          73 :   PetscFunctionBegin;
     196          73 :   PetscCall(PetscMalloc2(nt,&coeffs,nt,&coeffs2));
     197          73 :   switch (nep->scheme) {
     198          14 :   case NEP_REFINE_SCHEME_SCHUR:
     199          14 :     if (ini) {
     200           2 :       PetscCall(PetscCalloc1(1,&fctx));
     201           2 :       PetscCall(MatGetSize(A[0],&m0,&n0));
     202           2 :       PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A[0]),PETSC_DECIDE,PETSC_DECIDE,m0,n0,fctx,T));
     203           2 :       PetscCall(MatShellSetOperation(*T,MATOP_MULT,(void(*)(void))MatMult_FS));
     204          12 :     } else PetscCall(MatShellGetContext(*T,&fctx));
     205          14 :     M=fctx->M1;
     206          14 :     break;
     207          18 :   case NEP_REFINE_SCHEME_MBE:
     208          18 :     M=*T;
     209          18 :     break;
     210          41 :   case NEP_REFINE_SCHEME_EXPLICIT:
     211          41 :     M=*Mt;
     212          41 :     break;
     213             :   }
     214          73 :   if (ini) PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&M));
     215          60 :   else PetscCall(MatCopy(A[0],M,DIFFERENT_NONZERO_PATTERN));
     216         292 :   for (i=0;i<nt;i++) PetscCall(FNEvaluateFunction(ctx->fn[i],nep->eigr[idx],coeffs+i));
     217          73 :   if (coeffs[0]!=1.0) PetscCall(MatScale(M,coeffs[0]));
     218         219 :   for (i=1;i<nt;i++) PetscCall(MatAXPY(M,coeffs[i],A[i],(ini)?nep->mstr:SUBSET_NONZERO_PATTERN));
     219         292 :   for (i=0;i<nt;i++) PetscCall(FNEvaluateDerivative(ctx->fn[i],nep->eigr[idx],coeffs2+i));
     220          73 :   st = 0;
     221         145 :   for (i=0;i<nt && PetscAbsScalar(coeffs2[i])==0.0;i++) st++;
     222          73 :   PetscCall(MatMult(A[st],v,w));
     223          73 :   if (coeffs2[st]!=1.0) PetscCall(VecScale(w,coeffs2[st]));
     224         147 :   for (i=st+1;i<nt;i++) {
     225          74 :     PetscCall(MatMult(A[i],v,t));
     226          74 :     PetscCall(VecAXPY(w,coeffs2[i],t));
     227             :   }
     228             : 
     229          73 :   switch (nep->scheme) {
     230          41 :   case NEP_REFINE_SCHEME_EXPLICIT:
     231          41 :     comm = PetscObjectComm((PetscObject)A[0]);
     232          41 :     PetscCallMPI(MPI_Comm_rank(comm,&rank));
     233          41 :     PetscCallMPI(MPI_Comm_size(comm,&size));
     234          41 :     PetscCall(MatGetSize(M,&mg,NULL));
     235          41 :     PetscCall(MatGetOwnershipRange(M,&m0,&m1));
     236          41 :     if (ini) {
     237           8 :       PetscCall(MatCreate(comm,T));
     238           8 :       PetscCall(MatGetLocalSize(M,&ml,NULL));
     239           8 :       if (rank==size-1) ml++;
     240           8 :       PetscCall(MatSetSizes(*T,ml,ml,mg+1,mg+1));
     241           8 :       PetscCall(MatSetFromOptions(*T));
     242           8 :       *Mt = M;
     243           8 :       *P  = *T;
     244             :     }
     245             : 
     246             :     /* Set values */
     247          41 :     PetscCall(VecGetArrayRead(w,&array));
     248        2961 :     for (i=m0;i<m1;i++) {
     249        2920 :       PetscCall(MatGetRow(M,i,&ncols,&cols,&vals));
     250        2920 :       PetscCall(MatSetValues(*T,1,&i,ncols,cols,vals,INSERT_VALUES));
     251        2920 :       PetscCall(MatRestoreRow(M,i,&ncols,&cols,&vals));
     252        2920 :       PetscCall(MatSetValues(*T,1,&i,1,&mg,array+i-m0,INSERT_VALUES));
     253             :     }
     254          41 :     PetscCall(VecRestoreArrayRead(w,&array));
     255          41 :     PetscCall(VecConjugate(v));
     256          41 :     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A[0]),&size));
     257          41 :     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A[0]),&rank));
     258          41 :     if (size>1) {
     259          22 :       if (rank==size-1) {
     260          11 :         PetscCall(PetscMalloc1(nep->n,&cols2));
     261        1111 :         for (i=0;i<nep->n;i++) cols2[i]=i;
     262             :       }
     263          22 :       PetscCall(VecScatterBegin(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
     264          22 :       PetscCall(VecScatterEnd(ctx->nst,v,ctx->nv,INSERT_VALUES,SCATTER_FORWARD));
     265          22 :       PetscCall(VecGetArrayRead(ctx->nv,&array));
     266          22 :       if (rank==size-1) {
     267          11 :         PetscCall(MatSetValues(*T,1,&mg,nep->n,cols2,array,INSERT_VALUES));
     268          11 :         PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
     269             :       }
     270          22 :       PetscCall(VecRestoreArrayRead(ctx->nv,&array));
     271             :     } else {
     272          19 :       PetscCall(PetscMalloc1(m1-m0,&cols2));
     273        1839 :       for (i=0;i<m1-m0;i++) cols2[i]=m0+i;
     274          19 :       PetscCall(VecGetArrayRead(v,&array));
     275          19 :       PetscCall(MatSetValues(*T,1,&mg,m1-m0,cols2,array,INSERT_VALUES));
     276          19 :       PetscCall(MatSetValues(*T,1,&mg,1,&mg,&zero,INSERT_VALUES));
     277          19 :       PetscCall(VecRestoreArrayRead(v,&array));
     278             :     }
     279          41 :     PetscCall(VecConjugate(v));
     280          41 :     PetscCall(MatAssemblyBegin(*T,MAT_FINAL_ASSEMBLY));
     281          41 :     PetscCall(MatAssemblyEnd(*T,MAT_FINAL_ASSEMBLY));
     282          41 :     PetscCall(PetscFree(cols2));
     283             :     break;
     284          14 :   case NEP_REFINE_SCHEME_SCHUR:
     285          14 :     fctx->M2 = ctx->w;
     286          14 :     fctx->M3 = v;
     287          14 :     fctx->m3 = 1.0+PetscConj(nep->eigr[idx])*nep->eigr[idx];
     288          14 :     fctx->M4 = PetscConj(nep->eigr[idx]);
     289          14 :     fctx->M1 = M;
     290          14 :     if (ini) PetscCall(MatDuplicate(M,MAT_COPY_VALUES,P));
     291          12 :     else PetscCall(MatCopy(M,*P,SAME_NONZERO_PATTERN));
     292          14 :     if (fctx->M4!=0.0) {
     293          14 :       PetscCall(VecConjugate(v));
     294          14 :       PetscCall(VecPointwiseMult(t,v,w));
     295          14 :       PetscCall(VecConjugate(v));
     296          14 :       PetscCall(VecScale(t,-fctx->m3/fctx->M4));
     297          14 :       PetscCall(MatDiagonalSet(*P,t,ADD_VALUES));
     298             :     }
     299             :     break;
     300          18 :   case NEP_REFINE_SCHEME_MBE:
     301          18 :     *T = M;
     302          18 :     *P = M;
     303          18 :     break;
     304             :   }
     305          73 :   PetscCall(PetscFree2(coeffs,coeffs2));
     306          73 :   PetscFunctionReturn(PETSC_SUCCESS);
     307             : }
     308             : 
     309          13 : PetscErrorCode NEPNewtonRefinementSimple(NEP nep,PetscInt *maxits,PetscReal tol,PetscInt k)
     310             : {
     311          13 :   PetscInt            i,n,its,idx=0,*idx_sc,*its_sc,color,*fail_sc;
     312          13 :   PetscMPIInt         rank,size;
     313          13 :   Mat                 Mt=NULL,T=NULL,P=NULL;
     314          13 :   MPI_Comm            comm;
     315          13 :   Vec                 r,v,dv,rr=NULL,dvv=NULL,t[2];
     316          13 :   const PetscScalar   *array;
     317          13 :   PetscScalar         *array2,deig=0.0,tt[2],ttt;
     318          13 :   PetscReal           norm,error;
     319          13 :   PetscBool           ini=PETSC_TRUE,sc_pend,solved=PETSC_FALSE;
     320          13 :   NEPSimpNRefctx      *ctx;
     321          13 :   NEP_REFINE_MATSHELL *fctx=NULL;
     322          13 :   KSPConvergedReason  reason;
     323             : 
     324          13 :   PetscFunctionBegin;
     325          13 :   PetscCall(PetscLogEventBegin(NEP_Refine,nep,0,0,0));
     326          13 :   PetscCall(NEPSimpleNRefSetUp(nep,&ctx));
     327          13 :   its = (maxits)?*maxits:NREF_MAXIT;
     328          13 :   if (!nep->refineksp) PetscCall(NEPRefineGetKSP(nep,&nep->refineksp));
     329          13 :   if (nep->npart==1) PetscCall(BVGetColumn(nep->V,0,&v));
     330           8 :   else v = ctx->v;
     331          13 :   PetscCall(VecDuplicate(v,&ctx->w));
     332          13 :   PetscCall(VecDuplicate(v,&r));
     333          13 :   PetscCall(VecDuplicate(v,&dv));
     334          13 :   PetscCall(VecDuplicate(v,&t[0]));
     335          13 :   PetscCall(VecDuplicate(v,&t[1]));
     336          13 :   if (nep->npart==1) {
     337           5 :     PetscCall(BVRestoreColumn(nep->V,0,&v));
     338           5 :     PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
     339           8 :   } else PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
     340          13 :   PetscCallMPI(MPI_Comm_size(comm,&size));
     341          13 :   PetscCallMPI(MPI_Comm_rank(comm,&rank));
     342          13 :   PetscCall(VecGetLocalSize(r,&n));
     343          13 :   PetscCall(PetscMalloc3(nep->npart,&idx_sc,nep->npart,&its_sc,nep->npart,&fail_sc));
     344          34 :   for (i=0;i<nep->npart;i++) fail_sc[i] = 0;
     345          34 :   for (i=0;i<nep->npart;i++) its_sc[i] = 0;
     346          13 :   color = (nep->npart==1)?0:nep->refinesubc->color;
     347             : 
     348             :   /* Loop performing iterative refinements */
     349          13 :   while (!solved) {
     350         232 :     for (i=0;i<nep->npart;i++) {
     351         142 :       sc_pend = PETSC_TRUE;
     352         142 :       if (its_sc[i]==0) {
     353          29 :         idx_sc[i] = idx++;
     354          29 :         if (idx_sc[i]>=k) {
     355             :           sc_pend = PETSC_FALSE;
     356          21 :         } else PetscCall(NEPSimpleNRefScatterEigenvector(nep,ctx,i,idx_sc[i]));
     357             :       }  else { /* Gather Eigenpair from subcommunicator i */
     358         113 :         PetscCall(NEPSimpleNRefGatherEigenpair(nep,ctx,i,idx_sc[i],&fail_sc[i]));
     359             :       }
     360         317 :       while (sc_pend) {
     361         196 :         if (!fail_sc[i]) PetscCall(NEPComputeError(nep,idx_sc[i],NEP_ERROR_RELATIVE,&error));
     362         196 :         if (error<=tol || its_sc[i]>=its || fail_sc[i]) {
     363          83 :           idx_sc[i] = idx++;
     364          83 :           its_sc[i] = 0;
     365          83 :           fail_sc[i] = 0;
     366          83 :           if (idx_sc[i]<k) PetscCall(NEPSimpleNRefScatterEigenvector(nep,ctx,i,idx_sc[i]));
     367             :         } else {
     368         113 :           sc_pend = PETSC_FALSE;
     369         113 :           its_sc[i]++;
     370             :         }
     371         196 :         if (idx_sc[i]>=k) sc_pend = PETSC_FALSE;
     372             :       }
     373             :     }
     374             :     solved = PETSC_TRUE;
     375         188 :     for (i=0;i<nep->npart&&solved;i++) solved = PetscNot(idx_sc[i]<k);
     376          90 :     if (idx_sc[color]<k) {
     377             : #if !defined(PETSC_USE_COMPLEX)
     378          73 :       PetscCheck(nep->eigi[idx_sc[color]]==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Simple Refinement not implemented in real scalar for complex eigenvalues");
     379             : #endif
     380          73 :       if (nep->npart==1) PetscCall(BVGetColumn(nep->V,idx_sc[color],&v));
     381          40 :       else v = ctx->v;
     382          73 :       PetscCall(NEPSimpleNRefSetUpSystem(nep,ctx,ctx->A,idx_sc[color],&Mt,&T,&P,ini,t[0],v));
     383          73 :       PetscCall(NEP_KSPSetOperators(nep->refineksp,T,P));
     384          73 :       if (ini) {
     385          13 :         PetscCall(KSPSetFromOptions(nep->refineksp));
     386          13 :         if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT) {
     387           8 :           PetscCall(MatCreateVecs(T,&dvv,NULL));
     388           8 :           PetscCall(VecDuplicate(dvv,&rr));
     389             :         }
     390             :         ini = PETSC_FALSE;
     391             :       }
     392          73 :       switch (nep->scheme) {
     393          41 :       case NEP_REFINE_SCHEME_EXPLICIT:
     394          41 :         PetscCall(MatMult(Mt,v,r));
     395          41 :         PetscCall(VecGetArrayRead(r,&array));
     396          41 :         if (rank==size-1) {
     397          30 :           PetscCall(VecGetArray(rr,&array2));
     398          30 :           PetscCall(PetscArraycpy(array2,array,n));
     399          30 :           array2[n] = 0.0;
     400          30 :           PetscCall(VecRestoreArray(rr,&array2));
     401          11 :         } else PetscCall(VecPlaceArray(rr,array));
     402          41 :         PetscCall(KSPSolve(nep->refineksp,rr,dvv));
     403          41 :         PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
     404          41 :         if (reason>0) {
     405          41 :           if (rank != size-1) PetscCall(VecResetArray(rr));
     406          41 :           PetscCall(VecRestoreArrayRead(r,&array));
     407          41 :           PetscCall(VecGetArrayRead(dvv,&array));
     408          41 :           PetscCall(VecPlaceArray(dv,array));
     409          41 :           PetscCall(VecAXPY(v,-1.0,dv));
     410          41 :           PetscCall(VecNorm(v,NORM_2,&norm));
     411          41 :           PetscCall(VecScale(v,1.0/norm));
     412          41 :           PetscCall(VecResetArray(dv));
     413          41 :           if (rank==size-1) nep->eigr[idx_sc[color]] -= array[n];
     414          41 :           PetscCall(VecRestoreArrayRead(dvv,&array));
     415           0 :         } else fail_sc[color] = 1;
     416             :         break;
     417          18 :       case NEP_REFINE_SCHEME_MBE:
     418          18 :         PetscCall(MatMult(T,v,r));
     419             :         /* Mixed block elimination */
     420          18 :         PetscCall(VecConjugate(v));
     421          18 :         PetscCall(KSPSolveTranspose(nep->refineksp,v,t[0]));
     422          18 :         PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
     423          18 :         if (reason>0) {
     424          18 :           PetscCall(VecConjugate(t[0]));
     425          18 :           PetscCall(VecDot(ctx->w,t[0],&tt[0]));
     426          18 :           PetscCall(KSPSolve(nep->refineksp,ctx->w,t[1]));
     427          18 :           PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
     428          18 :           if (reason>0) {
     429          18 :             PetscCall(VecDot(t[1],v,&tt[1]));
     430          18 :             PetscCall(VecDot(r,t[0],&ttt));
     431          18 :             tt[0] = ttt/tt[0];
     432          18 :             PetscCall(VecAXPY(r,-tt[0],ctx->w));
     433          18 :             PetscCall(KSPSolve(nep->refineksp,r,dv));
     434          18 :             PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
     435          18 :             if (reason>0) {
     436          18 :               PetscCall(VecDot(dv,v,&ttt));
     437          18 :               tt[1] = ttt/tt[1];
     438          18 :               PetscCall(VecAXPY(dv,-tt[1],t[1]));
     439          18 :               deig = tt[0]+tt[1];
     440             :             }
     441             :           }
     442          18 :           PetscCall(VecConjugate(v));
     443          18 :           PetscCall(VecAXPY(v,-1.0,dv));
     444          18 :           PetscCall(VecNorm(v,NORM_2,&norm));
     445          18 :           PetscCall(VecScale(v,1.0/norm));
     446          18 :           nep->eigr[idx_sc[color]] -= deig;
     447          18 :           fail_sc[color] = 0;
     448             :         } else {
     449           0 :           PetscCall(VecConjugate(v));
     450           0 :           fail_sc[color] = 1;
     451             :         }
     452             :         break;
     453          14 :       case NEP_REFINE_SCHEME_SCHUR:
     454          14 :         fail_sc[color] = 1;
     455          14 :         PetscCall(MatShellGetContext(T,&fctx));
     456          14 :         if (fctx->M4!=0.0) {
     457          14 :           PetscCall(MatMult(fctx->M1,v,r));
     458          14 :           PetscCall(KSPSolve(nep->refineksp,r,dv));
     459          14 :           PetscCall(KSPGetConvergedReason(nep->refineksp,&reason));
     460          14 :           if (reason>0) {
     461          14 :             PetscCall(VecDot(dv,v,&deig));
     462          14 :             deig *= -fctx->m3/fctx->M4;
     463          14 :             PetscCall(VecAXPY(v,-1.0,dv));
     464          14 :             PetscCall(VecNorm(v,NORM_2,&norm));
     465          14 :             PetscCall(VecScale(v,1.0/norm));
     466          14 :             nep->eigr[idx_sc[color]] -= deig;
     467          14 :             fail_sc[color] = 0;
     468             :           }
     469             :         }
     470             :         break;
     471             :       }
     472         176 :       if (nep->npart==1) PetscCall(BVRestoreColumn(nep->V,idx_sc[color],&v));
     473             :     }
     474             :   }
     475          13 :   PetscCall(VecDestroy(&t[0]));
     476          13 :   PetscCall(VecDestroy(&t[1]));
     477          13 :   PetscCall(VecDestroy(&dv));
     478          13 :   PetscCall(VecDestroy(&ctx->w));
     479          13 :   PetscCall(VecDestroy(&r));
     480          13 :   PetscCall(PetscFree3(idx_sc,its_sc,fail_sc));
     481          13 :   PetscCall(VecScatterDestroy(&ctx->nst));
     482          13 :   if (nep->npart>1) {
     483           8 :     PetscCall(VecDestroy(&ctx->vg));
     484           8 :     PetscCall(VecDestroy(&ctx->v));
     485          32 :     for (i=0;i<nep->nt;i++) PetscCall(MatDestroy(&ctx->A[i]));
     486          24 :     for (i=0;i<nep->npart;i++) PetscCall(VecScatterDestroy(&ctx->scatter_id[i]));
     487           8 :     PetscCall(PetscFree2(ctx->A,ctx->scatter_id));
     488             :   }
     489          13 :   if (fctx && nep->scheme==NEP_REFINE_SCHEME_SCHUR) {
     490           2 :     PetscCall(MatDestroy(&P));
     491           2 :     PetscCall(MatDestroy(&fctx->M1));
     492           2 :     PetscCall(PetscFree(fctx));
     493             :   }
     494          13 :   if (nep->scheme==NEP_REFINE_SCHEME_EXPLICIT) {
     495           8 :     PetscCall(MatDestroy(&Mt));
     496           8 :     PetscCall(VecDestroy(&dvv));
     497           8 :     PetscCall(VecDestroy(&rr));
     498           8 :     PetscCall(VecDestroy(&ctx->nv));
     499           8 :     if (nep->npart>1) {
     500          24 :       for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&ctx->fn[i]));
     501           6 :       PetscCall(PetscFree(ctx->fn));
     502             :     }
     503             :   }
     504          13 :   if (nep->scheme==NEP_REFINE_SCHEME_MBE) {
     505           3 :     if (nep->npart>1) {
     506           8 :       for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&ctx->fn[i]));
     507           2 :       PetscCall(PetscFree(ctx->fn));
     508             :     }
     509             :   }
     510          13 :   PetscCall(MatDestroy(&T));
     511          13 :   PetscCall(PetscFree(ctx));
     512          13 :   PetscCall(PetscLogEventEnd(NEP_Refine,nep,0,0,0));
     513          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     514             : }