(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



HouseU[Col_]:=Module[{x,\[Gamma],\[Tau]=Max[Abs[Col]]},
(* Watkins 3rd 3.2.35*)
(* returns \[Gamma] and u *)
If[\[Tau]==0.0,Return[{0,Col}]];
x=Col/\[Tau];
\[Tau]=Sign[x[[1]]]Norm[x];
x[[1]]+=\[Tau];
\[Gamma]=x[[1]]/\[Tau];
x=x/x[[1]];
{\[Gamma],x}]


HessenbergReducedSymmetric[n_]:=Module[{A=RandomReal[{-1,1},{n,n}],Id=IdentityMatrix[n],\[Gamma],u,i},A=A+Transpose[A];
Do[
{\[Gamma],u}=HouseU[A[[i+1;;-1,i]]];
U=Join[ConstantArray[0,i],u];
Q=Id-Outer[Times,\[Gamma] U,U];
A=Chop[Q.A.Transpose[Q]],{i,1,n-2}];
A]


HessenbergDiagonalRatio[A_]:=Module[{},
Abs[Diagonal[A,-1]/Drop[Diagonal[A,0],-1]]
]
CreateBulge[AIn_,\[Rho]_]:=Module[
{A=AIn,i=1,\[Gamma],v,u,U,Q,Id =IdentityMatrix[Length[AIn]]},
v=A[[i;;-1,i]];
v[[1]]-=\[Rho];
{\[Gamma],u}=HouseU[v];
U=Join[ConstantArray[0,i-1],u];
Q=Id-Outer[Times,\[Gamma] U,U];
A=Q.A.Transpose[Q];
A]
ChaseBulge[AIn_]:=Module[
{A=AIn,\[Gamma],u,U,Q,j,Id =IdentityMatrix[Length[AIn]],i=1},
Do[
{\[Gamma],u}=HouseU[A[[i+j;;i+j+1,j]]];
U=Join[ConstantArray[0,j],u,ConstantArray[0,n-j-2]];
Q=Id-\[Gamma] Outer[Times,U,U];
A=Chop[Q.A.Transpose[Q]],
{j,1,n-2}
];
A
]
FrancisIteration[AIn_,\[Rho]_]:=Module[
{A=AIn},
A=CreateBulge[A,\[Rho]];
ChaseBulge[A]
]


WilkinsonShift[ASub_]:=Sort[Eigenvalues[ASub],Abs[#1-A[[-1,-1]]]<Abs[#2-A[[-1,-1]]]&][[1]]
WilkinsonShift[{{2,1},{1,1.8}}]



