Page 268 - MATLAB an introduction with applications

Numerical Methods ——— 253

 72 3 –1
 28 5 1 
Example E4.28: Transform the matrix [A]=   35 12 9  into tridiagonal form using Householder
 
  –1 1 9 7 
reduction.
Solution:
The following program is used.
%Input – A is an nxn symmetric matrix
A = [7 2 3 –1;2 8 5 1;3 5 12 9;–1 1 9 7];

%Output – T is a tridiagonal matrix
[n,n] = size(A);
for k = 1:n–2
s = norm(A(k+1:n,k));
if (A(k+1,k)<0)
s = –s;
end
r = sqrt(2*s*(A(k+1,k)+s));
W(1:k) = zeros(1,k);
W(k+1) = (A(k+1,k)+s)/r;
W(k+2:n) = A(k+2:n,k)’/r;
V(1:k) = zeros(1,k);
V(k+1:n) = A(k+1:n,k+1:n)*W(k+1:n)’;
c = W(k+1:n)*V(k+1:n)’;
Q(1:k) = zeros(1,k);
Q(k+1:n) = V(k+1:n)–c*W(k+1:n);
A(k+2:n,k) = zeros(n–k–1,1);
A(k,k+2:n) = zeros(1,n–k–1);
A(k+1,k) = –s;
A(k,k+1) = –s;
A(k+1:n,k+1:n) = A(k+1:n,k+1:n)–2*W(k+1:n)’*Q(k+1:n)–2*Q(k+1:n)’*W(k+1:n);
end
T = A;
fprintf(‘Matrix in tridiagonal form is\n’);
disp(T)

The output is as follows:
The matrix in tridiagonal form is
7.0000 –3.7417 0 0
–3.7417 10.6429 9.1309 0
0 9.1309 10.5942 4.7716
0 0 4.7716 5.7629