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}]:>RandomChoice[{-1,1}]RandomReal[{6,23}]
},{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.10425
A=FrancisIteration[AOriginal,ρ];
{Diagonal[A],Diagonal[A,-1]}

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

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

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

RayleighShift[A_]:= Module[
{},
A[[-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]
]

x=SparseArray[{1}->1.0,{Length[A]}];
Normal[x]

FrancisIterationMultiple[AIn_,ρList_]:= Module[
{A=AIn,x,H},
x=SparseArray[{1}->1.0,{Length[AIn]}];
Do[ x=AIn.x;x[[1]]-=ρList[[i]],{i,Length[ρList]}];
H=Householder[x,1];
A=H.A.Transpose[H];
AOut=A;
Do[
x=A[[All,i]];
H=Householder[x,i+1];
A=Chop[H.A.Transpose[H]],
{i,1,n-2}];
A
]

Eigenvalues[A[[-2;;-1,-2;;-1]]]

A=AOriginal;
Table[
{ρ1,ρ2}=Eigenvalues[A[[-2;;-1,-2;;-1]]];
A=FrancisIterationMultiple[A,{ρ1,ρ2}];
MatrixPlot[A,
PlotLabel->MatrixForm[A[[-2;;-1,-2;;-1]]]
],
{12}]

Timing Comparisons

The following timings are quite different!

n=256;
A=RandomReal[{-1,1},{n,n}];
ASym=A + Transpose[A];
Timing[Eigenvalues[A];]
Timing[Eigenvalues[ASym];]
Timing[Eigensystem[A];]
Timing[Eigensystem[ASym];]
Timing[SingularValueList[A];]
Timing[SingularValueList[ASym];]
Timing[SingularValueDecomposition[A];]
Timing[SingularValueDecomposition[ASym];]

Lets see how these change with problem size.

Data=Table[
n=Floor[50 (1.5^p)];
A=RandomReal[{-1,1},{n,n}];
ASym=A + Transpose[A];
{
{n,Timing[Eigenvalues[A];][[1]]},
{n,Timing[Eigenvalues[ASym];][[1]]},
{n,Timing[Eigensystem[A];][[1]]},
{n,Timing[Eigensystem[ASym];][[1]]},
{n,Timing[SingularValueList[A];][[1]]},
{n,Timing[SingularValueList[ASym];][[1]]},
{n,Timing[SingularValueDecomposition[A];][[1]]},
{n,Timing[SingularValueDecomposition[ASym];][[1]]}
},
{p,1,5}
];

ListLogLogPlot[
Table[Data[[All,i]],{i, 8}],
Joined->True,
Mesh->All,
Frame->True,
GridLines->Automatic,
PlotStyle->Table[Hue[i/8],{i,0,8}],
FrameLabel->{"n","Sec"}
]

Looks as though the difference is in the constant not the power.

Data[[-1]]