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