eqlbrm41.f

      function eqlbrm_f_g_h(b, astar, bigb)
      external bigb
      eqlbrm_f_g_h = bigb(astar, b)
      end

      function eqlbrm_f_g(a, astar, biga, bigb, bstar, n)
      external biga
      bstar = argmax_2(astar, bigb, bstar, n)
      eqlbrm_f_g = biga(a, bstar)
      end

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

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

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

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

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

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

      subroutine deriv2_1_1(astar, biga, bigb, 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, biga, bigb,
     +bstar, bstar_g3, n, y, yprime)
      end

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

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

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

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

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

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

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

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

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

      function deriv1_1(astar, biga, bigb, 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, biga, bigb, bstar,
     +bstar_g1, n, y, deriv1_1)
      end

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

      function deriv1_2(astar, bigb, x)
      x_g2 = 1.0
      astar_g2 = 0.0
      call deriv1_2_adf2_g2(x, x_g2, astar, astar_g2, bigb, 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(gmbiga, gmbigb, astar, bstar, n)
      print *, astar, bstar
      end