function [ xmin, fmin ] = multiplier_method( f, g, h, x )
% �æ�”–@
  alpha = 10; beta = 1/4; eps = 1e-5;
  gs = size(g(x)); lambda = zeros(gs); slackr = 10*ones(gs);
  hs = size(h(x));     mu = zeros(hs); slacks = 10*ones(hs);
  c = 1e+24;
  options = optimset('TolFun',1e-20);
  
  for iteration = 1:1000
    L = @(x) augmented_Lagrange_function(f,g,h,x,lambda,mu,slackr,slacks);
    x = fminsearch(L, x, options);
    G = abs(g(x)); H = abs(max(h(x), -mu./slacks));
    Gmax = max(G); Hmax = max(H);
    if (isempty(G)); Gmax = 0.00; end
    if (isempty(H)); Hmax = 0.00; end
    if (Gmax < c && Hmax < c)
      c = max(Gmax, Hmax);
      if (c <= eps)
        xmin = x; fmin = f(x);
        return;
      end
      lambda = lambda + slackr.*g(x);
      mu = max(mu + slacks.*h(x), 0);
    end
    
    for i = find(G <= beta*c)
      slackr(i) = alpha*slackr(i);
    end
    for j = find(H <= beta*c)
      slacks(j) = alpha*slacks(j);
    end
  end
end
