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

Generated by: LCOV version 1.14