%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	Matlab program:: Tests the results of SDC.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Read in the data file called "sdebug.m"
l=0;w=0;a=0;wi=0;wr=0;m=0;v1=0;aeig=0;n=0;d0=0;errest=0;
sdebug
str=sprintf('Done reading sdebug.m'); disp(str);

% reshape 
n=m;
wr = l(:,1); wi = l(:,2);
a = reshape(a,n,n)'; 
i = sqrt(-1);
[v1,d0] = eig(a);
for j = 1:m
   aeig(j) = d0(j,j);
   w(j) = wr(j) + wi(j)*i;
end
if (uppert)
   t = reshape(t,n,n)'; 
end

% Print error estimate
str=sprintf('Overall error estimate: %g', errest); disp(str);

% make sure t is upper triangular
if (uppert)
   fail = 0;
   for j = 1:m-2
      for k = j+2:m
         if (abs(t(k,j)) > 1.0e-10)
            tstr=sprintf('%s\n(%d, %d)', tstr,k,j);
            fail = 1;
         end
      end
   end
   if (fail == 1)
      str2=sprintf('T is not upper triangular! See tstr.'); 
      disp(str2);
   else
      str2=sprintf('T is upper triangular.'); 
      disp(str2);
   end
end

% use matlab solver and compare eigenvalues
maxdiff = 1.0e-15;
maxdiff = 0.0e0;
for j = 1:m
   mindiff = 1000.0;
   for k = 1:m
      diff = abs( aeig(k) - w(j) );
      if (diff < mindiff)
         mindiff = diff;
      end
   end
   if (mindiff > maxdiff)
      if (mindiff > 1.0e-8)
         str=sprintf('SDC bad eigenvalue : %g (%d)', w(j), j); disp(str);
      end
      maxdiff = mindiff;
   end
end
str=sprintf('Max difference in eigenvalues vs matlab = %g', maxdiff); disp(str);

% Compare eigenvalues of A and T
if (uppert)
   v1=0;d0=0;teig=0;
   i = sqrt(-1);
   [v1,d0] = eig(t);
   for j = 1:m
      teig(j) = d0(j,j);
   end
   maxdiff = 1.0e-15;
   for j = 1:m
      mindiff = 1000.0;
      for j = 1:m
         diff = abs( aeig(j) - teig(j) );
         if (diff < mindiff)
            mindiff = diff;
         end
      end
      if (mindiff > maxdiff)
         maxdiff = mindiff;
      end
   end
   str=sprintf('Max difference in eigenvalues of A and T = %g', maxdiff); 
   disp(str);
end

if (svects | evects)
   q = reshape(q,n,n)'; 
   if (evects)
      ov = reshape(v,n,n)';
      nv = zeros(n);
      j = 1;
      while (j <= n)
	 if (wi(j) == 0)
	    nv(:,j) = ov(:,j);
	    ei(j) = wr(j);
	    j=j+1;
	 else
	    nv(:,j) = ov(:,j) + i*ov(:,j+1);
	    nv(:,j+1) = ov(:,j) - i*ov(:,j+1);
	    ei(j) = wr(j) + i*wi(j); ei(j+1) = wr(j+1) + i*wi(j+1);
	    j=j+2;
	 end
      end
      v=nv;
   end
end

%%% Some simple tests: %%%

% is Q orthogonal?

if (svects)
   str=sprintf('|| q ||2 = %g \n|| q^t ||2 = %g',norm(q), norm(q'));
   disp(str);
   str=sprintf('|| q^t * q || = %g \n|| q * q^t || = %g',norm(q'*q), norm(q*q'));
   disp(str);
end

% Is q'*a*q = t ? 

if (uppert & svects)
   str=sprintf('|| q^t * a * q - t || = %g', norm(q'*a*q - t)); disp(str);
   str=sprintf('|| a || - || t || = %g', norm(a) - norm(t)); disp(str);
end


% A*v = lambda*v ?
if (evects)
   for j = 1:m
      nrm(j) = norm(a*v(:,j) - ei(j)*v(:,j));
   end
   str=sprintf('Maximum of || A*v - lambda*v || =  %g', max(nrm));disp(str);
end
