eqlbrm42.f

      function eqlbrm_f_g_h(b, astar)
      eqlbrm_f_g_h = gmbigb(astar, b)
      end

      function eqlbrm_f_g(a, astar, bstar, n)
      bstar = argmax_2(astar, bstar, n)
      eqlbrm_f_g = gmbiga(a, bstar)
      end

      function eqlbrm_f(astar, bstar, n)
      eqlbrm_f = argmax_1(astar, bstar, astar, n) - astar
      end

C     astar & bstar: guesses in, optimized values out
      subroutine eqlbrm(astar, bstar, n)
      astar = root_2(bstar, astar, n)
      end

      function root_1_1(astar, bstar, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_1_1(astar, bstar, n, x, y, yprime)
 1669 x = x - y / yprime
      root_1_1 = x
      end

      function root_1_2(astar, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_1_2(astar, x, y, yprime)
 1669 x = x - y / yprime
      root_1_2 = x
      end

      function root_2(bstar, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_2(bstar, n, x, y, yprime)
 1669 x = x - y / yprime
      root_2 = x
      end

      subroutine deriv2_1_1_adf3(x, astar, bstar, n, y)
      y = argmax_fprime_1(x, astar, bstar, n)
      end

      subroutine deriv2_1_1(astar, bstar, n, x, y, yprime)
      x_g3 = 1.0
      astar_g3 = 0.0
      bstar_g3 = 0.0
      call deriv2_1_1_adf3_g3(x, x_g3, astar, astar_g3, bstar, bstar_g3,
     +n, y, yprime)
      end

      subroutine deriv2_1_2_adf4(x, astar, y)
      y = argmax_fprime_2(x, astar)
      end

      subroutine deriv2_1_2(astar, x, y, yprime)
      x_g4 = 1.0
      astar_g4 = 0.0
      call deriv2_1_2_adf4_g4(x, x_g4, astar, astar_g4, y, yprime)
      end

      subroutine deriv2_2_adf5(x, bstar, n, y)
      y = eqlbrm_f(x, bstar, n)
      end

      subroutine deriv2_2(bstar, n, x, y, yprime)
      x_g5 = 1.0
      bstar_g5 = 1.0
      call deriv2_2_adf5_g5(x, x_g5, bstar, bstar_g5, n, y, yprime)
      end

      function argmax_fprime_1(x, astar, bstar, n)
      argmax_fprime_1 = deriv1_1(astar, bstar, n, x)
      end

      function argmax_fprime_2(x, astar)
      argmax_fprime_2 = deriv1_2(astar, x)
      end

      function argmax_1(astar, bstar, x0, n)
      argmax_1 = root_1_1(astar, bstar, x0, n)
      end

      function argmax_2(astar, x0, n)
      argmax_2 = root_1_2(astar, x0, n)
      end

      subroutine deriv1_1_adf1(x, astar, bstar, n, y)
      y = eqlbrm_f_g(x, astar, bstar, n)
      end

      function deriv1_1(astar, bstar, n, x)
      x_g1 = 1.0
      astar_g1 = 0.0
      bstar_g1 = 0.0
      call deriv1_1_adf1_g1(x, x_g1, astar, astar_g1, bstar, bstar_g1,
     +n, y, deriv1_1)
      end

      subroutine deriv1_2_adf2(x, astar, y)
      y = eqlbrm_f_g_h(x, astar)
      end

      function deriv1_2(astar, x)
      x_g2 = 1.0
      astar_g2 = 0.0
      call deriv1_2_adf2_g2(x, x_g2, astar, astar_g2, y, deriv1_2)
      end

      function gmbiga(a, b)
      price = 20 - 0.1*a - 0.1*b
      costs = a*(10 - 0.05*a)
      gmbiga = a*price - costs
      end

      function gmbigb(a, b)
      price = 20 - 0.1*b - 0.0999*a
      costs = b*(10.005 - 0.05*b)
      gmbigb = b*price - costs
      end

      program main
      read *, astar
      read *, bstar
      read *, n
      call eqlbrm(astar, bstar, n)
      print *, astar, bstar
      end