function xmin = feasibleDircMin(fun,x0,A,b)
%FEASIBLEMIN feasible direction method with equality constraint.607 hw2
%  xmin = feasibleDircMin(fun,x0,A,b) solves the problem
%               min  fun(x)  
%   ( P )       s.t. A*x = b
%         with initial guess x0.
%  It returns the minimizer xmin and displays the iterations.
%  
%  For equality constraints,we can transform constraint problem ( P )
%  to equivalent unconstraint problem 
%    ( P')      min  fun(xbar + Z * v)
%               where  xbar is one sulotion of  A*x = b;
%                      Z is one base of the null space of A,
%  then solve ( P') using desired methods for unconstraint problem.
%  Here we use fminunc.
%
%   Examples:
%      A=[1,1,1];b=3;x0=[-1;5;-1];
%      feasibleDircMin('myf',x0,A,b)

% Reference: 
% [1]http://www.cs.umd.edu/users/oleary/a607/607constrfdhand.pdf
% [2]S. G. Nash ,A Sofer:Linear and Nonlinear Programming  pages 505-507

%Author: Geping Liu 10-14-06 (geping@math.umd.edu)

global xbar Z funHandler;

% Check constraint qualification
[n,m] = size(A');
if (m >= n) || (rank(A) ~= m),
    disp('Error:constraint qualification not satisfied');
    xmin =[];
    return;
end

% Transform ( P ) to ( P')
funHandler = fun;

[Q,R] = qr(A');
xbar = A\b;       % one particular solution of A*x=b;
%%%%%%%%%%%%%%%%% Correction to next line 10-31-06
% Changed from Z = Q(:,n-m:end);
Z = Q(:,m+1:end); % one base of null space of A
%%%%%%%%%%%%%%%%% 
v0 = Z \( x0 - xbar);

% Solve ( P')
options = optimset('Display','iter');
vmin =fminunc(@F,v0,options);

% Transform back to x variable
xmin = xbar + Z * vmin; 

end


function z = F(v)
%F(v) := fun(xbar + Z * v)

global xbar Z funHandler;

z = feval(funHandler, xbar + Z * v);

end