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          47 : 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          47 :   PetscInt       i;
      25          47 :   PetscReal      *ca,*cb,*cg;
      26          47 :   PetscScalar    g,a;
      27          47 :   PetscBLASInt   k_,ldt_,ldtj_,ldtp_;
      28             : 
      29          47 :   PetscFunctionBegin;
      30          47 :   if (idx==0) {
      31          45 :     for (i=0;i<k;i++) {
      32          30 :       PetscCall(PetscArrayzero(Tj+i*ldtj,k));
      33          30 :       Tj[i+i*ldtj] = 1.0;
      34             :     }
      35             :   } else {
      36          32 :     PetscCall(PetscBLASIntCast(ldt,&ldt_));
      37          32 :     PetscCall(PetscBLASIntCast(ldtj,&ldtj_));
      38          32 :     PetscCall(PetscBLASIntCast(ldtp,&ldtp_));
      39          32 :     PetscCall(PetscBLASIntCast(k,&k_));
      40          32 :     ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
      41          94 :     for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
      42          32 :     if (idx>1) {
      43          49 :       for (i=0;i<k;i++) PetscCall(PetscArraycpy(Tj+i*ldtj,Tpp+i*ldtpp,k));
      44             :     }
      45          32 :     a = 1/ca[idx-1];
      46          32 :     g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
      47          32 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,Tp,&ldtp_,&g,Tj,&ldtj_));
      48          94 :     for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
      49             :   }
      50          47 :   PetscFunctionReturn(PETSC_SUCCESS);
      51             : }
      52             : 
      53             : /*
      54             :    PEPEvaluateBasis - evaluate the polynomial basis on a given parameter sigma
      55             : */
      56        4409 : PetscErrorCode PEPEvaluateBasis(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
      57             : {
      58        4409 :   PetscInt   nmat=pep->nmat,k;
      59        4409 :   PetscReal  *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
      60             : 
      61        4409 :   PetscFunctionBegin;
      62       18507 :   if (ivals) for (k=0;k<nmat;k++) ivals[k] = 0.0;
      63        4409 :   vals[0] = 1.0;
      64        4409 :   vals[1] = (sigma-b[0])/a[0];
      65             : #if !defined(PETSC_USE_COMPLEX)
      66        4409 :   if (ivals) ivals[1] = isigma/a[0];
      67             : #endif
      68       10657 :   for (k=2;k<nmat;k++) {
      69        6248 :     vals[k] = ((sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2])/a[k-1];
      70        6248 :     if (ivals) vals[k] -= isigma*ivals[k-1]/a[k-1];
      71             : #if !defined(PETSC_USE_COMPLEX)
      72        6248 :     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        4409 :   PetscFunctionReturn(PETSC_SUCCESS);
      76             : }
      77             : 
      78             : /*
      79             :    PEPEvaluateBasisDerivative - evaluate the derivative of the polynomial basis on a given parameter sigma
      80             : */
      81         246 : PetscErrorCode PEPEvaluateBasisDerivative(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
      82             : {
      83         246 :   PetscInt       nmat=pep->nmat,k;
      84         246 :   PetscReal      *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
      85             : 
      86         246 :   PetscFunctionBegin;
      87         246 :   PetscCall(PEPEvaluateBasis(pep,sigma,isigma,vals,ivals));
      88         766 :   for (k=nmat-1;k>0;k--) {
      89         520 :     vals[k] = vals[k-1];
      90         520 :     if (ivals) ivals[k] = ivals[k-1];
      91             :   }
      92         246 :   vals[0] = 0.0;
      93         246 :   vals[1] = vals[1]/a[0];
      94             : #if !defined(PETSC_USE_COMPLEX)
      95         246 :   if (ivals) ivals[1] = ivals[1]/a[0];
      96             : #endif
      97         520 :   for (k=2;k<nmat;k++) {
      98         274 :     vals[k] += (sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2];
      99         274 :     if (ivals) vals[k] -= isigma*ivals[k-1];
     100         274 :     vals[k] /= a[k-1];
     101             : #if !defined(PETSC_USE_COMPLEX)
     102         274 :     if (ivals) {
     103         274 :       ivals[k] += (sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2];
     104         274 :       ivals[k] /= a[k-1];
     105             :     }
     106             : #endif
     107             :   }
     108         246 :   PetscFunctionReturn(PETSC_SUCCESS);
     109             : }