*************************************************************************
*               gulf research and development function
* more, garbow, and hillstrom, acm toms vol. 7 no. 1 (march 1981) 17-41
*************************************************************************

      subroutine getfun( x, n, f, m, ftf, fj, lfj, g, mode)

      implicit double precision (a-h,o-z)

      integer            n, m, lfj, mode

      double precision   x(n), f(m), ftf, fj(lfj,n), g(n)

      integer            nprob, nprobs, nstart, nstrts
      common /PROBLM/    nprob, nprobs, nstart, nstrts

      integer            nout
      common /IOUNIT/    nout

      logical            lf, lj

      integer            na, nb, nc, nd, nt, nh

      integer            i, j

      double precision   ai, av, bi, ci, d1, d2, d3, ei, ti, yi
      double precision   x1, x2, x3, x1inv

      intrinsic          abs, dble, exp, log

      double precision   ddot

      double precision   two3rd
      common /PARAM1/    two3rd
      save   /PARAM1/

      double precision   zero, one, point1, twnty5, fifty
      parameter         (zero = 0.d0, one = 1.d0)
      parameter         (point1 = .01d0, twnty5 = 25.d0, fifty = 50.d0)

*=======================================================================

      if (mode .eq.  0)  goto    20
      if (mode .eq. -1)  goto    10
      if (mode .eq. -2)  goto    30

      x1 = x(1)
      x2 = x(2)
      x3 = x(3)
      if (x1 .eq. zero) then
        write( 6, 1)
 1      format(/' +++ singularity in gulf function evaluation'/)
        stop
      endif
      x1inv = one / x1

      na = mode / 1000
      nt = mode - na*1000
      nb = nt / 100
      nh = nt - nb*100
      nc = nh / 10
      nd = nh - nc*10

      lf = (na .ne. 0) .or. (nb .ne. 0) .or. (nd .ne. 0)
      lj = (nc .ne. 0) .or. (nd .ne. 0)

      if (lf .and. lj)  goto 300
      if (lf)           goto 100
      if (lj)           goto 200

*-----------------------------------------------------------------------

   10 continue

      nprobs = 1
      nstrts = 1

      n     =  3
      m     = 10

      two3rd = 2.d0 / 3.d0

      if (nout .gt. 0)  write( nout, 9999) n, m

      return

*-----------------------------------------------------------------------

   20 continue

      x(1) =  5.d0
      x(2) =  2.5d0
      x(3) =  0.15d0

      return

*-----------------------------------------------------------------------

   30 continue

      x(1) =  50.d0
      x(2) =  25.d0
      x(3) =  1.5d0

      ftf = zero

      return

*-----------------------------------------------------------------------

 100  continue

      do 110 i = 1, m
        ti       = dble(i)*point1
        yi       = twnty5 + (-fifty*log(ti)) ** two3rd
        ai       = yi - x2
        f(i)     = exp(-((abs(ai)**x3)/x1)) - ti
 110  continue

      if (nb .ne. 0)  ftf = ddot( m, f, 1, f, 1)

      return
c
 200  continue
c
      do 210 i = 1, m
        ti = dble(i)*point1
        yi = twnty5 + (-fifty*log(ti)) ** two3rd
        ai = yi - x2
        av = abs(ai)
        bi = av ** x3
        ci = bi*x1inv
        ei = exp(-ci)
        d1 =   ci * x1inv
        d2 =   x3 * x1inv * av**(x3 - one)
        d3 = - log(av) * ci
        fj(i,1) = d1 * ei
        if (ai .ge. zero)  fj(i,2) =  d2 * ei
        if (ai .lt. zero)  fj(i,2) = -d2 * ei
        fj(i,3) = d3 * ei
 210  continue
      return

 300  continue

      do 310 i = 1, m
        ti = dble(i)*point1
        yi = twnty5 + (-fifty*log(ti)) ** two3rd
        ai = yi - x2
        av = abs(ai)
        bi = av**x3
        ci = bi*x1inv
        ei = exp(-ci)
        d1 =   ci * x1inv
        d2 =   x3 * x1inv * av**(x3 - one)
        d3 = - log(av) * ci
        f(i) = ei - ti
        fj(i,1) = d1 * ei
        if (ai .ge. zero)  fj(i,2) =  d2 * ei
        if (ai .lt. zero)  fj(i,2) = -d2 * ei
        fj(i,3) = d3 * ei
 310  continue

      if (nb .ne. 0)  ftf = ddot( m, f, 1, f, 1)

      if (nd .eq. 0)  return

      do 400 j = 1, n
        g(j) = ddot( m, fj( 1, j), 1, f, 1)
 400  continue
      return

9999  format(/'1',70('=')//,
     *' gulf research and development function (more et al.)'//,
     *'        number of variables =', i4, '  (      3       )'/,
     *'        number of functions =', i4, '  ( >= 3, <= 100 )'//,
     *        ' ',70('=')/)
      end

************************************************************************
************************************************************************

      subroutine dfjdxk ( k, x, n, dfj, ldfj, m, nonzro)

      implicit double precision (a-h,o-z)

      integer            k, n, ldfj, m, nonzro(n)

      double precision   x(n), dfj(ldfj,n)

      intrinsic          abs, dble, exp, log

      double precision   two3rd
      common /PARAM1/    two3rd
      save   /PARAM1/

      double precision   zero, one, point1, twnty5, fifty
      parameter         (zero = 0.d0, one = 1.d0, two = 2.d0)
      parameter         (point1 = .01d0, twnty5 = 25.d0, fifty = 50.d0)

*=======================================================================

      x1 = x(1)
      x2 = x(2)
      x3 = x(3)

      x1inv = one/x1

      nonzro(1) = 1
      nonzro(2) = 1
      nonzro(3) = 1

      goto ( 100, 200, 300), k

 100  continue

      do 110 i = 1, m
        ti = dble(i)*point1
        yi = twnty5 + (-fifty*log(ti)) ** two3rd
        ai = yi - x2
        av = abs(ai)
        b1 = av**(x3-one)
        bi = av**x3
        gi = log(av)
        ci = bi * x1inv
        ei = exp(-ci)
        d1 =   ci * x1inv
        d2 =   x3 * x1inv * b1
        d3 = - log(av) * ci
        d11 = - two * ci * x1inv * x1inv
        d12 = - x3 * b1 * x1inv * x1inv
        d13 =   gi * ci * x1inv
        dfj(i,1) = (d11 + d1*d1) * ei
        if (ai .ge. zero)  dfj(i,2) =  (d12 + d2*d2) * ei
        if (ai .lt. zero)  dfj(i,2) = -(d12 + d1*d2) * ei
        dfj(i,3) = (d13 + d1*d3) * ei
 110  continue

      return

 200  continue

      do 210 i = 1, m
        ti = dble(i)*point1
        yi = twnty5 + (-fifty*log(ti)) ** two3rd
        ai = yi - x2
        av = abs(ai)
        b2 = av**(x3-two)
        b1 = av**(x3-one)
        bi = av**x3
        gi = log(av)
        ci = bi * x1inv
        ei = exp(-ci)
        d1 =   ci * x1inv
        d2 =   x3 * x1inv * b1
        d3 = - log(av) * ci
        d21 = - x3 * b1 * x1inv * x1inv
        d22 = - x3 * (x3 - one) * b2 * x1inv
        d23 =   b1 * x1inv * (one + x3 * gi)
        if (ai .ge. zero)  dfj(i,1) =  (d21 + d2*d1) * ei
        if (ai .lt. zero)  dfj(i,1) = -(d21 + d2*d1) * ei
        dfj(i,2) =  (d22 + d2*d2) * ei
        if (ai .ge. zero)  dfj(i,3) =  (d23 + d2*d3) * ei
        if (ai .lt. zero)  dfj(i,3) = -(d23 + d2*d3) * ei
 210  continue

      return

 300  continue

      do 310 i = 1, m
        ti = dble(i)*point1
        yi = twnty5 + (-fifty*log(ti)) ** two3rd
        ai = yi - x2
        av = abs(ai)
        b1 = av**(x3-one)
        bi = av**x3
        gi = log(av)
        ci = bi * x1inv
        ei = exp(-ci)
        d1 =   ci * x1inv
        d2 =   x3 * x1inv * b1
        d3 = - log(av) * ci
        d31 =   gi * ci * x1inv
        d32 =   b1 * x1inv * (one + x3 * gi)
        d33 = - gi * gi * ci
        dfj(i,1) =  (d31 + d3*d1) * ei
        if (ai .ge. zero)  dfj(i,2) =  (d32 + d3*d2) * ei
        if (ai .lt. zero)  dfj(i,2) = -(d32 + d3*d2) * ei
        dfj(i,3) =  (d33 + d3*d3) * ei
 310  continue

      return

      end

************************************************************************
************************************************************************

      subroutine dfkdij ( k, x, n, hess, lhess, linear)

      implicit double precision (a-h,o-z)

      logical            linear

      integer            k, n, lhess

      double precision   x(n), hess(lhess,n)

      intrinsic          abs, dble, exp, log

      double precision   two3rd
      common /PARAM1/    two3rd
      save   /PARAM1/

      double precision   zero, one, point1, twnty5, fifty
      parameter         (zero = 0.d0, one = 1.d0, two = 2.d0)
      parameter         (point1 = .01d0, twnty5 = 25.d0, fifty = 50.d0)

*=======================================================================

      linear = .false.

      x1 = x(1)
      x2 = x(2)
      x3 = x(3)

      tk = dble(k)*point1
      yk = twnty5 + (-fifty*log(tk)) ** two3rd
      ak = yk - x2
      av = abs(ak)
      b2 = av**(x3 - two)
      b1 = av**(x3 - one)
      bk = av**x3
      if (x1 .eq. zero)  then
        write( 6, 1)
   1    format( ' second derivative does not exist')
        stop
      endif
      x1inv = one / x1
      ck = bk * x1inv
      ek = exp(-ck)
      gk = log(av)

      d1 =   ck * x1inv
      d2 =   x3 * x1inv * b1
      d3 = - gk * ck

      d11 = - two * ck * x1inv * x1inv
      d12 = - x3 * b1 * x1inv * x1inv
      d13 =   gk * ck * x1inv
      d22 = - x3 * (x3 - one) * b2 * x1inv
      d23 =   b1 * x1inv * (one + x3 * gk)
      d33 = - gk * gk * ck
      if (ak .ge. zero)  goto 100
        d2  = -d2
        d12 = -d12
        d23 = -d23

 100  continue

      hess( 1, 1) = (d11 + d1*d1) * ek
      hess( 1, 2) = (d12 + d1*d2) * ek
      hess( 1, 3) = (d13 + d1*d3) * ek
      hess( 2, 2) = (d22 + d2*d2) * ek
      hess( 2, 3) = (d23 + d2*d3) * ek
      hess( 3, 3) = (d33 + d3*d3) * ek
      hess( 2, 1) = hess( 1, 2)
      hess( 3, 1) = hess( 1, 3)
      hess( 3, 2) = hess( 2, 3)

      return
      end