parameter (n=5)
      real*8 A(n,n), S(n,n), w(14*n), shift
      integer iw(5*n), i, j
 
      do j=1, n
         do i=j, n
            A(i,j) = 1.0D0/dfloat(i+j)
         enddo
      enddo

      print*, 'Enter shift'
      read*, shift

      do j=1, n
         A(j,j) = A(j,j) -shift
      enddo

      call DSYPHI(n,A,S,w,iw)
      
      do i=1, n
         print*, (S(i,j),j=1,i)
      enddo
 
      end

********************************************************************
*     This is our routine for calculating varphi(A) for a symmetric
*     matrix A.
*     Please consult the following paper for details.
*
*     Ya Yan Lu, Computing a Matrix Function for Exponential 
*                Integrators , preprint
*
*     available at: http://math.cityu.edu.hk/~mayylu
*
*     It involves the following steps:
*     1. reduce A to tridiagonal matrix T by Householder reflections.
*     2. calculate the LARGEST eigenvalue lambda_1
*     3. Evaluate \varphi(T) using three different formulas
*     4. evaluate \varphi(A) = Q varphi(T) Q'
*
*     n -- input integer, A is nxn symmetric matrix
*     A -- the symmetric matrix given by its lower triangular part
*     S -- varphi(A), also symmetric and given for its lower triangular part 
*     w -- work space, double precision
*     iw -- work space, integer
******************************************************************
      subroutine DSYPHI(n,A,S,w,iw)
      integer n
      integer iw(5*n)
      real*8 A(n,n), S(n,n), w(14*n)
*
      call DSYPHW(n,A,S,w(1), w(n+1), w(2*n+1),w(3*n+1),iw)
      return
      end

********************************************************************
*
      subroutine DSYPHW(n,A,S,d,e,tau,w,iw)
      implicit none
      integer n
      integer iw(5*n)
      real*8 A(n,n),S(n,n),d(n),e(n-1),tau(n), w(11*n)
      real*8 lams,vl,vu,abstol
      integer lwork,info,mfound,nsplit,il,iu

      lwork = 11*n

*     Reduce the matrix A to tridiagonal form by LAPACK
*     routine xSYTRD.

      call DSYTRD( 'L', n, A, n, d, e, tau, w, lwork, info )
*      print*, 'Best lwork = ', w(1)
*      print*, 'DSYTRD info = ', info

*     Find the largest and smallest eigenvalues by LAPACK
*     routine xSTEBZ.

      vl = 0.0D0
      vu = 0.0D0
      il = n
      iu = n
      abstol = 0.0D0
*      abstol = 2*DLAMCH('S')

      call DSTEBZ( 'I', 'B', n, vl, vu, il, iu, abstol, d, e,
     $     mfound, nsplit, w, iw, iw(n+1), w(n+1), iw(2*n+1),
     $     info )


*      print*, 'DSTEBZ info = ', info
      lams = w(1)
      print*, 'The largest eigenvalue is ', lams

*     Now, we calculate the varphi(T) based on lambda_1

      if ( lams .LT. -1.0D0 ) then
         print*, 'calling DSTPH1'
         call DSTPH1(n,d,e,S,w,lams,iw)
      else if (lams. LE. 0.0D0) then
         print*, 'calling DSTPH2'
         call DSTPH2(n,d,e,S,w,lams,iw)
      else
         print*, 'calling DSTPH3'
         call DSTPH3(n,d,e,S,w,lams,iw)
      endif

*      print*, 'Computing varphiT(T) done!'

*     Given symmetric matrix S (lower triangular), we calculate QSQ^t
*     Q is given as in the xSYTRD, stored in elementary reflectors.

      call APPLIQ(n,A,tau,S, w, w(n+1))
*      print*, 'AppliQ done!'

      return
      end

************************************************************
*     This is our program for calculating varphi(T) based on
*     varphi(T) = T^{-1}(exp(T)-I), we use rational Chebyshev
*     Pade approximation for exp(T) = exp(lam1) exp(T - lam1*I).
*     The coefficients here are for x positive
*
*     exp(-x) = alpha(0) + Re [ SUM alpha(j)/(x - theta(j))]
*
*     Thus for x < 0, 
*     exp(x)  = alpha(0) - Re [ SUM alpha(j)/(x + theta(j)) ].
*     This is just the opposite sign as used in the paper.
*     
*     The diagonal of T is the vector d,
*     the upper /lower diagonal of T is the vector e
*
*     T = diag(d(1:n)) + diag(e(1:n-1),1)+diag(e(1:n-1),-1)
*
*************************************************************     

      subroutine DSTPH1(n,d,e,B,w,lam1,iw)
      implicit none
      integer n, i, j, k, info, iw(n)
      real*8 d(n), e(n-1), B(n,n), alpha0, w(11*n), lam1, s
      complex*16 alpha(7), theta(7), sc, t

      alpha0  = 1.832174378247001445178D-14
      alpha(1)=( 5.575032388013989775304D+01, 2.0429467998037060112D+02)
      alpha(2)=(-9.386654897738416781398D+01,-9.1287299537244701102D+01)
      alpha(3)=( 4.699646418199852101538D+01, 1.1616718259322739191D+01)
      alpha(4)=(-9.614224197620382854756D+00, 2.6419587674922023894D+00)
      alpha(5)=( 7.527200775595744942100D-01,-6.7036694058767631296D-01)
      alpha(6)=(-1.887805062128241188701D-02, 3.4369583916841158280D-02)
      alpha(7)=( 1.430857612700760251239D-04,-2.8722086698954584504D-04)
      theta(1)=(-5.623142572742583221495D+00,-1.1940690463436614062D+00)
      theta(2)=(-5.089345060577111615918D+00,-3.5888240290260982584D+00)
      theta(3)=(-3.993369710574816069647D+00,-6.0048316422335547350D+00)
      theta(4)=(-2.269783829227002748129D+00,-8.4617379730382146007D+00)
      theta(5)=( 2.087586382547148287361D-01,-1.0991260561898794116D+01)
      theta(6)=( 3.703275049428656060120D+00,-1.3656371871480423129D+01)
      theta(7)=( 8.897773186474949120586D+00,-1.6630982619899044992D+01)

      s = dexp(lam1)*alpha0-1.0D0
      do j=1, n
         B(j,j) = s
      enddo
      do j=1,n-1
         do i=j+1,n
            B(i,j) = 0.0D0
         enddo
      enddo

      do k=1, 7
         sc = - dexp(lam1)*alpha(k)
         t = lam1 - theta(k)
         call BUPDAT(n,d,e, sc, t, B,w,w(2*n+1),
     #        w(4*n+1),w(6*n+1), w(8*n+1), iw)
      enddo

*     Finally, we calculate B:= T^{-1} B, but only the lower 
*     triangular part of B is really available. Besides, we need to
*     switch the sign.

      do j=2, n
         do i=1, j-1
            B(i,j) = - B(j,i)
         enddo
      enddo
      do j=1, n
         do i=j, n
            B(i,j) = - B(i,j)
         enddo
      enddo
      do j=1, n-1
         d(j) = -d(j)
         e(j) = -e(j)
      enddo
      d(n) = -d(n)
         
      call DPTTRF( n, d, e, info )
      call DPTTRS( n, n, d, e, B, n, info )  
      return
      end

************************************************************
*     This is our program for calculating varphi(T), when the 
*     largest eigenvalue of T is between [-1, 0], we simply use
*     the Chebyshev rational approximation for itself.
*
*     [ exp(-x) -1 ]/ x  = alpha(0) + 2 Re SUM alpha(j)/(x - theta(j))
*
*     Now, for x <0, 
*
*     [exp(x)-1]/x = -alpha(0) + 2 Re SUM alpha(j)/(x + theta(j))
*
*     The diagonal of T is the vector d,
*     the upper /lower diagonal of T is the vector e
*
*     T = diag(d(1:n)) + diag(e(1:n-1),1)+diag(e(1:n-1),-1)
*
*************************************************************     

      subroutine DSTPH2(n,d,e,B,w,lam1,iw)
      implicit none
      integer n, i, j, iw(n)
      real*8 d(n), e(n-1), B(n,n), w(11*n), lam1, a0
      complex*16 a(7), t(7), sc, tc

      a0 =  -6.8944296265527394984D-16
      a(1)=(-1.6598679663720768703D+01,-3.9025784287223383670D+01)
      a(2)=( 2.2963504666229092280D+01, 9.0186818220061090091D+00)
      a(3)=(-7.5350149609204679786D+00, 3.0951732326685966968D+00)
      a(4)=( 6.5440260116974146874D-01,-1.2832270822767467541D+00)
      a(5)=( 1.7992885377582909731D-02, 1.2021513848300774960D-01)
      a(6)=(-2.2224782352681356103D-03,-3.1546051373084948534D-03)
      a(7)=( 1.6950103692838164789D-05, 1.8407950619535128862D-05)

      t(1)=(-6.5586170606958520061D+00,-1.2541312162940416924D+00)
      t(2)=(-6.0329668674314355458D+00,-3.7686693138308950662D+00)
      t(3)=(-4.9527072954283340179D+00,-6.3037280204340004157D+00)
      t(4)=(-3.2515207076218489674D+00,-8.8794008802441251574D+00)
      t(5)=(-8.0133602893611439276D-01,-1.1529259279403978988D+01)
      t(6)=( 2.6587124072174283827D+00,-1.4320672417208411550D+01)
      t(7)=( 7.8095944003956373966D+00,-1.7439142275890278426D+01)

      do j=1,n
         do i=j,n
            B(i,j) = 0.0D0
         enddo
         B(j,j) = -a0
      enddo

      do j=1, 7
         tc =  - t(j)
         sc = 2.0D0*a(j)
         call BUPDAT(n,d,e, sc, tc, B,w,w(2*n+1),
     #        w(4*n+1),w(6*n+1), w(8*n+1), iw)
      enddo
      return
      end

************************************************************
*     This is our program for calculating varphi(T), when T
*     has a positive eigenvalue. From the approximation for 
*     exp(x), for x < 0:
*
*     exp(x)  = alpha(0) - Re [ SUM alpha(j)/(x + theta(j)) ].
*
*     which is just the opposite sign as used in the paper.
*
*     we approximate varphi(x) by 
*     
*     exp(lam1) Re [ -(alpha(j)/(lam1-theta(j))) 1/( x - lam1 + theta(j))
*
*     The diagonal of T is the vector d,
*     the upper /lower diagonal of T is the vector e
*
*     T = diag(d(1:n)) + diag(e(1:n-1),1)+diag(e(1:n-1),-1)
*
*************************************************************     

      subroutine DSTPH3(n,d,e,B,w,lam1,iw)
      implicit none
      integer n, i, j, iw(n)
      real*8 d(n), e(n-1), B(n,n), w(11*n), lam1
      complex*16 alpha(7), theta(7), sc, t

      alpha(1)=( 5.575032388013989775304D+01, 2.0429467998037060112D+02)
      alpha(2)=(-9.386654897738416781398D+01,-9.1287299537244701102D+01)
      alpha(3)=( 4.699646418199852101538D+01, 1.1616718259322739191D+01)
      alpha(4)=(-9.614224197620382854756D+00, 2.6419587674922023894D+00)
      alpha(5)=( 7.527200775595744942100D-01,-6.7036694058767631296D-01)
      alpha(6)=(-1.887805062128241188701D-02, 3.4369583916841158280D-02)
      alpha(7)=( 1.430857612700760251239D-04,-2.8722086698954584504D-04)
      theta(1)=(-5.623142572742583221495D+00,-1.1940690463436614062D+00)
      theta(2)=(-5.089345060577111615918D+00,-3.5888240290260982584D+00)
      theta(3)=(-3.993369710574816069647D+00,-6.0048316422335547350D+00)
      theta(4)=(-2.269783829227002748129D+00,-8.4617379730382146007D+00)
      theta(5)=( 2.087586382547148287361D-01,-1.0991260561898794116D+01)
      theta(6)=( 3.703275049428656060120D+00,-1.3656371871480423129D+01)
      theta(7)=( 8.897773186474949120586D+00,-1.6630982619899044992D+01)

      do j=1,n
         do i=j,n
            B(i,j) = 0.0D0
         enddo
      enddo
      
      do j=1, 7
         t = lam1 - theta(j)
         sc = - dexp(lam1)*alpha(j)/t
         call BUPDAT(n,d,e, sc, t, B,w,w(2*n+1),
     #        w(4*n+1),w(6*n+1), w(8*n+1), iw)
      enddo
      return
      end

**********************************************************************
*     Here we calculate B := B + alpha(k) * ( T - theta(k) I)^{-1}
*
      subroutine BUPDAT(n,d,e,alp,tht,B,zd,zl,zu,du2,bb,ipiv)
      implicit none
      integer n
      real*8 d(n), e(n-1), B(n,n)
      complex*16 alp, tht, zd(n), zl(n-1), zu(n-1), du2(n-2), bb(n)
      integer ipiv(n), info, nrhs, ldb, i, j

      do i=1, n
         zd(i) = d(i) - tht
      enddo
      do i=1, n-1
         zl(i) = e(i)
         zu(i) = e(i)
      enddo

      call ZGTTRF( n, zl, zd, zu, du2, ipiv, info )

      nrhs = 1
      ldb = n

      do j=1, n
         do i=1, n
            bb(i)=0.0D0
         enddo
         bb(j) = alp
         call ZGTTRS( 'N', n, nrhs, zl, zd, zu, du2, ipiv, bb, ldb,
     $     info )
         do i=j,n
            B(i,j) = B(i,j) + bb(i)
         enddo
      enddo

      return
      end

******************************************************************
*     This program requires 2n^3 flops.
*
*     SSYTRD reduces a real symmetric matrix A to real symmetric
*     tridiagonal form T by an orthogonal similarity transformation:
*     Q**T * A * Q = T, where Q is is represented as a product of 
*     elementary reflectors.
*
*     It takes: 4/3 n^3 operations.
*
*     If UPLO = 'L', the matrix Q is represented as a product of elementary
*     reflectors
*
*     Q = H(1) H(2) . . . H(n-1).
*
*     Each H(i) has the form
*
*     H(i) = I - tau * v * v'
*
*     where tau is a real scalar, and v is a real vector with
*     v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
*     and tau in TAU(i).
*
*     The contents of A on exit are illustrated by the following examples
*     with n = 5:
*
*     if UPLO = 'U':                       if UPLO = 'L':
*
*     (  d   e   v2  v3  v4 )              (  d                  )
*     (      d   e   v3  v4 )              (  e   d              )
*     (          d   e   v4 )              (  v1  e   d          )
*     (              d   e  )              (  v1  v2  e   d      )
*     (                  d  )              (  v1  v2  v3  e   d  )
*
*     where d and e denote diagonal and off-diagonal elements of T, and vi
*     denotes an element of the vector defining H(i).

      subroutine APPLIQ(n,A,tau,S,v,w)
      implicit none
      integer n
      real*8 A(n,n), tau(n-1), S(n,n), v(n), w(n)
      integer i, j, k
      real*8 vtw

      do k= n-2, 1, -1
*         print*, k

*     Redefine the vector v, as sqrt(tau)*v, H = I - v*v'

         v(k+1) = sqrt(tau(k))
         do i=k+2, n
            v(i) = A(i,k)*v(k+1)
         enddo

*     Calculate the lower left corner, dimension is: (n-k) x k
*     Total operations:  4 (n-k) k

         do j=1, k
            w(j) = 0.0D0
            do i=k+1, n
               w(j) = w(j) + v(i)*S(i,j)
            enddo
         enddo

         do j=1, k
            do i=k+1, n
               S(i,j) = S(i,j) - v(i)*w(j)
            enddo
         enddo

*     Calculate the lower right trailing block, only the lower
*     triangular part, domension (n-k) x (n-k)
*     Total operations: 4 (n-k)^2

*     Operations:  2 (n-k)^2
         do i=k+1, n
            w(i) = 0.0D0
            do j=k+1, i
               w(i) = w(i) + S(i,j)*v(j)
            enddo
            do j=i+1, n
               w(i) = w(i) + S(j,i)*v(j)
            enddo
         enddo

         vtw = 0.0D0
         do i=k+1, n
            vtw = vtw + v(i)*w(i)
         enddo
         
         vtw = 0.5D0*vtw
         do i=k+1,n
            w(i) = w(i) - vtw*v(i)
         enddo

*     Operations: 4/2 (n-k)^2 = 2(n-k)^2

         do j= k+1, n
            do i=j, n
               S(i,j) = S(i,j) - w(i)*v(j)-v(i)*w(j)
            enddo
         enddo
               
      enddo
      return
      end