A5 Scripts

P1

%P1 Tests HouseC
clc
randn('seed',0)
i = sqrt(-1);

% Simple Small "Easy" Example
disp('Example 1:')
format long
n = 5;
x = randn(n,1) + i*randn(n,1);
v = HouseC(x)
disp(sprintf('Norm(v,2) = %20.15f',norm(v)))
disp(' ')
disp('(I - 2*v*v'')x =')
disp(' ')
disp( x - (2*(v'*x))*v)

disp(' ')

% Simple Small "Tricky" Example
disp('Example 2:')
n = 5; x = randn(n,1) + i*randn(n,1); x(2:n) = .00000001*x(2:n);
v = HouseC(x)
disp(sprintf('Norm(v,2) = %20.15f',norm(v)))
disp(' ')
disp('(I - 2*v*v'')x =')
disp(' ')
disp( x - (2*(v'*x))*v)

disp(' ')
% Simple Small Real
disp('Example 3:')
n = 3; x = randn(n,1);
v = HouseC(x)
disp(sprintf('Norm(v,2) = %20.15f',norm(v)))
disp(' ')
disp('(I - 2*v*v'')x =')
disp(' ')
disp( x - (2*(v'*x))*v)

disp(' ')


% Simple Small Pure Imaginary Example
disp('Example 4:')
n = 3; x = i*randn(n,1);
v = HouseC(x)
disp(sprintf('Norm(v,2) = %20.15f',norm(v)))
disp(' ')
disp('(I - 2*v*v'')x =')
disp(' ')
disp( x - (2*(v'*x))*v)

P2

% P2: Tests TriDiag

clc
randn('seed',0);

% Example 1

disp('A real problem...')
disp(' ')
n = 5;
A = randn(n,n);
A = A+A';
T = TriDiag(A)
dA = sort(eig(A)); 
dT = sort(eig(real(T)));
disp('     Eigenvalues of A       Eigenvalues of T')
disp('----------------------------------------------------')
for k=1:n
    disp(sprintf('  %20.15f   %20.15f',dA(k),dT(k)))
end

% Example 2
disp(' ')
disp('--------------------------')
disp(' ')
disp('A tridiagonal problem...')
disp(' ')
n = 5;
A = randn(n,n) + sqrt(-1)*randn(n,n);
A = tril(triu(A+A',-1),1);
T = TriDiag(A)
dA = sort(eig(A)); 
dT = sort(eig(real(T)));
disp('     Eigenvalues of A       Eigenvalues of T')
disp('----------------------------------------------------')
for k=1:n
    disp(sprintf('  %20.15f   %20.15f',dA(k),dT(k)))
end

% Example 3
disp(' ')
disp('--------------------------')
disp(' ')
disp('A general problem...')
disp(' ')
n = 5;
A = randn(n,n) + sqrt(-1)*randn(n,n);
A = A+A';
T = TriDiag(A)
dA = sort(eig(A)); 
dT = sort(eig(real(T)));
disp('     Eigenvalues of A       Eigenvalues of T')
disp('----------------------------------------------------')
for k=1:n
    disp(sprintf('  %20.15f   %20.15f',dA(k),dT(k)))
end

% Example 4
disp(' ')
disp('--------------------------')
disp(' ')
disp('A low rank problem...')

n = 5;
X = randn(n,2) + sqrt(-1)*randn(n,2);
A = X*X';
T = TriDiag(A)
dA = sort(eig(A)); 
dT = sort(eig(real(T)));
disp('     Eigenvalues of A       Eigenvalues of T')
disp('----------------------------------------------------')
for k=1:n
    disp(sprintf('  %20.15f   %20.15f',dA(k),dT(k)))
end

P3

% P3: Tests BorderDiag

clc
randn('seed',0)

% Example 1
disp('Example 1. (Nice)');
n = 5;
d = -sort(-randn(n-1,1));
v = randn(n-1,1);
alpha = randn(1,1);
A = [diag(d) v; v' alpha ]
[Q_True,lambda_True] = schur(A);
lambdaVec = BorderedDiag(d,v,alpha);
disp(' ')
disp('      Eigenvalues           Eigenvalues')
disp('       Via Schur          via BorderedDiag')
disp('-----------------------------------------------')
e1 = sort(diag(lambda_True));
e2 = sort(lambdaVec);
for k=1:n
    disp(sprintf('%20.16f    %20.16f',e1(k),e2(k)))
end
disp(' ')

% Example 2
disp('Example 2. (|d(2)-d(3)| = .000001)');
n = 5;
d = -sort(-randn(n-1,1));
d(3) = d(2) - .000001;
v = randn(n-1,1);
alpha = randn(1,1);
A = [diag(d) v; v' alpha ]
[Q_True,lambda_True] = schur(A);
lambdaVec = BorderedDiag(d,v,alpha);
disp(' ')
disp('      Eigenvalues           Eigenvalues')
disp('       Via Schur          via BorderedDiag')
disp('-----------------------------------------------')
e1 = sort(diag(lambda_True));
e2 = sort(lambdaVec);
for k=1:n
    disp(sprintf('%20.16f    %20.16f',e1(k),e2(k)))
end
disp(' ')


% Example 3
disp('Example 3. (v(2) = .000001)');
n = 5;
d = -sort(-randn(n-1,1));
v = randn(n-1,1);
v(2) = .000001;
alpha = randn(1,1);
A = [diag(d) v; v' alpha ]
[Q_True,lambda_True] = schur(A);
lambdaVec = BorderedDiag(d,v,alpha);
disp(' ')
disp('      Eigenvalues           Eigenvalues')
disp('       Via Schur          via BorderedDiag')
disp('-----------------------------------------------')
e1 = sort(diag(lambda_True));
e2 = sort(lambdaVec);
for k=1:n
    disp(sprintf('%20.16f    %20.16f',e1(k),e2(k)))
end

P4

% P4: Tests BorderDiagQ

clc
randn('seed',0)

% Example 1
disp('Example 1. (Nice)');
n = 5;
d = -sort(-randn(n-1,1));
v = randn(n-1,1);
alpha = randn(1,1);
A = [diag(d) v; v' alpha ]
[Q_True,lambda_True] = schur(A);
disp('Q via Schur...')
disp(' ')
disp(Q_True)
disp('Eigenvalues of A...')
disp(' ')
disp(diag(lambda_True)')
disp(' ')
disp('Q via BorderedDiagQ ...')
Q = BorderedDiagQ(d,v,alpha,diag(lambda_True));
disp(' ')
disp(Q)
disp(sprintf('norm(Q''*Q-I)            = %10.3e',norm(Q'*Q-eye(n,n))))
disp(sprintf('norm(A*Q-Q*lambda_True) = %10.3e',norm(A*Q-Q*lambda_True)))
disp(' ')
disp('-------------------------')
disp(' ')

% Example 2
disp('Example 2. (|d(2)-d(3)| = .000001)');
n = 5;
d = -sort(-randn(n-1,1));
d(3) = d(2) - .000001;
v = randn(n-1,1);
alpha = randn(1,1);
A = [diag(d) v; v' alpha ]
[Q_True,lambda_True] = schur(A);
disp('Q via Schur...')
disp(' ')
disp(Q_True)
disp('Eigenvalues of A...')
disp(' ')
disp(diag(lambda_True)')
disp(' ')
disp('Q via BorderedDiagQ ...')
Q = BorderedDiagQ(d,v,alpha,diag(lambda_True));
disp(' ')
disp(Q)
disp(sprintf('norm(Q''*Q-I)            = %10.3e',norm(Q'*Q-eye(n,n))))
disp(sprintf('norm(A*Q-Q*lambda_True) = %10.3e',norm(A*Q-Q*lambda_True)))
disp(' ')
disp('-------------------------')
disp(' ')

% Example 3
disp('Example 3. (v(2) = .000001)');
n = 5;
d = -sort(-randn(n-1,1));
v = randn(n-1,1);
v(2) = .000001;
alpha = randn(1,1);
A = [diag(d) v; v' alpha ]
[Q_True,lambda_True] = schur(A);
disp('Q via Schur...')
disp(' ')
disp(Q_True)
disp('Eigenvalues of A...')
disp(' ')
disp(diag(lambda_True)')
disp(' ')
disp('Q via BorderedDiagQ ...')
Q = BorderedDiagQ(d,v,alpha,diag(lambda_True));
disp(' ')
disp(Q)
disp(sprintf('norm(Q''*Q-I)            = %10.3e',norm(Q'*Q-eye(n,n))))
disp(sprintf('norm(A*Q-Q*lambda_True) = %10.3e',norm(A*Q-Q*lambda_True)))
disp(' ')
disp('-------------------------')
disp(' ')