Symmetric Matrix Eigenvalue Problems - Numerical Methods in Engineering with MATLAB - page 364

utilized to store the vectors u that are neededinthe computation of the transformation
matrix P .

functionA=householder(A)
%HousholderreductionofAtotridiagonalformA=[c\d\c].
%Extractcanddbyd=diag(A),c=diag(A,1).
%USAGE:A=householder(A)

n = size(A,1);
fork=1:n-2
u = A(k+1:n,k);
uMag = sqrt(dot(u,u));
ifu(1)<0;uMag=-uMag;end
u(1) = u(1) + uMag;
A(k+1:n,k) = u; % Save u in lower part of A.

H = dot(u,u)/2;
v = A(k+1:n,k+1:n)*u/H;
g = dot(u,v)/(2*H);
v=v-g*u;
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - v*u' - u*v';
A(k,k+1) = -uMag;

end

householderP

The function householderP returns the accumulated transformationmatrix P . There
is no need to call it ifonly the eigenvalues aretobecomputed. Note that the input
parameter A is not the original matrix, but the matrix returnedby householder .

functionP=householderP(A)
% Computes transformation matrix P after
% householder reduction of A is carried out.
%USAGE:P=householderP(A).

n = size(A,1);
P = eye(n);
fork=1:n-2
u = A(k+1:n,k);
H = dot(u,u)/2;
v = P(1:n,k+1:n)*u/H;