%  [u, s, v, Gz] = zerodir( G, z, epp );
%           or
%  [u, s, v, Gz] = zerodir( G, z );
%
%  This function computes the singular value decomposition of a given
%  system at s = z. The outputs are the output directions U, the singular
%  values S along the diagonal, the input directions V and the transfer 
%  function G(z).
%
%  The last ny-r columns of u, corresponds to the singular output directions.
%  ny - number of outputs, r - is the rank of G(z), if z is an zero then r<ny.
%  This should is a safe way to compute the directions of a zero. However, 
%  the zeros must be known in advance and if more than one singular value is 
%  zero or close to zero it can be difficult to pick out the corrosponding
%  direction.
%
%  If G(s) has a zero on the same location as a pole, then the matrix (zI-A)
%  becomes singular (This is also true when the directions of the pole and 
%  the zero are not the same). The function then adds small singular values 
%  to the zero singular values of (zI-A). The singular values of (zI-A)
%  which are zero are set to the EPP times the smallest non-zero singular 
%  value. If the parameter EPP is not supplied in the call to the function
%  the default value is 1e-8.
%
%  The method with generalized eigenvalues does not suffer from this weakness.
%
%  See also:
%              OZDE, IZDE ans SZEROS.
%
%  Written by: Kjetil Havre 20/4 - 1995.
%  e-mail:     havre@kjemi.unit.no or skoge@kjemi.unit.no
%

%
%  Malfunctions: When the system has a zero at the same location of a pole,
%                the SVD function promts warning message about singular
%                matrix for G(z). The reason is that the matrix (zI-A)
%                becomes singular. This is corrected by adding small 
%                singular values for the zero singular values. See 
%                description above.
%

function [u, s, v, Gz] = zerodir( G, z, epp )

   if nargin < 3
      epp = 1e-8;
   end


   [A,B,C,D] = unpck( G );
   
   % May expect problems if Iz - A is singular. Have to look into this.
   H = z*eye(size(A)) - A;
   [U, S, V]  = svd( H );
   Si = zeros(size(S));
   nr = rank(S);
   [nxr, nx] = size( A );
   if nr < max(size(A))
      disp( ['Warning: Matrix (zI-A) is singular adds small singular values' ] );
      disp( ['         of size: ', num2str(S(nr,nr)*epp), ' to singular part.' ] );
      disp( ['         Result may be inaccurate, the result should be compared to'] );
      disp( ['         the result obtained with generalized eigenvalue methode.'] );
      for i=nr+1:nx,
         S(i,i) = S(nr,nr)*epp;
      end
   end
   for i=1:nx
      Si(i,i) = 1/S(i,i);
   end
   Hi = V*Si*U';
   Gz = C*Hi*B+D;
   [u, s, v] = svd( Gz );
end










