%
%  This is a test driver for MOR via 
%  Shift-Invert Arnoldi.   It is applied
%  to a CD-Player Model of order n = 120
%  It should produce a reduced model of
%  order k = 20.  Bode plots comparing the
%  output of the reduced vs full system are
%  produced.   Also, eigenvalue distribution
%  of the reduced model is plotted.
%

%clear;

%
%  This will input matrices A, B, C for CD-Player
%
   load CD.mat;
%
%  Set k = requested order of the reduced model
%
   k = 20;
% 
%  A positive real number s1 is input and
%  
%    shift = s1*i  = s1*sqrt(-1)   is constructed 
%                                  and displayed
%
   s1 = input('Shift sigma')

   shift = s1*i

   disp(' strike a key ')
pause

   n = size(A,1); I = eye(n);

%
%  Run the Arnoldi process for k-steps  to
%  construct reduced basis V
%

   [V,H,f] = ArnoldiC(inv(A - s1*I),k,B);
%
%  Construct the k-th order reduced model  
%
%  NOTE:   You would NEVER compute inv(A - s1*I) 
%          in practice.    Instead, you would factor
%          [L,U] = lu(A - s1*I ) and solve triangular
%          systems to get effect of inv(A - s1*I)
%

   Ak = V'*A*V;
   Bk = V' * B;
   Ck = C* V;

 sys = ss(A,B,C,0);
 sys_r = ss(Ak,Bk,Ck,0);

%
%  Display Bode plots from Matlab sigma function
%
   figure(1)
   clf
   sigma(sys)
   hold
   sigma(sys_r,'r')
   hold
   legend('full system','reduced system')
   title('CD Player')
   
   n = size(A,1);
   k =  size(Ak,1);
%
%  Display Bode plots from our own routine 
%
   
   figure(2)
   clf
   e1 = 0; e2 = 7;
   textstring = [ 'Freq-Response CD-Player :  n = ', num2str(n),  ' k = ', num2str(k)];
   Mysigma_log(A,B,C,D,e1,e2,textstring,'b');
   hold
   Mysigma_log(Ak,Bk,Ck,D,e1,e2,textstring,'r-.');
   legend('Original', 'Reduced')
   
   Dimension_Full_System = size(A,1) 
   Dimension_Reduced_System =  size(Ak,1)
%
%  Display eigenvalue plots for original and reduced systems
%

   t = eig(A); s = eig(Ak);
   figure(3)
   plot(real(t), imag(t), '*' ,'linewidth',2.5);
   hold on
   
   plot(real(s), imag(s), 'ro','linewidth',2.5);
   
   hold off