%
%  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;
   test = input('1 for match inf or 2 for match 0')
%
%  Now construct the reduced basis V  from k-step Arnoldi
%

%
%  Use Arnoldi with A  to match moments at infty
%

   if (test == 1),
      [V,H,f] = ArnoldiC(A,k,B);
      Ak = H;
   end
%
%  Use Arnoldi with inv(A)  to match moments at 0 
%
   
   if (test == 2),
      [V,H,f] = ArnoldiC(inv(A),k,B);
      Ak = V'*A*V;
   end
%  NOTE:   You would NEVER compute inv(A) 
%          in practice.    Instead, you would factor
%          [L,U] = lu(A) and solve triangular
%          systems to get effect of inv(A)
%

%
%  Construct the reduced order model
%
   
   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)