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             :    Common subroutines for several PEP solvers
      12             : */
      13             : 
      14             : #include <slepc/private/pepimpl.h>    /*I "slepcpep.h" I*/
      15             : #include <slepcblaslapack.h>
      16             : 
      17             : /*
      18             :   Computes T_j=phy_idx(T) for a matrix T.
      19             :     Tp (if j>0) and Tpp (if j>1) are the evaluation
      20             :     of phy_(j-1) and phy(j-2)respectively.
      21             : */
      22          68 : PetscErrorCode PEPEvaluateBasisMat(PEP pep,PetscInt k,PetscScalar *T,PetscInt ldt,PetscInt idx,PetscScalar *Tpp,PetscInt ldtpp,PetscScalar *Tp,PetscInt ldtp,PetscScalar *Tj,PetscInt ldtj)
      23             : {
      24          68 :   PetscInt       i;
      25          68 :   PetscReal      *ca,*cb,*cg;
      26          68 :   PetscScalar    g,a;
      27          68 :   PetscBLASInt   k_,ldt_,ldtj_,ldtp_;
      28             : 
      29          68 :   PetscFunctionBegin;
      30          68 :   if (idx==0) {
      31          61 :     for (i=0;i<k;i++) {
      32          39 :       PetscCall(PetscArrayzero(Tj+i*ldtj,k));
      33          39 :       Tj[i+i*ldtj] = 1.0;
      34             :     }
      35             :   } else {
      36          46 :     PetscCall(PetscBLASIntCast(ldt,&ldt_));
      37          46 :     PetscCall(PetscBLASIntCast(ldtj,&ldtj_));
      38          46 :     PetscCall(PetscBLASIntCast(ldtp,&ldtp_));
      39          46 :     PetscCall(PetscBLASIntCast(k,&k_));
      40          46 :     ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
      41         126 :     for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
      42          46 :     if (idx>1) {
      43          65 :       for (i=0;i<k;i++) PetscCall(PetscArraycpy(Tj+i*ldtj,Tpp+i*ldtpp,k));
      44             :     }
      45          46 :     a = 1/ca[idx-1];
      46          46 :     g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
      47          46 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,Tp,&ldtp_,&g,Tj,&ldtj_));
      48         126 :     for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
      49             :   }
      50          68 :   PetscFunctionReturn(PETSC_SUCCESS);
      51             : }
      52             : 
      53             : /*
      54             :    PEPEvaluateBasis - evaluate the polynomial basis on a given parameter sigma
      55             : */
      56        6097 : PetscErrorCode PEPEvaluateBasis(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
      57             : {
      58        6097 :   PetscInt   nmat=pep->nmat,k;
      59        6097 :   PetscReal  *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
      60             : 
      61        6097 :   PetscFunctionBegin;
      62       19245 :   if (ivals) for (k=0;k<nmat;k++) ivals[k] = 0.0;
      63        6097 :   vals[0] = 1.0;
      64        6097 :   vals[1] = (sigma-b[0])/a[0];
      65             : #if !defined(PETSC_USE_COMPLEX)
      66             :   if (ivals) ivals[1] = isigma/a[0];
      67             : #endif
      68       15971 :   for (k=2;k<nmat;k++) {
      69        9874 :     vals[k] = ((sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2])/a[k-1];
      70        9874 :     if (ivals) vals[k] -= isigma*ivals[k-1]/a[k-1];
      71             : #if !defined(PETSC_USE_COMPLEX)
      72             :     if (ivals) ivals[k] = ((sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2])/a[k-1];
      73             : #endif
      74             :   }
      75        6097 :   PetscFunctionReturn(PETSC_SUCCESS);
      76             : }
      77             : 
      78             : /*
      79             :    PEPEvaluateBasisDerivative - evaluate the derivative of the polynomial basis on a given parameter sigma
      80             : */
      81         783 : PetscErrorCode PEPEvaluateBasisDerivative(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
      82             : {
      83         783 :   PetscInt       nmat=pep->nmat,k;
      84         783 :   PetscReal      *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
      85             : 
      86         783 :   PetscFunctionBegin;
      87         783 :   PetscCall(PEPEvaluateBasis(pep,sigma,isigma,vals,ivals));
      88        3141 :   for (k=nmat-1;k>0;k--) {
      89        2358 :     vals[k] = vals[k-1];
      90        2358 :     if (ivals) ivals[k] = ivals[k-1];
      91             :   }
      92         783 :   vals[0] = 0.0;
      93         783 :   vals[1] = vals[1]/a[0];
      94             : #if !defined(PETSC_USE_COMPLEX)
      95             :   if (ivals) ivals[1] = ivals[1]/a[0];
      96             : #endif
      97        2358 :   for (k=2;k<nmat;k++) {
      98        1575 :     vals[k] += (sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2];
      99        1575 :     if (ivals) vals[k] -= isigma*ivals[k-1];
     100        1575 :     vals[k] /= a[k-1];
     101             : #if !defined(PETSC_USE_COMPLEX)
     102             :     if (ivals) {
     103             :       ivals[k] += (sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2];
     104             :       ivals[k] /= a[k-1];
     105             :     }
     106             : #endif
     107             :   }
     108         783 :   PetscFunctionReturn(PETSC_SUCCESS);
     109             : }