A5 Scripts

P1

% P1
clc
randn('seed',0)

% n = 4 tests

disp('n = 4 examples')
disp(' ')
disp('norm(Q''*Q - I)   norm(Q''*A*Q - D)    norm(Q - I,''fro'')  ')
disp('------------------------------------------------------------')
I = eye(4,4);
for k=1:10
   A = randn(4,4); A = A-A';
   [Q,D] = skew34(A);
   eQ = norm(Q'*Q - I);
   eD = norm(Q'*A*Q - D);
   eI = norm(Q-I,'fro');
   disp(sprintf('  %10.3e        %10.3e           %10.3e',eQ,eD,eI))
end

disp(' ')
disp('n = 3 examples')
disp(' ')
disp('norm(Q''*Q - I)   norm(Q''*A*Q - D)    norm(Q - I,''fro'')  ')
disp('------------------------------------------------------------')
I = eye(3,3);
for k=1:10
   A = randn(3,3); A = A-A';
   [Q,D] = skew34(A);
   eQ = norm(Q'*Q - I);
   eD = norm(Q'*A*Q - D);
   eI = norm(Q-I,'fro');
   disp(sprintf('  %10.3e        %10.3e           %10.3e',eQ,eD,eI))
end

P2

% P2
% Tests SkewJacobi

randn('seed',0)
clc
tol = .000000000000001;
itmax = 25;

for n = 20:21
   
   A = randn(n,n); A = A - A'; A = A/norm(A,'fro');

   [Q,D,SweepHistory] = SkewJacobi(A,'skew34',tol,itmax);
   disp(' ')
   disp(sprintf('n = %2d  with skew34 ...',n))
   disp(' ')
   disp(' Sweep      BlockOff(A)')
   disp('--------------------------')
   for k = 1:length(SweepHistory)
      disp(sprintf('   %2d       % 10.3e',k,SweepHistory(k)))
   end
   disp(' ')
   disp(sprintf('norm(Q''*Q - I)   = %10.3e',norm(Q'*Q-eye(n,n))))
   disp(sprintf('norm(Q''*A*Q - D) = %10.3e',norm(Q'*A*Q - D)))

   [Q,D,SweepHistory] = SkewJacobi(A,'schur',tol,itmax);
   disp(' ')
   disp(sprintf('n = %2d  with Schur ...',n))
   disp(' ')
   disp(' Sweep      BlockOff(A)')
   disp('--------------------------')
   for k = 1:length(SweepHistory)
      disp(sprintf('   %2d       % 10.3e',k,SweepHistory(k)))
   end
   disp(' ')
   disp(sprintf('norm(Q''*Q - I)   = %10.3e',norm(Q'*Q-eye(n,n))))
   disp(sprintf('norm(Q''*A*Q - D) = %10.3e',norm(Q'*A*Q - D)))
   
end

P3

% P3
% Tests Special Schur

% Correctness

clc
randn('seed',0)
n = 10;
p = 3;
X = randn(n,p); B = randn(p,p); B = B+B';
[Q,D] = SpecialSchur(X,B);
[Q1,D1] = schur(eye(n,n)+X*B*X');
d = sort(diag(D));
d1 = sort(diag(D1));
error = norm(Q'*(eye(n,n)+X*B*X')*Q - D);
disp(sprintf('norm(Q''(I+XBX'')Q - D) = %10.3e',error))
errorQ = norm(Q'*Q - eye(n,n));
disp(sprintf('norm(Q''Q - eye(n,n) ) = %10.3e',errorQ))
disp(' ')
disp(' k           D(k,k)                D1(k,k)')
disp('-----------------------------------------------------')
for k=1:n
   disp(sprintf('%2.0d   %20.16f   %20.16f',k,d(k),d1(k)))
end

% Efficiency

nRepeat = 10; n = 300; p = 5;
X = randn(n,p);
B = randn(p,p); B = B+B';
tic 
for k=1:nRepeat
   [Q,D] = SpecialSchur(X,B);
end
t1 = toc;

tic 
for k=1:nRepeat
   [Q,D] = schur(eye(n,n)+X*B*X');
end
t2 = toc;

disp(' ')
disp(sprintf('n = %2.0d,  p = %2.0d',n,p))
disp(sprintf('Schur Time / SpecialSchur Time = %6.2f',t2/t1))
   

P4

% P4
% Tests GenEval

clc
randn('seed',0)
n = 10;

% Random ill-conditioned A and B ...

[Z,T] = qr(randn(n,n)); A = Z*diag((0:n-1)+.000001)*Z';
[Z,T] = qr(randn(n,n)); B = Z*diag((0:n-1)+.000001)*Z';

[X,lambda] = GenEval(A,B);

disp('  j            lambda(j)           norm( (A - lambda(j)*B)*X(:j))/norm(X(:,j))')
disp('------------------------------------------------------------------------------')
for j=1:n
   error = norm((A - lambda(j)*B)*X(:,j))/norm(X(:,j));
   disp(sprintf(' %2d  %27.15f                  %10.3e',j,lambda(j),error))
end