LCOV - code coverage report
Current view: top level - nep/impls/slp - slp-twosided.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 336 357 94.1 %
Date: 2024-03-29 00:55:19 Functions: 17 18 94.4 %
Legend: Lines: hit not hit

          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             :    Two-sided variant of the NEPSLP solver.
      12             : */
      13             : 
      14             : #include <slepc/private/nepimpl.h>         /*I "slepcnep.h" I*/
      15             : #include "slp.h"
      16             : 
      17             : typedef struct _n_nep_def_ctx *NEP_NEDEF_CTX;
      18             : 
      19             : struct _n_nep_def_ctx {
      20             :   PetscInt    n;
      21             :   PetscBool   ref;
      22             :   PetscScalar *eig;
      23             :   BV          V,W;
      24             : };
      25             : 
      26             : typedef struct {   /* context for two-sided solver */
      27             :   Mat         Ft;
      28             :   Mat         Jt;
      29             :   Vec         w;
      30             : } NEP_SLPTS_MATSHELL;
      31             : 
      32             : typedef struct {   /* context for non-equivalence deflation */
      33             :   NEP_NEDEF_CTX defctx;
      34             :   Mat           F;
      35             :   Mat           P;
      36             :   Mat           J;
      37             :   KSP           ksp;
      38             :   PetscBool     isJ;
      39             :   PetscScalar   lambda;
      40             :   Vec           w[2];
      41             : } NEP_NEDEF_MATSHELL;
      42             : 
      43         468 : static PetscErrorCode MatMult_SLPTS_Right(Mat M,Vec x,Vec y)
      44             : {
      45         468 :   NEP_SLPTS_MATSHELL *ctx;
      46             : 
      47         468 :   PetscFunctionBegin;
      48         468 :   PetscCall(MatShellGetContext(M,&ctx));
      49         468 :   PetscCall(MatMult(ctx->Jt,x,ctx->w));
      50         468 :   PetscCall(MatSolve(ctx->Ft,ctx->w,y));
      51         468 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54         496 : static PetscErrorCode MatMult_SLPTS_Left(Mat M,Vec x,Vec y)
      55             : {
      56         496 :   NEP_SLPTS_MATSHELL *ctx;
      57             : 
      58         496 :   PetscFunctionBegin;
      59         496 :   PetscCall(MatShellGetContext(M,&ctx));
      60         496 :   PetscCall(MatMultTranspose(ctx->Jt,x,ctx->w));
      61         496 :   PetscCall(MatSolveTranspose(ctx->Ft,ctx->w,y));
      62         496 :   PetscFunctionReturn(PETSC_SUCCESS);
      63             : }
      64             : 
      65          10 : static PetscErrorCode MatDestroy_SLPTS(Mat M)
      66             : {
      67          10 :   NEP_SLPTS_MATSHELL *ctx;
      68             : 
      69          10 :   PetscFunctionBegin;
      70          10 :   PetscCall(MatShellGetContext(M,&ctx));
      71          10 :   PetscCall(VecDestroy(&ctx->w));
      72          10 :   PetscCall(PetscFree(ctx));
      73          10 :   PetscFunctionReturn(PETSC_SUCCESS);
      74             : }
      75             : 
      76             : #if defined(PETSC_HAVE_CUDA)
      77             : static PetscErrorCode MatCreateVecs_SLPTS(Mat M,Vec *left,Vec *right)
      78             : {
      79             :   NEP_SLPTS_MATSHELL *ctx;
      80             : 
      81             :   PetscFunctionBegin;
      82             :   PetscCall(MatShellGetContext(M,&ctx));
      83             :   if (right) PetscCall(VecDuplicate(ctx->w,right));
      84             :   if (left) PetscCall(VecDuplicate(ctx->w,left));
      85             :   PetscFunctionReturn(PETSC_SUCCESS);
      86             : }
      87             : #endif
      88             : 
      89          10 : static PetscErrorCode NEPSLPSetUpEPSMat(NEP nep,Mat F,Mat J,PetscBool left,Mat *M)
      90             : {
      91          10 :   Mat                Mshell;
      92          10 :   PetscInt           nloc,mloc;
      93          10 :   NEP_SLPTS_MATSHELL *shellctx;
      94             : 
      95          10 :   PetscFunctionBegin;
      96             :   /* Create mat shell */
      97          10 :   PetscCall(PetscNew(&shellctx));
      98          10 :   shellctx->Ft = F;
      99          10 :   shellctx->Jt = J;
     100          10 :   PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
     101          10 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell));
     102          10 :   if (left) PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLPTS_Left));
     103           5 :   else PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLPTS_Right));
     104          10 :   PetscCall(MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_SLPTS));
     105             : #if defined(PETSC_HAVE_CUDA)
     106             :   PetscCall(MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_SLPTS));
     107             : #endif
     108          10 :   *M = Mshell;
     109          10 :   PetscCall(MatCreateVecs(nep->function,&shellctx->w,NULL));
     110          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     111             : }
     112             : 
     113             : /* Functions for deflation */
     114           5 : static PetscErrorCode NEPDeflationNEDestroy(NEP_NEDEF_CTX defctx)
     115             : {
     116           5 :   PetscFunctionBegin;
     117           5 :   if (!defctx) PetscFunctionReturn(PETSC_SUCCESS);
     118           5 :   PetscCall(PetscFree(defctx->eig));
     119           5 :   PetscCall(PetscFree(defctx));
     120           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     121             : }
     122             : 
     123           5 : static PetscErrorCode NEPDeflationNECreate(NEP nep,BV V,BV W,PetscInt sz,NEP_NEDEF_CTX *defctx)
     124             : {
     125           5 :   NEP_NEDEF_CTX  op;
     126             : 
     127           5 :   PetscFunctionBegin;
     128           5 :   PetscCall(PetscNew(&op));
     129           5 :   *defctx = op;
     130           5 :   op->n   = 0;
     131           5 :   op->ref = PETSC_FALSE;
     132           5 :   PetscCall(PetscCalloc1(sz,&op->eig));
     133           5 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     134           5 :   PetscCall(PetscObjectStateIncrease((PetscObject)W));
     135           5 :   op->V = V;
     136           5 :   op->W = W;
     137           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     138             : }
     139             : 
     140          57 : static PetscErrorCode NEPDeflationNEComputeFunction(NEP nep,Mat M,PetscScalar lambda)
     141             : {
     142          57 :   NEP_NEDEF_MATSHELL *matctx;
     143             : 
     144          57 :   PetscFunctionBegin;
     145          57 :   PetscCall(MatShellGetContext(M,&matctx));
     146          57 :   if (lambda==matctx->lambda) PetscFunctionReturn(PETSC_SUCCESS);
     147          41 :   PetscCall(NEPComputeFunction(nep,lambda,matctx->F,matctx->P));
     148          41 :   if (matctx->isJ) PetscCall(NEPComputeJacobian(nep,lambda,matctx->J));
     149          41 :   if (matctx->ksp) PetscCall(NEP_KSPSetOperators(matctx->ksp,matctx->F,matctx->P));
     150          41 :   matctx->lambda = lambda;
     151          41 :   PetscFunctionReturn(PETSC_SUCCESS);
     152             : }
     153             : 
     154         497 : static PetscErrorCode MatMult_NEPDeflationNE(Mat M,Vec x,Vec r)
     155             : {
     156         497 :   Vec                t,tt;
     157         497 :   PetscScalar        *h,*alpha,lambda,*eig;
     158         497 :   PetscInt           i,k;
     159         497 :   NEP_NEDEF_MATSHELL *matctx;
     160             : 
     161         497 :   PetscFunctionBegin;
     162         497 :   PetscCall(MatShellGetContext(M,&matctx));
     163         497 :   if (matctx->defctx->n && !matctx->defctx->ref) {
     164          89 :     k = matctx->defctx->n;
     165          89 :     lambda = matctx->lambda;
     166          89 :     eig = matctx->defctx->eig;
     167          89 :     t = matctx->w[0];
     168          89 :     PetscCall(VecCopy(x,t));
     169          89 :     PetscCall(PetscMalloc2(k,&h,k,&alpha));
     170         196 :     for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
     171          89 :     PetscCall(BVDotVec(matctx->defctx->V,t,h));
     172         196 :     for (i=0;i<k;i++) h[i] *= alpha[i];
     173          89 :     PetscCall(BVMultVec(matctx->defctx->W,-1.0,1.0,t,h));
     174          89 :     PetscCall(MatMult(matctx->isJ?matctx->J:matctx->F,t,r));
     175          89 :     if (matctx->isJ) {
     176         176 :       for (i=0;i<k;i++) h[i] *= (1.0/((lambda-eig[i])*(lambda-eig[i])))/alpha[i];
     177          80 :       tt = matctx->w[1];
     178          80 :       PetscCall(BVMultVec(matctx->defctx->W,-1.0,0.0,tt,h));
     179          80 :       PetscCall(MatMult(matctx->F,tt,t));
     180          80 :       PetscCall(VecAXPY(r,1.0,t));
     181             :     }
     182          89 :     PetscCall(PetscFree2(h,alpha));
     183         408 :   } else PetscCall(MatMult(matctx->isJ?matctx->J:matctx->F,x,r));
     184         497 :   PetscFunctionReturn(PETSC_SUCCESS);
     185             : }
     186             : 
     187         525 : static PetscErrorCode MatMultTranspose_NEPDeflationNE(Mat M,Vec x,Vec r)
     188             : {
     189         525 :   Vec                t,tt;
     190         525 :   PetscScalar        *h,*alphaC,lambda,*eig;
     191         525 :   PetscInt           i,k;
     192         525 :   NEP_NEDEF_MATSHELL *matctx;
     193             : 
     194         525 :   PetscFunctionBegin;
     195         525 :   PetscCall(MatShellGetContext(M,&matctx));
     196         525 :   t    = matctx->w[0];
     197         525 :   PetscCall(VecCopy(x,t));
     198         525 :   if (matctx->defctx->n && !matctx->defctx->ref) {
     199          89 :     PetscCall(VecConjugate(t));
     200          89 :     k = matctx->defctx->n;
     201          89 :     lambda = matctx->lambda;
     202          89 :     eig = matctx->defctx->eig;
     203          89 :     PetscCall(PetscMalloc2(k,&h,k,&alphaC));
     204         196 :     for (i=0;i<k;i++) alphaC[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
     205          89 :     PetscCall(BVDotVec(matctx->defctx->W,t,h));
     206         196 :     for (i=0;i<k;i++) h[i] *= alphaC[i];
     207          89 :     PetscCall(BVMultVec(matctx->defctx->V,-1.0,1.0,t,h));
     208          89 :     PetscCall(VecConjugate(t));
     209          89 :     PetscCall(MatMultTranspose(matctx->isJ?matctx->J:matctx->F,t,r));
     210          89 :     if (matctx->isJ) {
     211         176 :       for (i=0;i<k;i++) h[i] *= PetscConj(1.0/((lambda-eig[i])*(lambda-eig[i])))/alphaC[i];
     212          80 :       tt = matctx->w[1];
     213          80 :       PetscCall(BVMultVec(matctx->defctx->V,-1.0,0.0,tt,h));
     214          80 :       PetscCall(VecConjugate(tt));
     215          80 :       PetscCall(MatMultTranspose(matctx->F,tt,t));
     216          80 :       PetscCall(VecAXPY(r,1.0,t));
     217             :     }
     218          89 :     PetscCall(PetscFree2(h,alphaC));
     219         436 :   } else PetscCall(MatMultTranspose(matctx->isJ?matctx->J:matctx->F,t,r));
     220         525 :   PetscFunctionReturn(PETSC_SUCCESS);
     221             : }
     222             : 
     223         468 : static PetscErrorCode MatSolve_NEPDeflationNE(Mat M,Vec b,Vec x)
     224             : {
     225         468 :   PetscScalar        *h,*alpha,lambda,*eig;
     226         468 :   PetscInt           i,k;
     227         468 :   NEP_NEDEF_MATSHELL *matctx;
     228             : 
     229         468 :   PetscFunctionBegin;
     230         468 :   PetscCall(MatShellGetContext(M,&matctx));
     231         468 :   if (!matctx->ksp) {
     232           0 :     PetscCall(VecCopy(b,x));
     233           0 :     PetscFunctionReturn(PETSC_SUCCESS);
     234             :   }
     235         468 :   PetscCall(KSPSolve(matctx->ksp,b,x));
     236         468 :   if (matctx->defctx->n && !matctx->defctx->ref) {
     237          80 :     k = matctx->defctx->n;
     238          80 :     lambda = matctx->lambda;
     239          80 :     eig = matctx->defctx->eig;
     240          80 :     PetscCall(PetscMalloc2(k,&h,k,&alpha));
     241          80 :     PetscCall(BVDotVec(matctx->defctx->V,x,h));
     242         176 :     for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
     243         176 :     for (i=0;i<k;i++) h[i] *= alpha[i]/(1.0-alpha[i]);
     244          80 :     PetscCall(BVMultVec(matctx->defctx->W,1.0,1.0,x,h));
     245          80 :     PetscCall(PetscFree2(h,alpha));
     246             :   }
     247         468 :   PetscFunctionReturn(PETSC_SUCCESS);
     248             : }
     249             : 
     250         496 : static PetscErrorCode MatSolveTranspose_NEPDeflationNE(Mat M,Vec b,Vec x)
     251             : {
     252         496 :   PetscScalar        *h,*alphaC,lambda,*eig;
     253         496 :   PetscInt           i,k;
     254         496 :   NEP_NEDEF_MATSHELL *matctx;
     255             : 
     256         496 :   PetscFunctionBegin;
     257         496 :   PetscCall(MatShellGetContext(M,&matctx));
     258         496 :   if (!matctx->ksp) {
     259           0 :     PetscCall(VecCopy(b,x));
     260           0 :     PetscFunctionReturn(PETSC_SUCCESS);
     261             :   }
     262         496 :   PetscCall(KSPSolveTranspose(matctx->ksp,b,x));
     263         496 :   if (matctx->defctx->n && !matctx->defctx->ref) {
     264          80 :     PetscCall(VecConjugate(x));
     265          80 :     k = matctx->defctx->n;
     266          80 :     lambda = matctx->lambda;
     267          80 :     eig = matctx->defctx->eig;
     268          80 :     PetscCall(PetscMalloc2(k,&h,k,&alphaC));
     269          80 :     PetscCall(BVDotVec(matctx->defctx->W,x,h));
     270         176 :     for (i=0;i<k;i++) alphaC[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
     271         176 :     for (i=0;i<k;i++) h[i] *= alphaC[i]/(1.0-alphaC[i]);
     272          80 :     PetscCall(BVMultVec(matctx->defctx->V,1.0,1.0,x,h));
     273          80 :     PetscCall(PetscFree2(h,alphaC));
     274          80 :     PetscCall(VecConjugate(x));
     275             :   }
     276         496 :   PetscFunctionReturn(PETSC_SUCCESS);
     277             : }
     278             : 
     279          10 : static PetscErrorCode MatDestroy_NEPDeflationNE(Mat M)
     280             : {
     281          10 :   NEP_NEDEF_MATSHELL *matctx;
     282             : 
     283          10 :   PetscFunctionBegin;
     284          10 :   PetscCall(MatShellGetContext(M,&matctx));
     285          10 :   PetscCall(VecDestroy(&matctx->w[0]));
     286          10 :   PetscCall(VecDestroy(&matctx->w[1]));
     287          10 :   PetscCall(PetscFree(matctx));
     288          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     289             : }
     290             : 
     291           0 : static PetscErrorCode MatCreateVecs_NEPDeflationNE(Mat M,Vec *right,Vec *left)
     292             : {
     293           0 :   NEP_NEDEF_MATSHELL *matctx;
     294             : 
     295           0 :   PetscFunctionBegin;
     296           0 :   PetscCall(MatShellGetContext(M,&matctx));
     297           0 :   PetscCall(MatCreateVecs(matctx->F,right,left));
     298           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     299             : }
     300             : 
     301          10 : static PetscErrorCode NEPDeflationNEFunctionCreate(NEP_NEDEF_CTX defctx,NEP nep,Mat F,Mat P,Mat J,KSP ksp,PetscBool isJ,Mat *Mshell)
     302             : {
     303          10 :   NEP_NEDEF_MATSHELL *matctx;
     304          10 :   PetscInt           nloc,mloc;
     305             : 
     306          10 :   PetscFunctionBegin;
     307             :   /* Create mat shell */
     308          10 :   PetscCall(PetscNew(&matctx));
     309          10 :   PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
     310          10 :   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Mshell));
     311          10 :   matctx->F   = F;
     312          10 :   matctx->P   = P;
     313          10 :   matctx->J   = J;
     314          10 :   matctx->isJ = isJ;
     315          10 :   matctx->ksp = ksp;
     316          10 :   matctx->defctx = defctx;
     317          10 :   matctx->lambda = PETSC_MAX_REAL;
     318          10 :   PetscCall(MatCreateVecs(F,&matctx->w[0],NULL));
     319          10 :   PetscCall(VecDuplicate(matctx->w[0],&matctx->w[1]));
     320          10 :   PetscCall(MatShellSetOperation(*Mshell,MATOP_MULT,(void(*)(void))MatMult_NEPDeflationNE));
     321          10 :   PetscCall(MatShellSetOperation(*Mshell,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_NEPDeflationNE));
     322          10 :   PetscCall(MatShellSetOperation(*Mshell,MATOP_SOLVE,(void(*)(void))MatSolve_NEPDeflationNE));
     323          10 :   PetscCall(MatShellSetOperation(*Mshell,MATOP_SOLVE_TRANSPOSE,(void(*)(void))MatSolveTranspose_NEPDeflationNE));
     324          10 :   PetscCall(MatShellSetOperation(*Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_NEPDeflationNE));
     325          10 :   PetscCall(MatShellSetOperation(*Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_NEPDeflationNE));
     326          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     327             : }
     328             : 
     329          13 : static PetscErrorCode NEPDeflationNERecoverEigenvectors(NEP_NEDEF_CTX defctx,Vec u,Vec w,PetscScalar lambda)
     330             : {
     331          13 :   PetscScalar    *h,*alpha,*eig;
     332          13 :   PetscInt       i,k;
     333             : 
     334          13 :   PetscFunctionBegin;
     335          13 :   if (w) PetscCall(VecConjugate(w));
     336          13 :   if (defctx->n && !defctx->ref) {
     337           4 :     eig = defctx->eig;
     338           4 :     k = defctx->n;
     339           4 :     PetscCall(PetscMalloc2(k,&h,k,&alpha));
     340           9 :     for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
     341           4 :     PetscCall(BVDotVec(defctx->V,u,h));
     342           9 :     for (i=0;i<k;i++) h[i] *= alpha[i];
     343           4 :     PetscCall(BVMultVec(defctx->W,-1.0,1.0,u,h));
     344           4 :     PetscCall(VecNormalize(u,NULL));
     345           4 :     if (w) {
     346           4 :       PetscCall(BVDotVec(defctx->W,w,h));
     347           9 :       for (i=0;i<k;i++) alpha[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
     348           9 :       for (i=0;i<k;i++) h[i] *= alpha[i];
     349           4 :       PetscCall(BVMultVec(defctx->V,-1.0,1.0,w,h));
     350           4 :       PetscCall(VecNormalize(w,NULL));
     351             :     }
     352           4 :     PetscCall(PetscFree2(h,alpha));
     353             :   }
     354          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     355             : }
     356             : 
     357           9 : static PetscErrorCode NEPDeflationNELocking(NEP_NEDEF_CTX defctx,Vec u,Vec w,PetscScalar lambda)
     358             : {
     359           9 :   PetscInt       n;
     360             : 
     361           9 :   PetscFunctionBegin;
     362           9 :   n = defctx->n++;
     363           9 :   defctx->eig[n] = lambda;
     364           9 :   PetscCall(BVInsertVec(defctx->V,n,u));
     365           9 :   PetscCall(BVInsertVec(defctx->W,n,w));
     366           9 :   PetscCall(BVSetActiveColumns(defctx->V,0,defctx->n));
     367           9 :   PetscCall(BVSetActiveColumns(defctx->W,0,defctx->n));
     368           9 :   PetscCall(BVBiorthonormalizeColumn(defctx->V,defctx->W,n,NULL));
     369           9 :   PetscFunctionReturn(PETSC_SUCCESS);
     370             : }
     371             : 
     372          13 : static PetscErrorCode NEPDeflationNESetRefine(NEP_NEDEF_CTX defctx,PetscBool ref)
     373             : {
     374          13 :   PetscFunctionBegin;
     375          13 :   defctx->ref = ref;
     376          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     377             : }
     378             : 
     379           5 : PetscErrorCode NEPSolve_SLP_Twosided(NEP nep)
     380             : {
     381           5 :   NEP_SLP        *ctx = (NEP_SLP*)nep->data;
     382           5 :   Mat            mF,mJ,M,Mt;
     383           5 :   Vec            u,r,t,w;
     384           5 :   BV             X,Y;
     385           5 :   PetscScalar    sigma,lambda,mu,im=0.0,mu2,im2;
     386           5 :   PetscReal      resnorm,resl;
     387           5 :   PetscInt       nconv,nconv2,i;
     388           5 :   PetscBool      skip=PETSC_FALSE,lock=PETSC_FALSE;
     389           5 :   NEP_NEDEF_CTX  defctx=NULL;    /* Extended operator for deflation */
     390             : 
     391           5 :   PetscFunctionBegin;
     392             :   /* get initial approximation of eigenvalue and eigenvector */
     393           5 :   PetscCall(NEPGetDefaultShift(nep,&sigma));
     394           5 :   if (!nep->nini) PetscCall(BVSetRandomColumn(nep->V,0));
     395           5 :   PetscCall(BVSetRandomColumn(nep->W,0));
     396           5 :   lambda = sigma;
     397           5 :   if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
     398           5 :   PetscCall(BVDuplicate(nep->V,&X));
     399           5 :   PetscCall(BVDuplicate(nep->W,&Y));
     400           5 :   PetscCall(NEPDeflationNECreate(nep,X,Y,nep->nev,&defctx));
     401           5 :   PetscCall(BVGetColumn(nep->V,0,&t));
     402           5 :   PetscCall(VecDuplicate(t,&u));
     403           5 :   PetscCall(VecDuplicate(t,&w));
     404           5 :   PetscCall(BVRestoreColumn(nep->V,0,&t));
     405           5 :   PetscCall(BVCopyVec(nep->V,0,u));
     406           5 :   PetscCall(BVCopyVec(nep->W,0,w));
     407           5 :   PetscCall(VecDuplicate(u,&r));
     408           5 :   PetscCall(NEPDeflationNEFunctionCreate(defctx,nep,nep->function,nep->function_pre?nep->function_pre:nep->function,NULL,ctx->ksp,PETSC_FALSE,&mF));
     409           5 :   PetscCall(NEPDeflationNEFunctionCreate(defctx,nep,nep->function,nep->function,nep->jacobian,NULL,PETSC_TRUE,&mJ));
     410           5 :   PetscCall(NEPSLPSetUpEPSMat(nep,mF,mJ,PETSC_FALSE,&M));
     411           5 :   PetscCall(NEPSLPSetUpEPSMat(nep,mF,mJ,PETSC_TRUE,&Mt));
     412           5 :   PetscCall(EPSSetOperators(ctx->eps,M,NULL));
     413           5 :   PetscCall(MatDestroy(&M));
     414           5 :   PetscCall(EPSSetOperators(ctx->epsts,Mt,NULL));
     415           5 :   PetscCall(MatDestroy(&Mt));
     416             : 
     417             :   /* Restart loop */
     418          30 :   while (nep->reason == NEP_CONVERGED_ITERATING) {
     419          25 :     nep->its++;
     420             : 
     421             :     /* form residual,  r = T(lambda)*u (used in convergence test only) */
     422          25 :     PetscCall(NEPDeflationNEComputeFunction(nep,mF,lambda));
     423          25 :     PetscCall(MatMultTranspose(mF,w,r));
     424          25 :     PetscCall(VecNorm(r,NORM_2,&resl));
     425          25 :     PetscCall(MatMult(mF,u,r));
     426             : 
     427             :     /* convergence test */
     428          25 :     PetscCall(VecNorm(r,NORM_2,&resnorm));
     429          25 :     resnorm = PetscMax(resnorm,resl);
     430          25 :     PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
     431          25 :     nep->eigr[nep->nconv] = lambda;
     432          25 :     if (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol) {
     433           9 :       if (nep->errest[nep->nconv]<=ctx->deftol && !defctx->ref && nep->nconv) {
     434           4 :         PetscCall(NEPDeflationNERecoverEigenvectors(defctx,u,w,lambda));
     435           4 :         PetscCall(VecConjugate(w));
     436           4 :         PetscCall(NEPDeflationNESetRefine(defctx,PETSC_TRUE));
     437           4 :         PetscCall(MatMultTranspose(mF,w,r));
     438           4 :         PetscCall(VecNorm(r,NORM_2,&resl));
     439           4 :         PetscCall(MatMult(mF,u,r));
     440           4 :         PetscCall(VecNorm(r,NORM_2,&resnorm));
     441           4 :         resnorm = PetscMax(resnorm,resl);
     442           4 :         PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
     443           4 :         if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
     444           5 :       } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
     445             :     }
     446          16 :     if (lock) {
     447           9 :       lock = PETSC_FALSE;
     448           9 :       skip = PETSC_TRUE;
     449           9 :       PetscCall(NEPDeflationNERecoverEigenvectors(defctx,u,w,lambda));
     450           9 :       PetscCall(NEPDeflationNELocking(defctx,u,w,lambda));
     451           9 :       PetscCall(NEPDeflationNESetRefine(defctx,PETSC_FALSE));
     452           9 :       PetscCall(BVInsertVec(nep->V,nep->nconv,u));
     453           9 :       PetscCall(BVInsertVec(nep->W,nep->nconv,w));
     454           9 :       PetscCall(VecConjugate(w));
     455           9 :       nep->nconv = nep->nconv + 1;
     456             :     }
     457          25 :     PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
     458          25 :     if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));
     459             : 
     460          25 :     if (nep->reason == NEP_CONVERGED_ITERATING) {
     461          20 :       if (!skip) {
     462             :         /* evaluate T(lambda) and T'(lambda) */
     463          16 :         PetscCall(NEPDeflationNEComputeFunction(nep,mF,lambda));
     464          16 :         PetscCall(NEPDeflationNEComputeFunction(nep,mJ,lambda));
     465          16 :         PetscCall(EPSSetInitialSpace(ctx->eps,1,&u));
     466          16 :         PetscCall(EPSSetInitialSpace(ctx->epsts,1,&w));
     467             : 
     468             :         /* compute new eigenvalue correction mu and eigenvector approximation u */
     469          16 :         PetscCall(EPSSolve(ctx->eps));
     470          16 :         PetscCall(EPSSolve(ctx->epsts));
     471          16 :         PetscCall(EPSGetConverged(ctx->eps,&nconv));
     472          16 :         PetscCall(EPSGetConverged(ctx->epsts,&nconv2));
     473          16 :         if (!nconv||!nconv2) {
     474           0 :           PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", inner iteration failed, stopping solve\n",nep->its));
     475           0 :           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
     476           0 :           break;
     477             :         }
     478          16 :         PetscCall(EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL));
     479          16 :         for (i=0;i<nconv2;i++) {
     480          16 :           PetscCall(EPSGetEigenpair(ctx->epsts,i,&mu2,&im2,w,NULL));
     481          16 :           if (SlepcAbsEigenvalue(mu-mu2,im-im2)/SlepcAbsEigenvalue(mu,im)<nep->tol*1000) break;
     482             :         }
     483          16 :         if (i==nconv2) {
     484           0 :           PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", inner iteration failed, stopping solve\n",nep->its));
     485           0 :           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
     486           0 :           break;
     487             :         }
     488             : 
     489          16 :         mu = 1.0/mu;
     490          16 :         PetscCheck(PetscAbsScalar(im)<PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex eigenvalue approximation - not implemented in real scalars");
     491             :       } else {
     492           4 :         nep->its--;  /* do not count this as a full iteration */
     493             :         /* use second eigenpair computed in previous iteration */
     494           4 :         PetscCall(EPSGetConverged(ctx->eps,&nconv));
     495           4 :         if (nconv>=2 && nconv2>=2) {
     496           4 :           PetscCall(EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL));
     497           4 :           PetscCall(EPSGetEigenpair(ctx->epsts,1,&mu2,&im2,w,NULL));
     498           4 :           mu = 1.0/mu;
     499             :         } else {
     500           0 :           PetscCall(BVSetRandomColumn(nep->V,nep->nconv));
     501           0 :           PetscCall(BVSetRandomColumn(nep->W,nep->nconv));
     502           0 :           PetscCall(BVCopyVec(nep->V,nep->nconv,u));
     503           0 :           PetscCall(BVCopyVec(nep->W,nep->nconv,w));
     504           0 :           mu = lambda-sigma;
     505             :         }
     506             :         skip = PETSC_FALSE;
     507             :       }
     508             :       /* correct eigenvalue */
     509          20 :       lambda = lambda - mu;
     510             :     }
     511             :   }
     512           5 :   PetscCall(VecDestroy(&u));
     513           5 :   PetscCall(VecDestroy(&w));
     514           5 :   PetscCall(VecDestroy(&r));
     515           5 :   PetscCall(MatDestroy(&mF));
     516           5 :   PetscCall(MatDestroy(&mJ));
     517           5 :   PetscCall(BVDestroy(&X));
     518           5 :   PetscCall(BVDestroy(&Y));
     519           5 :   PetscCall(NEPDeflationNEDestroy(defctx));
     520           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     521             : }

Generated by: LCOV version 1.14