LCOV - code coverage report
Current view: top level - nep/interface - neprefine.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 425 429 99.1 %
Date: 2021-08-02 00:35:43 Functions: 6 6 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14