C     =================================================================
C     File: bukm-prob.f
C     =================================================================

C     ******************************************************************
C     ******************************************************************

      subroutine inip(n,x,l,u,m,lambda,equatn,linear,coded,checkder)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer m,n
      logical checkder

C     ARRAY ARGUMENTS
      logical coded(10),equatn(*),linear(*)
      double precision l(*),lambda(*),u(*),x(*)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision beta,ctol,inbeta,inMtot,infixc,Mtot

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision adv(nassetsmax),c(nassetsmax),eps(nassetsmax),
     +        pi(nassetsmax),theta(nscenmax,nassetsmax),
     +        thetave(nassetsmax)

C     LOCAL SCALARS
      integer i,j

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /tracosts/ adv,c,eps,pi,beta,Mtot
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     ------------------------------------------------------------------
C     SET PROBLEM DATA
C     ------------------------------------------------------------------

C     AZN.L BARC.L KGF.L LLOY.L MKS.L TSCO.L VOD.L
C     (March 19, 2004 to February 12, 2008)

C     Constants

      nassets = 8
      nscen   = 1000
      r       = int( 0.99d0 * nscen )
      ctol    = 0.05d0  ! Proportion between 0 and 1 (0.05 means 5%)

      beta    = inbeta  ! Use 1 for linear transactional costs
      Mtot    = inMtot  ! 10^4 for a small investor, 10^8 for a large investor

      do i = 1,nassets - 1
          c(i) = infixc ! 0.01 for a small investor, 0.001 for a large one
      end do

      c(nassets) = 0.0d0

C     Scenarios and objective function coefficients:

      open(20,file='theta-010days.txt')

      do i = 1,nscen
          read(20,*) (theta(i,j),j=1,nassets)
      end do

      read(20,*) (thetave(j),j=1,nassets)

      close(20)

C     Pi:
      pi(1)      =   1.9550000000000000D+03
      pi(2)      =   4.5575000000000000D+02
      pi(3)      =   1.3630000000000001D+02
      pi(4)      =   4.2150000000000000D+02
      pi(5)      =   4.2000000000000000D+02
      pi(6)      =   4.0600000000000000D+02
      pi(7)      =   1.7769999999999999D+02
      pi(8)      =   1.0000000000000000D+00
 
C     ADV:
      adv(1)     =   6.9676110000000000D+06
      adv(2)     =   3.5557812799999997D+07
      adv(3)     =   1.8997996800000001D+07
      adv(4)     =   2.9665247500000000D+07
      adv(5)     =   1.3074735400000000D+07
      adv(6)     =   3.2649118399999999D+07
      adv(7)     =   2.4735484609999999D+08
      adv(8)     =   1.0000000000000000D+00
 
C     Spreads:
      eps(1)     =   2.0000000000000000D+00
      eps(2)     =   7.5000000000000000D-01
      eps(3)     =   5.9999999999999998D-01
      eps(4)     =   2.5000000000000000D-01
      eps(5)     =   5.0000000000000000D-01
      eps(6)     =   2.5000000000000000D-01
      eps(7)     =   2.0000000000000001D-01
      eps(8)     =   0.0000000000000000D+00

C     ------------------------------------------------------------------
C     SET MODEL DATA
C     ------------------------------------------------------------------

C     Number of variables

      n = nassets

C     Initial point

      do i = 1,n - 1
          x(i) = 0.5d0 / ( n - 1 )
      end do

      x(n) = 0.5d0

C     Lower and upper bounds

      do i = 1,n
          l(i) = 0.0d0
          u(i) = 1.0d+20
      end do

C     Number of constraints (equalities plus inequalities)

      m = 1 + r

C     Lagrange multipliers approximation, equatn and linear arrays.

      do i = 1,m
          lambda(i) = 0.0d0
      end do

      equatn(1) = .true.
      do i = 2,m
          equatn(i) = .false.
      end do

      if ( beta .eq. 1.0d0 ) then
          do i = 1,m
              linear(i) = .true.
          end do

      else
          linear(1) = .true.
          do i = 2,m
              linear(i) = .false.
          end do
      end if

C     Indicate which subroutines did you code.

      coded( 1) = .true.  ! evalf    
      coded( 2) = .true.  ! evalg   
      coded( 3) = .true.  ! evalh   
      coded( 4) = .true.  ! evalc   
      coded( 5) = .true.  ! evaljac
      coded( 6) = .true.  ! evalhc  
      coded( 7) = .false. ! evalfc
      coded( 8) = .false. ! evalgjac
      coded( 9) = .false. ! evalhl
      coded(10) = .false. ! evalhlp

C     Set checkder

      checkder = .false.

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalf(n,x,f,flag)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer j
      double precision t

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      external t

      flag = 0

      f = 0.0d0

      do j = 1,nassets
          f = f - thetave(j) * x(j) 
          if ( usetc ) f = f + t(j,x(j))
      end do

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalg(n,x,g,flag)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer flag,n

C     ARRAY ARGUMENTS
      double precision g(n),x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer j
      double precision dtdx

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      external dtdx

      flag = 0

      do j = 1,nassets
          g(j) = - thetave(j) 
          if ( usetc ) g(j) = g(j) + dtdx(j,x(j))
      end do
 
      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalh(n,x,hlin,hcol,hval,hnnz,flag)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer flag,n,hnnz

C     ARRAY ARGUMENTS
      integer hcol(*),hlin(*)
      double precision hval(*),x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer j
      double precision dtdxx

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      external dtdxx

      flag = 0

      if ( usetc ) then
          hnnz = n

          do j = 1,nassets
              hlin(j) = j
              hcol(j) = j
              hval(j) = dtdxx(j,x(j))
          end do

      else
          hnnz = 0
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalc(n,x,ind,c,flag)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer ind,flag,n
      double precision c

C     ARRAY ARGUMENTS
      double precision x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer i,j
      double precision t

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      external t

      flag = 0

      if ( ind .eq. 1 ) then
          c = - 1.0d0

          do i = 1,nassets
              c = c + x(i)
          end do

      else if ( 2 .le. ind .and. ind .le. r + 1 ) then
          i = find(ind - 1)

          c = 1.0d0 - ctol
          do j = 1,nassets
              c = c - theta(i,j) * x(j) 
              if ( usetc ) c = c + t(j,x(j))
          end do
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljac(n,x,ind,jcvar,jcval,jcnnz,flag)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer flag,ind,jcnnz,n

C     ARRAY ARGUMENTS
      integer jcvar(n)
      double precision x(n),jcval(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer i,j
      double precision dtdx

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      external dtdx

      flag = 0

      if ( ind .eq. 1 ) then
          jcnnz = nassets

          do j = 1,nassets
              jcvar(j) = j
              jcval(j) = 1.0d0
          end do

      else if ( 2 .le. ind .and. ind .le. r + 1 ) then
          i = find(ind - 1)

          jcnnz = nassets

          do j = 1,nassets
              jcvar(j) = j
              jcval(j) = - theta(i,j)
              if ( usetc ) jcval(j) = jcval(j) + dtdx(j,x(j))
          end do
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalhc(n,x,ind,hclin,hccol,hcval,hcnnz,flag)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer flag,hcnnz,ind,n

C     ARRAY ARGUMENTS
      integer hccol(*),hclin(*)
      double precision hcval(*),x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer j
      double precision dtdxx

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      external dtdxx

      flag = 0

      if ( ind .eq. 1 ) then
          hcnnz = 0

      else if ( 2 .le. ind .and. ind .le. r + 1 ) then

          if ( usetc ) then
              hcnnz = n

              do j = 1,nassets
                  hclin(j) = j
                  hccol(j) = j
                  hcval(j) = dtdxx(j,x(j))
              end do
          else
              hcnnz = 0
          end if
      end if

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalfc(n,x,f,m,c,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,m,n
      double precision f

C     ARRAY ARGUMENTS
      double precision c(m),x(n)

      flag = - 1

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,jcnnz,m,n

C     ARRAY ARGUMENTS
      integer jcfun(*),jcvar(*)
      double precision g(n),jcval(*),x(n)

      flag = - 1

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalhl(n,x,m,lambda,scalef,scalec,hllin,hlcol,hlval,
     +hlnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,hlnnz,m,n
      double precision scalef

C     ARRAY ARGUMENTS
      integer hlcol(*),hllin(*)
      double precision hlval(*),lambda(m),scalec(m),x(n)

      flag = - 1

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalhlp(n,x,m,lambda,sf,sc,p,hp,goth,flag)

      implicit none

C     SCALAR ARGUMENTS
      logical goth
      integer flag,m,n
      double precision sf

C     ARRAY ARGUMENTS
      double precision hp(n),lambda(m),p(n),sc(m),x(n)

      flag = - 1

      end

C     ******************************************************************
C     ******************************************************************

      subroutine endp(n,x,l,u,m,lambda,equatn,linear)

      implicit none

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision l(n),lambda(m),u(n),x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer cmin,i,j
      double precision tmp

C     LOCAL ARRAYS
      double precision c(nscenmax)

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      double precision t
      external t

C     Testing feasibility of the VaR constraint

      do i = 1,nscen
          c(i) = 1.0d0 - ctol
          do j = 1,nassets
              c(i) = c(i) - theta(i,j) * x(j) 
              if ( usetc ) c(i) = c(i) + t(j,x(j))
          end do
      end do

      do i = 1,nscen
          cmin = i
          do j = i + 1,nscen
              if ( c(j) .lt. c(cmin) ) then
                  cmin = j
              end if
          end do
          tmp     = c(i)
          c(i)    = c(cmin)
          c(cmin) = tmp
      end do

      do i = 1,nscen
          write(*,*) i,c(i)
      end do

      write(*,*) 'Value of the r-th smallest constraint: ',c(r)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compfind(n,x)

C     PARAMETERS
      integer nassetsmax,nscenmax
      parameter ( nassetsmax =    8 )
      parameter ( nscenmax   = 1000 )

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      double precision x(n)

C     COMMON SCALARS
      logical usetc
      integer algtype,nassets,nscen,r
      double precision ctol,inbeta,inMtot,infixc

C     COMMON ARRAYS
      integer find(nscenmax)
      double precision theta(nscenmax,nassetsmax),thetave(nassetsmax)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision fval(nscenmax)

C     COMMON BLOCKS
      common /probdata/ theta,thetave,find,ctol,nassets,nscen,r
      common /metadata/ algtype,usetc,inbeta,inMtot,infixc

C     EXTERNAL FUNCTIONS
      double precision t
      external t

      do i = 1,nscen
          find(i) = i
          fval(i) = 1.0d0 - ctol
          do j = 1,nassets
              fval(i) = fval(i) - theta(i,j) * x(j) 
              if ( usetc ) fval(i) = fval(i) + t(j,x(j))
          end do
      end do

      call rmin(nscen,r,fval,find)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine rmin(nscen,r,fval,find)

C     SCALAR ARGUMENTS
      integer nscen,r

C     ARRAY ARGUMENTS
      integer find(nscen)
      double precision fval(nscen)

C     LOCAL SCALARS
      integer classify,ileft,ind,iright,rtrial
      double precision drand,seed

C     EXTERNAL FUNCTIONS
      external drand

      seed = 123456.0d0

      ileft  = 1
      iright = nscen

 10   continue

      ind = ileft + int( drand(seed) * ( iright - ileft + 1 ) )
      ind = max( ileft, min( ind, iright ) )

      call iswap(find(ind),find(ileft))
      call dswap(fval(ind),fval(ileft))

      rtrial = classify(n,ileft,iright,fval,find)

      if ( rtrial .lt. r ) then
          ileft  = rtrial + 1
          go to 10

      else if ( rtrial .gt. r ) then

          if ( rtrial .eq. iright ) then
              call iswap(find(ileft),find(iright))
              call dswap(fval(ileft),fval(iright))
              rtrial = rtrial - 1
          end if

          iright = rtrial
          go to 10
      end if

      end
 
C     ******************************************************************
C     ******************************************************************

      integer function classify(nscen,ileft,iright,fval,find)

C     SCALAR ARGUMENTS
      integer nscen,ileft,iright

C     ARRAY ARGUMENTS
      integer find(nscen)
      double precision fval(nscen)

C     LOCAL SCALARS
      integer i,ilast

      i     = ileft + 1
      ilast = iright

 10   if ( i .le. ilast ) then
          if ( fval(i) .le. fval(ileft) ) then
              i = i + 1
          else
              call iswap(find(i),find(ilast))
              call dswap(fval(i),fval(ilast))
              ilast = ilast - 1
          end if
          go to 10
      end if

      classify = i - 1

      return

      end

C     ******************************************************************
C     ******************************************************************

      subroutine dswap(a,b)

C     SCALAR ARGUMENTS
      double precision a,b

C     LOCAL SCALARS
      double precision c

      c = a
      a = b
      b = c

      end

C     ******************************************************************
C     ******************************************************************

      subroutine iswap(a,b)

C     SCALAR ARGUMENTS
      integer a,b

C     LOCAL SCALARS
      integer c

      c = a
      a = b
      b = c

      end


C     ******************************************************************
C     ******************************************************************

      logical function ieqsets(n,a,b)

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      integer a(n),b(n)

C     LOCAL SCALARS
      logical found
      integer i,j

      ieqsets = .true.

      i = 1

 10   if ( ieqsets .and. i .le. n ) then

          found = .false.

          j = 1

 20       if ( .not. found .and. j .le. n ) then

              if ( a(i) .eq. b(j) ) found = .true.

              j = j + 1
              go to 20

          end if

          if ( .not. found ) ieqsets = .false.

          i = i + 1
          go to 10

      end if   

      return

      end

C     ******************************************************************
C     ******************************************************************

      double precision function t(j,xj)

C     PARAMETERS
      integer nassetsmax
      parameter ( nassetsmax =    8 )

C     SCALAR ARGUMENTS
      integer j
      double precision xj

C     COMMON SCALARS
      double precision beta,Mtot

C     COMMON ARRAYS
      double precision adv(nassetsmax),c(nassetsmax),eps(nassetsmax),
     +        pi(nassetsmax)

C     LOCAL SCALARS
      double precision h,g,c1,c2,c3,v1,v2

C     COMMON BLOCKS
      common /tracosts/ adv,c,eps,pi,beta,Mtot

      if ( xj .le. 1.0d-16 ) then
          t = 0.0d0
          return
      end if

      c1 = eps(j) / pi(j)

      c2 = 1.0d0 / ( ( adv(j) / Mtot ) * pi(j) ) ** beta

      c3 = 0.5d0 * c1 + 100.0d0 ** beta * c1

      v1 = xj ** beta

      h  = c3 * c2 * v1

      g  = c1 * c2 * 10.0d0 ** beta * v1

      v2 = h + g + c(j)

      t  = v2 * xj

      return

      end

C     ******************************************************************
C     ******************************************************************

      double precision function dtdx(j,xj)

C     PARAMETERS
      integer nassetsmax
      parameter ( nassetsmax =    8 )

C     SCALAR ARGUMENTS
      integer j
      double precision xj

C     COMMON SCALARS
      double precision beta,Mtot

C     COMMON ARRAYS
      double precision adv(nassetsmax),c(nassetsmax),eps(nassetsmax),
     +        pi(nassetsmax)

C     LOCAL SCALARS
      double precision h,g,c1,c2,c3,c4,c5,v1,v2,dhdx,dgdx,dv1dx,dv2dx

C     COMMON BLOCKS
      common /tracosts/ adv,c,eps,pi,beta,Mtot

      if ( xj .le. 1.0d-16 ) then
          dtdx = 0.0d0
          return
      end if

      c1    = eps(j) / pi(j)

      c2    = 1.0d0 / ( ( adv(j) / Mtot ) * pi(j) ) ** beta

      c3    = 0.5d0 * c1 + 100.0d0 ** beta * c1

      c4    = c3 * c2

      c5    = c1 * c2 * 10.0d0 ** beta

      v1    = xj ** beta

      dv1dx = beta * xj ** ( beta - 1.0d0 )

      h     = c4 * v1

      dhdx  = c4 * dv1dx

      g     = c5 * v1

      dgdx  = c5 * dv1dx

      v2    = h + g + c(j)

      dv2dx = dhdx + dgdx

      dtdx  = dv2dx * xj + v2 

      return

      end

C     ******************************************************************
C     ******************************************************************

      double precision function dtdxx(j,xj)

C     PARAMETERS
      integer nassetsmax
      parameter ( nassetsmax =    8 )

C     SCALAR ARGUMENTS
      integer j
      double precision xj

C     COMMON SCALARS
      double precision beta,Mtot

C     COMMON ARRAYS
      double precision adv(nassetsmax),c(nassetsmax),eps(nassetsmax),
     +        pi(nassetsmax)

C     LOCAL SCALARS
      double precision h,g,c1,c2,c3,c4,c5,dhdx,dgdx,dv1dx,dv2dx,dhdxx,
     +        dgdxx,dv1dxx,dv2dxx

C     COMMON BLOCKS
      common /tracosts/ adv,c,eps,pi,beta,Mtot

      if ( xj .le. 1.0d-16 ) then
          dtdxx = 0.0d0
          return
      end if

      c1     = eps(j) / pi(j)

      c2     = 1.0d0 / ( ( adv(j) / Mtot ) * pi(j) ) ** beta

      c3     = 0.5d0 * c1 + 100.0d0 ** beta * c1

      c4     = c3 * c2

      c5     = c1 * c2 * 10.0d0 ** beta

      v1     = xj ** beta

      dv1dx  = beta * xj ** ( beta - 1.0d0 )

      dv1dxx = beta * ( beta - 1.0d0 ) * xj ** ( beta - 2.0d0 )

      h      = c4 * v1

      dhdx   = c4 * dv1dx

      dhdxx  = c4 * dv1dxx

      g      = c5 * v1

      dgdx   = c5 * dv1dx

      dgdxx  = c5 * dv1dxx

      v2     = h + g + c(j)

      dv2dx  = dhdx + dgdx

      dv2dxx = dhdxx + dgdxx

      dtdx   = dv2dx * xj + v2 

      dtdxx  = dv2dxx * xj + 2.0d0 * dv2dx

      return

      end
