LCOV - code coverage report
Current view: top level - nep/interface - nepresolv.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 68 69 98.6 %
Date: 2021-08-02 00:35:43 Functions: 3 3 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             :    NEP routines related to resolvent T^{-1}(z) = sum_i (z-lambda_i)^{-1} x_i y_i'
      12             : */
      13             : 
      14             : #include <slepc/private/nepimpl.h>       /*I "slepcnep.h" I*/
      15             : 
      16             : typedef struct {
      17             :   NEP              nep;
      18             :   RG               rg;
      19             :   PetscScalar      omega;
      20             :   PetscScalar      *nfactor;         /* normalization factors y_i'*T'(lambda_i)*x_i */
      21             :   PetscBool        *nfactor_avail;
      22             :   PetscScalar      *dots;            /* inner products y_i'*v */
      23             :   PetscBool        *dots_avail;
      24             :   PetscObjectId    vid;
      25             :   PetscObjectState vstate;
      26             : } NEP_RESOLVENT_MATSHELL;
      27             : 
      28           4 : static PetscErrorCode MatMult_Resolvent(Mat M,Vec v,Vec r)
      29             : {
      30           4 :   PetscErrorCode         ierr;
      31           4 :   NEP_RESOLVENT_MATSHELL *ctx;
      32           4 :   NEP                    nep;
      33           4 :   PetscInt               i,inside=1;
      34           4 :   PetscScalar            alpha;
      35           4 :   Vec                    x,y,z,w;
      36             : 
      37           4 :   PetscFunctionBegin;
      38           4 :   ierr = MatShellGetContext(M,(void**)&ctx);CHKERRQ(ierr);
      39           4 :   nep = ctx->nep;
      40           4 :   w = nep->work[0];
      41           4 :   z = nep->work[1];
      42           4 :   if (((PetscObject)v)->id != ctx->vid || ((PetscObject)v)->state != ctx->vstate) {
      43           2 :     ierr = PetscArrayzero(ctx->dots_avail,ctx->nep->nconv);CHKERRQ(ierr);
      44           2 :     ierr = PetscObjectGetId((PetscObject)v,&ctx->vid);CHKERRQ(ierr);
      45           2 :     ierr = PetscObjectStateGet((PetscObject)v,&ctx->vstate);CHKERRQ(ierr);
      46             :   }
      47           4 :   ierr = VecSet(r,0.0);CHKERRQ(ierr);
      48          36 :   for (i=0;i<nep->nconv;i++) {
      49          32 :     if (ctx->rg) {
      50           0 :       ierr = RGCheckInside(ctx->rg,1,&nep->eigr[i],&nep->eigi[i],&inside);CHKERRQ(ierr);
      51             :     }
      52          32 :     if (inside>=0) {
      53          32 :       ierr = BVGetColumn(nep->V,i,&x);CHKERRQ(ierr);
      54          32 :       ierr = BVGetColumn(nep->W,i,&y);CHKERRQ(ierr);
      55          32 :       ierr = NEPApplyJacobian(nep,nep->eigr[i],x,z,w,NULL);CHKERRQ(ierr);
      56          32 :       if (!ctx->dots_avail[i]) {
      57          16 :         ierr = VecDot(v,y,&ctx->dots[i]);CHKERRQ(ierr);
      58          16 :         ctx->dots_avail[i] = PETSC_TRUE;
      59             :       }
      60          32 :       if (!ctx->nfactor_avail[i]) {
      61           8 :         ierr = VecDot(w,y,&ctx->nfactor[i]);CHKERRQ(ierr);
      62           8 :         ctx->nfactor_avail[i] = PETSC_TRUE;
      63             :       }
      64          32 :       alpha = ctx->dots[i]/(ctx->nfactor[i]*(ctx->omega-nep->eigr[i]));
      65          32 :       ierr = VecAXPY(r,alpha,x);CHKERRQ(ierr);
      66          32 :       ierr = BVRestoreColumn(nep->V,i,&x);CHKERRQ(ierr);
      67          32 :       ierr = BVRestoreColumn(nep->W,i,&y);CHKERRQ(ierr);
      68             :     }
      69             :   }
      70           4 :   PetscFunctionReturn(0);
      71             : }
      72             : 
      73           1 : static PetscErrorCode MatDestroy_Resolvent(Mat M)
      74             : {
      75           1 :   PetscErrorCode         ierr;
      76           1 :   NEP_RESOLVENT_MATSHELL *ctx;
      77             : 
      78           1 :   PetscFunctionBegin;
      79           1 :   if (M) {
      80           1 :     ierr = MatShellGetContext(M,(void**)&ctx);CHKERRQ(ierr);
      81           1 :     ierr = PetscFree4(ctx->nfactor,ctx->nfactor_avail,ctx->dots,ctx->dots_avail);CHKERRQ(ierr);
      82           1 :     ierr = PetscFree(ctx);CHKERRQ(ierr);
      83             :   }
      84           1 :   PetscFunctionReturn(0);
      85             : }
      86             : 
      87             : /*@
      88             :    NEPApplyResolvent - Applies the resolvent T^{-1}(z) to a given vector.
      89             : 
      90             :    Collective on nep
      91             : 
      92             :    Input Parameters:
      93             : +  nep   - eigensolver context obtained from NEPCreate()
      94             : .  rg    - optional region
      95             : .  omega - value where the resolvent must be evaluated
      96             : -  v     - input vector
      97             : 
      98             :    Output Parameter:
      99             : .  r     - result vector
     100             : 
     101             :    Notes:
     102             :    The resolvent T^{-1}(z) = sum_i (z-lambda_i)^{-1}*x_i*y_i' is evaluated at
     103             :    z=omega and the matrix-vector multiplication r = T^{-1}(omega)*v is computed.
     104             :    Vectors x_i and y_i are right and left eigenvectors, respectively, normalized
     105             :    so that y_i'*T'(lambda_i)*x_i=1. The sum contains only eigenvectors that have
     106             :    been previously computed with NEPSolve(), and if a region rg is given then only
     107             :    those corresponding to eigenvalues inside the region are considered.
     108             : 
     109             :    Level: intermediate
     110             : 
     111             : .seealso: NEPGetLeftEigenvector(), NEPSolve()
     112             : @*/
     113           4 : PetscErrorCode NEPApplyResolvent(NEP nep,RG rg,PetscScalar omega,Vec v,Vec r)
     114             : {
     115           4 :   PetscErrorCode         ierr;
     116           4 :   NEP_RESOLVENT_MATSHELL *ctx;
     117             : 
     118           4 :   PetscFunctionBegin;
     119           4 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     120           8 :   PetscValidLogicalCollectiveScalar(nep,omega,3);
     121           4 :   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
     122           4 :   PetscValidHeaderSpecific(r,VEC_CLASSID,5);
     123           4 :   NEPCheckSolved(nep,1);
     124             : 
     125           4 :   ierr = PetscLogEventBegin(NEP_Resolvent,nep,0,0,0);CHKERRQ(ierr);
     126           4 :   if (!nep->resolvent) {
     127           1 :     ierr = PetscNew(&ctx);CHKERRQ(ierr);
     128           1 :     ctx->nep = nep;
     129           1 :     ierr = PetscCalloc4(nep->nconv,&ctx->nfactor,nep->nconv,&ctx->nfactor_avail,nep->nconv,&ctx->dots,nep->nconv,&ctx->dots_avail);CHKERRQ(ierr);
     130           1 :     ierr = MatCreateShell(PetscObjectComm((PetscObject)nep),nep->nloc,nep->nloc,nep->n,nep->n,ctx,&nep->resolvent);CHKERRQ(ierr);
     131           1 :     ierr = MatShellSetOperation(nep->resolvent,MATOP_MULT,(void(*)(void))MatMult_Resolvent);CHKERRQ(ierr);
     132           1 :     ierr = MatShellSetOperation(nep->resolvent,MATOP_DESTROY,(void(*)(void))MatDestroy_Resolvent);CHKERRQ(ierr);
     133             :   } else {
     134           3 :     ierr = MatShellGetContext(nep->resolvent,(void**)&ctx);CHKERRQ(ierr);
     135             :   }
     136           4 :   ierr = NEPComputeVectors(nep);CHKERRQ(ierr);
     137           4 :   ierr = NEPSetWorkVecs(nep,2);CHKERRQ(ierr);
     138           4 :   ctx->rg    = rg;
     139           4 :   ctx->omega = omega;
     140           4 :   ierr = MatMult(nep->resolvent,v,r);CHKERRQ(ierr);
     141           4 :   ierr = PetscLogEventEnd(NEP_Resolvent,nep,0,0,0);CHKERRQ(ierr);
     142           4 :   PetscFunctionReturn(0);
     143             : }
     144             : 

Generated by: LCOV version 1.14