%   [I, O, err] = pdsvd(G, P, epp)
%
%   The ``Pole-Direction-SVD'' function computes the input and
%   output pole directions through the use of Singular Value 
%   Decomposition of G(p) where p is a pole. Note that this is an
%   ill-posed problem. Better way of calculating the pole directions 
%   are given in PDFST-function.
%   If largest singular value of G(p) is larger than EPP then the 
%   corrosponding input and output singular vectors are stored in
%   columns of I and O. Otherwise, zero vectors are stored in the 
%   corrosponding columns.
%   A warning message is printed if also the second singular 
%   value is grater than EPP.
%   
%   Inputs:  G   - system matrix.
%            P   - vector of poles.
%            epp - tolerance (optional) for singular values, 
%                  default value is 1/eps which most often is
%                  to large typical value around 1e14.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to pole P(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to pole P(i).
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No direction with gain larger 
%                               than EPP.
%                     err = -1, One basis vector is not enough
%                               to describe the subspace with 
%                               gains larger than EPP.
%
%   Written by: Kjetil Havre, 24/11-1995, e-mail: kjetil@ife.no
%
%   See also: IZDE, OZDE, ZDSVD and PDFST.
%
%   NOTE! This function should not be used for anything else than 
%         checking the results returned by PDFST, since it depends 
%         poor numerical properties. 
%         In the function PDFST the pole directions are calculated 
%         avoiding these poor numerical properties.
%
   function [I, O, err] = pdsvd(G, P, epp)

   I = []; O = []; err = 0;

   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = pdsvd(G, P, epp)' );
      return
   end
   
   if nargin < 2
      disp( 'Usage: [I,O] = pdsvd(G, P, epp)' );
      return
   end
   
   if nargin < 3
      epp = 1/eps;
   end
   
   for i=1:max(size(P))
      wp  = P(i)/sqrt(-1);
      Gp = frsp(G, wp);              % This may cause problems.
      [u,s,v] = svd(Gp(1:ny,1:nu));
      nr = min(size(s));
      if s(1, 1) > epp
         if nr > 1
            if s(2, 2) > epp
               disp(['Warning: Second singular value is larger than eep: ', num2str(epp),','])
               disp(['         S(2,2): ', num2str(s(nr-1,nr-1)),'.'])
	       err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      else
         disp(['Warning: Gain of G(p) is not larger than EPP: ', num2str(epp)]);
	 disp(['         for i=',int2str(i),', p=',num2str(P(i)),'.'])
         I(:,i) = zeros(nu,1);
	 O(:,i) = zeros(ny,1);
         err = 1;
      end
   end
end

   
