function [v,info] = sylvest(g,d,h,tol)
%
% SYLVEST  Solves a Sylvester equation
%
%  solves the Sylvester equation
%         v = d + g v h
%  for v where both g and h must be  upper block triangular.
%  The output info is zero on a successful return. 
%  The input tol indicates when an element of g or h should be considered
%  zero.
%   
%
%   REFERENCES:
%         Golub, Nash and Van Loan(1979)
%         E.W. Anderson, L.P. Hansen, E.R McGrattan and T.J Sargent(1996)
%
%   PROGRAM  WRITTEN BY:  E. W. Anderson    Nov 26, 1994  (mex file only)
%   MODIFIED:             EWA               Dec 26, 1998  (m-file code added)   


%Find out dimensions

[m,n] = size(d);


%Initialize some variables

v = zeros(m,n);
w = eye(m);
i = 1;
temp = [];
info=0;


%First handle the i = 1 case outside the loop

if i< n,
  if abs(h(i+1,i)) < tol,
    v(:,i)= (w - g*h(i,i))\d(:,i);
    i = i+1;
  else
    A = [w-g*h(i,i)     (-g*h(i+1,i));...
       -g*h(i,i+1)    w-g*h(i+1,i+1)];
    C = [d(:,i); d(:,i+1)];
    X = A\C;
    v(:,i) = X(1:m,:);
    v(:,i+1) = X(m+1:2*m,  :); 
    i = i+2; 
  end
end


%Handle the rest of the matrix with the possible exception of i=n


while i<n,
  b= i-1;
  temp = [temp g*v(:,size(temp,2)+1:b)]; 
  if abs(h(i+1,i)) < tol,
    v(:,i) = (w - g*h(i,i))\(d(:,i) + temp*h(1:b,i));
    i = i+1;
  else
    A = [w - g*h(i,i)    (-g*h(i+1,i));   ...
         -g*h(i,i+1)    w - g*h(i+1,i+1)];
    C = [d(:,i) + temp*h(1:b,i);         ...
         d(:,i+1) + temp*h(1:b,i+1)];         
    X = A\C;
    v(:,i) = X(1:m,:);
    v(:,i+1) = X(m+1:2*m, :);
    i = i+2;
  end
end



%Handle the i = n case if i=n was not in a 2-2 block

if i==n,
  b = i-1;
  temp = [temp g*v(:,size(temp,2)+1:b)];
  v(:,i) = (w-g*h(i,i))\(d(:,i) + temp*h(1:b,i));
end