function [out1, out2, out3, out4, out5] = guptri(A, B, epsu, gap, zero)
%
% NAME
%      guptri - Reduction to Generalized Upper Triangular Form
%  
% SYNOPSIS
%      guptri(A, B)
%      guptri(A, B, EPSU)
%      guptri(A, B, EPSU, GAP)
%      guptri(A, B, EPSU, GAP, ZERO)
%  
%      With the output formats:
%  
%                    kstr = guptri(...)
%                  [S, T] = guptri(...)
%            [S, T, kstr] = guptri(...)
%            [S, T, P, Q] = guptri(...)
%      [S, T, P, Q, kstr] = guptri(...)
%  
% DESCRIPTION
%      GUPTRI reduces the m by n pencil A - lambda B to generalized
%      upper  triangular  (guptri)  form  S  - lambda T via unitary
%      equivalence transformations P'*(A -  lambda   B)*Q.   guptri
%      form is described as follows:
%  
%          ( art  *   *   *   *  )      ( brt  *   *   *   *  )
%          (  0  azr  *   *   *  )      (  0  bzr  *   *   *  )
%      a = (  0   0  afn  *   *  ), b = (  0   0  bfn  *   *  )
%          (  0   0   0  ain  *  )      (  0   0   0  bin  *  )
%          (  0   0   0   0  alt )      (  0   0   0   0  blt )
%  
%      the diagonal blocks describe the  kronecker  canonical  form
%      (kcf) of the pencil a - lambda b as follows:
%  
%      art - lambda brt  has all right singular structure
%      azr - lambda bzr  has all jordan structure for 0 eigenvalue
%      afn - lambda bfn  has all jordan structure for finite
%                        nonzero eigenvalues
%      ain - lambda bin  has all jordan structure for infinite
%                        eigenvalue
%      alt - lambda blt  has all left singular structure
%  
%  
%      [S, T, P, Q, kstr] = guptri(A, B, EPSU,  GAP,  ZERO)  stores
%      the  transformation  matrices  in P and Q and the upper tri-
%      angular matrices P'*A*Q and P'*B*Q in S and T, respectively.
%      The  complete  kronecker  structure  is revealed by kstr, as
%      described in example below.
%  
%      EPSU is the relative uncertainty in data (should be at least
%      about  eps).  Used  to make rank decisions in singular value
%      decompositions. By default EPSU = 1E-8.
%  
%      GAP is used in singular value  decomposition  to  make  rank
%      descisions  by  searching for adjacent singular values whose
%      ratios exceeds GAP. GAP should at least be 1  and  nominally
%      1000.  By default GAP = 1000.
%  
%      ZERO is used as a logical variable which when it is true  (=
%      nonzero)  forces  guptri  to  zero out small singular values
%      during the reduction process. ZERO = 0 (false) will  give  a
%      reduced pencil made by true equivalence transformations, but
%      the pencil will not exactly have the structure described  by
%      kstr. By default ZERO = 1.
%  
% EXAMPLE
%  
%                     A = [0 1 0; 0 0 2];
%                     B = [0 0 0; 0 0 3];
%                     [S,T,P,Q,kstr]= guptri(A,B)
%  
%                     S =
%                          0  2  0
%                          0  0 -1
%  
%                     T =
%                          0  3  0
%                          0  0  0
%  
%                     P =
%                          0  1
%                          1  0
%  
%                     Q =
%                         -1  0  0
%                          0  0 -1
%                          0  1  0
%  
%                     kstr =
%                          1 -1  1  0 -1  1 -1
%                          0 -1  1  0 -1  1 -1
%  
%      The first part of kstr (the columns to the left of the first -1)
%      shall be interpreted as follows:
%        - The number of Li blocks is kstr(1,i+1) - kstr(2,i+1).
%        - The number of Ji(0) blocks is kstr(2,i) - kstr(1,i+1).
%      The second part of kstr contains the corresponding information
%      about Li' and Ni blocks and the third contains information
%      about the size of the regular part.
%      In the example, there is 1 L0 block, 1 N1 block, and a 1 by 1
%      regular part.
%  
% SEE ALSO
%      jnf
% 


% Parameter check:
% Check that number of input parameters are 2-5 and number of
% output parameters are 1-5.
% Initialize epsu, gap, and izero with default values if they 
% are not submitted and chach that they are scalar values 
% otherwise.
%
  if nargin < 2 
    error('Not enough input arguments.');
  elseif nargin > 5 
    error('Too many input arguments.');
  elseif nargout > 5
    error('Too many output arguments.');
  elseif min(size(A) == size(B))==0
    error('Matrix dimensions must agree.');
  elseif nargin > 3 & length(epsu) > 1
    error('Variable EPSU must be a scalar.');
  elseif nargin > 4 & length(gap) > 1
    error('Variable GAP must be a scalar.');
  elseif nargin > 5 & length(izero) > 1
    error('Variable IZERO must be a scalar.');
  elseif nargin < 5
    zero = 1;               % Use default value for izero
    if nargin < 4
      gap = 1000;           % Use default value for gap
      if nargin < 3
	epsu = sqrt(eps);   % Use default value for epsu
      end
    end
  end
  
  %
  % Error check for and initialization for other variables.
  %
  [m,n] = size(A);
  maxmn = max(m,n);
  minmn = min(m,n);
  wsize = 2*(2*(maxmn*maxmn) + m*n +minmn*minmn + 12*maxmn +3*minmn + 1);
  work = zeros(wsize,1);
  %
  % Store A and B as complex vectors with real and imaginary parts
  % interleaved (also if the matrices in fact are real).
  %
  Ac(1:2:2*m*n) = reshape(real(A),m*n,1);
  Ac(2:2:2*m*n) = reshape(imag(A),m*n,1);
  Bc(1:2:2*m*n) = reshape(real(B),m*n,1);
  Bc(2:2:2*m*n) = reshape(imag(B),m*n,1);
  
  [Pc, Qc, kstr, info] = mguptri(Ac,Bc,m,n,epsu,gap,zero,work,wsize);
  
  %
  %  Matrices S, T, P, Q constructed from real/complex interleved vectors.
  %
  S = reshape(Ac(1:2:2*m*n),m,n) + sqrt(-1)*reshape(Ac(2:2:2*m*n),m,n);
  T = reshape(Bc(1:2:2*m*n),m,n) + sqrt(-1)*reshape(Bc(2:2:2*m*n),m,n);
  P = reshape(Pc(1:2:2*m*m),m,m) + sqrt(-1)*reshape(Pc(2:2:2*m*m),m,m);
  Q = reshape(Qc(1:2:2*n*n),n,n) + sqrt(-1)*reshape(Qc(2:2:2*n*n),n,n);
  
  lkstr = max(max([find(kstr(1,:)),find(kstr(2,:)),find(kstr(3,:)),find(kstr(4,:))]));
  
  kstr = kstr(1:2,1:lkstr);
  if info ~= 0
    warning('INFO non-zero on return from Fortran-guptri (continuing)'); 
  end		
  
  
  if nargout > 1
    out1 = S;
    out2 = T;
  end
  if nargout > 3
    out3 = P;
    out4 = Q;
  end
  if nargout == 0 | nargout == 1
    out1 = kstr;
  elseif nargout == 3
    out3 = kstr;
  elseif nargout == 5
    out5 = kstr;
  end
  
