Assignment 6 Scripts

P1

% P1
% Tests Schur4
A11 = [3 16; -1 3]    
A22 = [1  1; -4 1]
A12 = [5  6; 7 8 ];
clc
for mu = [0 .9999  .99999999  1]
   T = [ A11  A12 ; zeros(2,2) A22+mu*(A11-A22)]
   [U,sep] = Schur4(T);
   S = U'*T*U
   disp(' ')
   disp(sprintf('sep = %10.3e',sep))
   evT11 = eig(T(1:2,1:2)); 
   evT22 = eig(T(3:4,3:4));
   evS11 = eig(S(1:2,1:2));
   evS22 = eig(S(3:4,3:4));
   EigenValuesT11_and_S22 = [evT11 evS22]
   EigenValuesT22_and_S11 = [evT22 evS11]
end

P2

% P2: Tests BlockHess

randn('seed',0)
n=12;
A = randn(n,n);
clc
for r=1:4
    [U1,S] = qr(randn(n,r),0);
    [U,H] = BlockHess(A,U1);
    errU = norm(U'*U - eye(n,n));
    errH = norm(U'*A*U-H);
    errU1 = norm(U1-U(:,1:r));
    disp(' ')
    disp(sprintf('r=%1d       errU = %8.3e   errH = %8.3e   errU1 = %8.3e',r,errU,errH,errU1))
    disp(' ')
    disp('H = ')
    for i=1:n
        disp(sprintf('  %5.2f',H(i,:)))
    end
    disp('---------------------------------------------------------------------------------')
end

GTD

  function X = SylvesterSolver(F,G,H)
% F is p-by-p upper triangular
% G is q-by-q upper triangular
% F and G have no eigenvalues in common
% H is p-by-q
% X is p-by-q so that FX-XG = H

[p,q] = size(H);
X = zeros(p,q);
I = eye(p,p);
for k=1:q
   X(:,k) = (F - G(k,k)*I)\(H(:,k) + X(:,1:k-1)*G(1:k-1,k));
end



function kappa = CondInf(T)
% T is a nonsingular n--by-n upper triangular matrix
% kappa is an estimate of its inf-norm condition

[n,n] = size(T);
d = zeros(n,1);
y = zeros(n,1);
d(n) = 1;
y(n) = d(n)/T(n,n);
s = T(1:n-1,n)*y(n);
for k=n-1:-1:1
    % Choose d(k) = 1 or -1. The goal is to maximize growth.
    % Note that
    %     y(k)     = (d(k) - s(k))/T(k,k)
    %     s(1:k-1) = s(1:k-1) + T(1:k-1,k)*y(k)
    
    ykplus = (1 - s(k))/T(k,k);
    splus = s(1:k-1) + T(1:k-1,k)*ykplus;
    ykminus = (-1 - s(k))/T(k,k);
    sminus = s(1:k-1) + T(1:k-1,k)*ykminus;
    if abs(ykplus) + norm(splus,1) > abs(ykminus) + norm(sminus,1)
        y(k) = ykplus;
        s(1:k-1) = splus;
    else
        y(k) = ykminus;
        s(1:k-1) = sminus;
    end
end
kappa = norm(T,inf)*norm(y,inf);