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             :    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 :   NEP_RESOLVENT_MATSHELL *ctx;
      31           4 :   NEP                    nep;
      32           4 :   PetscInt               i,inside=1;
      33           4 :   PetscScalar            alpha;
      34           4 :   Vec                    x,y,z,w;
      35             : 
      36           4 :   PetscFunctionBegin;
      37           4 :   PetscCall(MatShellGetContext(M,&ctx));
      38           4 :   nep = ctx->nep;
      39           4 :   w = nep->work[0];
      40           4 :   z = nep->work[1];
      41           4 :   if (((PetscObject)v)->id != ctx->vid || ((PetscObject)v)->state != ctx->vstate) {
      42           2 :     PetscCall(PetscArrayzero(ctx->dots_avail,ctx->nep->nconv));
      43           2 :     PetscCall(PetscObjectGetId((PetscObject)v,&ctx->vid));
      44           2 :     PetscCall(PetscObjectStateGet((PetscObject)v,&ctx->vstate));
      45             :   }
      46           4 :   PetscCall(VecSet(r,0.0));
      47          36 :   for (i=0;i<nep->nconv;i++) {
      48          32 :     if (ctx->rg) PetscCall(RGCheckInside(ctx->rg,1,&nep->eigr[i],&nep->eigi[i],&inside));
      49          32 :     if (inside>=0) {
      50          32 :       PetscCall(BVGetColumn(nep->V,i,&x));
      51          32 :       PetscCall(BVGetColumn(nep->W,i,&y));
      52          32 :       PetscCall(NEPApplyJacobian(nep,nep->eigr[i],x,z,w,NULL));
      53          32 :       if (!ctx->dots_avail[i]) {
      54          16 :         PetscCall(VecDot(v,y,&ctx->dots[i]));
      55          16 :         ctx->dots_avail[i] = PETSC_TRUE;
      56             :       }
      57          32 :       if (!ctx->nfactor_avail[i]) {
      58           8 :         PetscCall(VecDot(w,y,&ctx->nfactor[i]));
      59           8 :         ctx->nfactor_avail[i] = PETSC_TRUE;
      60             :       }
      61          32 :       alpha = ctx->dots[i]/(ctx->nfactor[i]*(ctx->omega-nep->eigr[i]));
      62          32 :       PetscCall(VecAXPY(r,alpha,x));
      63          32 :       PetscCall(BVRestoreColumn(nep->V,i,&x));
      64          32 :       PetscCall(BVRestoreColumn(nep->W,i,&y));
      65             :     }
      66             :   }
      67           4 :   PetscFunctionReturn(PETSC_SUCCESS);
      68             : }
      69             : 
      70           1 : static PetscErrorCode MatDestroy_Resolvent(Mat M)
      71             : {
      72           1 :   NEP_RESOLVENT_MATSHELL *ctx;
      73             : 
      74           1 :   PetscFunctionBegin;
      75           1 :   if (M) {
      76           1 :     PetscCall(MatShellGetContext(M,&ctx));
      77           1 :     PetscCall(PetscFree4(ctx->nfactor,ctx->nfactor_avail,ctx->dots,ctx->dots_avail));
      78           1 :     PetscCall(PetscFree(ctx));
      79             :   }
      80           1 :   PetscFunctionReturn(PETSC_SUCCESS);
      81             : }
      82             : 
      83             : /*@
      84             :    NEPApplyResolvent - Applies the resolvent T^{-1}(z) to a given vector.
      85             : 
      86             :    Collective
      87             : 
      88             :    Input Parameters:
      89             : +  nep   - eigensolver context obtained from NEPCreate()
      90             : .  rg    - optional region
      91             : .  omega - value where the resolvent must be evaluated
      92             : -  v     - input vector
      93             : 
      94             :    Output Parameter:
      95             : .  r     - result vector
      96             : 
      97             :    Notes:
      98             :    The resolvent T^{-1}(z) = sum_i (z-lambda_i)^{-1}*x_i*y_i' is evaluated at
      99             :    z=omega and the matrix-vector multiplication r = T^{-1}(omega)*v is computed.
     100             :    Vectors x_i and y_i are right and left eigenvectors, respectively, normalized
     101             :    so that y_i'*T'(lambda_i)*x_i=1. The sum contains only eigenvectors that have
     102             :    been previously computed with NEPSolve(), and if a region rg is given then only
     103             :    those corresponding to eigenvalues inside the region are considered.
     104             : 
     105             :    Level: intermediate
     106             : 
     107             : .seealso: NEPGetLeftEigenvector(), NEPSolve()
     108             : @*/
     109           4 : PetscErrorCode NEPApplyResolvent(NEP nep,RG rg,PetscScalar omega,Vec v,Vec r)
     110             : {
     111           4 :   NEP_RESOLVENT_MATSHELL *ctx;
     112             : 
     113           4 :   PetscFunctionBegin;
     114           4 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     115          16 :   PetscValidLogicalCollectiveScalar(nep,omega,3);
     116           4 :   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
     117           4 :   PetscValidHeaderSpecific(r,VEC_CLASSID,5);
     118           4 :   NEPCheckSolved(nep,1);
     119             : 
     120           4 :   PetscCall(PetscLogEventBegin(NEP_Resolvent,nep,0,0,0));
     121           4 :   if (!nep->resolvent) {
     122           1 :     PetscCall(PetscNew(&ctx));
     123           1 :     ctx->nep = nep;
     124           1 :     PetscCall(PetscCalloc4(nep->nconv,&ctx->nfactor,nep->nconv,&ctx->nfactor_avail,nep->nconv,&ctx->dots,nep->nconv,&ctx->dots_avail));
     125           1 :     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nep->nloc,nep->nloc,nep->n,nep->n,ctx,&nep->resolvent));
     126           1 :     PetscCall(MatShellSetOperation(nep->resolvent,MATOP_MULT,(void(*)(void))MatMult_Resolvent));
     127           1 :     PetscCall(MatShellSetOperation(nep->resolvent,MATOP_DESTROY,(void(*)(void))MatDestroy_Resolvent));
     128           3 :   } else PetscCall(MatShellGetContext(nep->resolvent,&ctx));
     129           4 :   PetscCall(NEPComputeVectors(nep));
     130           4 :   PetscCall(NEPSetWorkVecs(nep,2));
     131           4 :   ctx->rg    = rg;
     132           4 :   ctx->omega = omega;
     133           4 :   PetscCall(MatMult(nep->resolvent,v,r));
     134           4 :   PetscCall(PetscLogEventEnd(NEP_Resolvent,nep,0,0,0));
     135           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     136             : }