Householder[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=SparseArray[Band[{1,1}]->1,{n,n}];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-(2/((Norm[u])^2))Outer[Times,u,u]
]

Francis Algorithm

OK lets pretend that we have reduced a symmetric matrix to Hessenberg form.  This is a symmetric tridiagonal matrix.

n=9;
AOriginal=Normal[SparseArray[{
Band[{1,1}]:>RandomReal[{-2,2}],
Band[{1,2}]:>RandomReal[{-2,2}]
},{n,n}]];
AOriginal=AOriginal+Transpose[AOriginal];
MatrixPlot[AOriginal]

We want to implement the Francis algorithm with a shift ρ.  Here is Wednesdays manual version.

A=AOriginal;
Id =IdentityMatrix[n];
e1=Table[KroneckerDelta[1,i],{i,1,n}];
ρ=0.1;
B=A-ρ Id;
x=B[[All,1]];
H=Householder[x,1];
(*MatrixPlot[H.A]*)
A=H.A.Transpose[H];
MatrixPlot[A]

Pics=Table[
x=A[[All,i]];
H=Householder[x,i+1];
A=Chop[H.A.Transpose[H]];
MatrixPlot[A],
{i,1,n-2}];

Manipulate[
Pics[[i]],
{i,1,n-2,1}
]

FrancisIteration[AIn_,ρ_]:= Module[
{A=AIn,x,H},
x=AIn[[All,1]];
x[[1]]-=ρ;
H=Householder[x,1];
A=H.A.Transpose[H];
Do[
x=A[[All,i]];
H=Householder[x,i+1];
A=Chop[H.A.Transpose[H]],
{i,1,n-2}];
A
]

Lets make sure that this does not change the eigenvalues.

ρ=12.3;
A=FrancisIteration[AOriginal,ρ];
Eigenvalues[A]
Norm[Eigenvalues[AOriginal]-Eigenvalues[A]]

Did it make the matrix more triangular?

{Diagonal[AOriginal],Diagonal[AOriginal,-1]}
{Diagonal[A],Diagonal[A,-1]}

But it does if we pick ρ well.  The following is called a Rayleigh shift.

ρ=1.0452
A=FrancisIteration[AOriginal,ρ];
{Diagonal[A],Diagonal[A,-1]}

ρ=0.85934
A=FrancisIteration[A,ρ];
{Diagonal[A],Diagonal[A,-1]}

ρ=0.850418
A=FrancisIteration[A,ρ];
{Diagonal[A],Diagonal[A,-1]}

ρ=0.850416
A=FrancisIteration[A,ρ];
{Diagonal[A],Diagonal[A,-1]}

RayleighShift[A_]:= Module[
{},
AEnd[[-1,-1]]
]

A Wilkinson shift uses the eigenvalues of the bottom 2x2 matrix.

WilkinsonShift[A_]:= Module[
{AEnd=A[[-2,-2]],λ1,λ2},
{λ1,λ2}=Eigenvalues[AEnd];
If[ Abs[λ1-AEnd[[-1,-1]]]<Abs[λ2-AEnd[[-1,-1]]], λ1,λ2]
]