Jacobi Creators: General

J[{i_,j_},{p1_,p2_},m_]:= Module[
{n=Norm[{p1,p2}]},
{pp1,pp2}={p1,p2}/n;
(* (({{ p1, p2}, {-p2, P1}}))*)
 SparseArray[{
{i,i}->pp1,
{j,j}->pp1,
{i,j}->pp2,
{j,i}->-pp2,
Band[{1,1}]->1
},
{m,m}
]
]

J[{i_,j_,k_,l_},{p1_,p2_,p3_,p4_},m_]:= Module[
{n=Norm[{p1,p2,p3,p4}],pp1,pp2,pp3,pp4},
(* LinAlg+App 332(2001) 37-80 Fasbender p65 *)
(* (({{ p1,  p2,  p3,  p4}, {-p2,  p1,  p4, -p3}, {-p3, -p4,  p1,  p2}, {-p4,  p3, -p2,  p1}})) *)
{pp1,pp2,pp3,pp4}={p1,p2,p3,p4}/n;
 SparseArray[{
{i,i}->pp1,  {i,j}->pp2,{i,k}->pp3, {i,l}->pp4,
{j,i}->-pp2,{j,j}->pp1,{j,k}->pp4,{j,l}->-pp3,
{k,i}->-pp3,{k,j}->-pp4,{k,k}->pp1,{k,l}->pp2,
{l,i}->-pp4,{l,j}->pp3,{l,k}->-pp2,{l,l}->pp1,
Band[{1,1}]->1
},
{m,m}
]
]

J[
{i1_,i2_,i3_,i4_,i5_,i6_,i7_,i8_},
{p1_,p2_,p3_,p4_,p5_,p6_,p7_,p8_},
m_]:= Module[
{n=Norm[{p1,p2,p3,p4,p5,p6,p7,p8}],
pp1,pp2,pp3,pp4,pp5,pp6,pp7,pp8},
(* LinAlg+App 332(2001) 37-80 Fasbender p66 *)
(* (({{ p1,  p2,  p3,  p4,  p5,  p6,  p7,  p8}, {-p2,  p1,  p4, -p3,  p6, -p5, -p8,  p7}, {-p3, -p4,  p1,  p2,  p7,  p8, -p5, -p6}, {-p4,  p3, -p2,  p1,  p8, -p7,  p6, -p5}, {-p5, -p6, -p7, -p8,  p1,  p2,  p3,  p4}, {-p6,  p5, -p8,  p7, -p2,  p1,  p4, -p3}, {-p7,  p8,  p5, -p6, -p3, -p4,  p1,  p2}, {-p8, -p7,  p6, -p5, -p4,  p3, -p2,  p1}})) *)
{pp1,pp2,pp3,pp4,pp5,pp6,pp7,pp8}={p1,p2,p3,p4,p5,p6,p7,p8}/n;
 SparseArray[{
{i,i}->pp1,  {i,j}->pp2,{i,k}->pp3, {i,l}->pp4,
{j,i}->-pp2,{j,j}->pp1,{j,k}->pp4,{j,l}->-pp3,
{k,i}->-pp3,{k,j}->-pp4,{k,k}->pp1,{k,l}->pp2,
{l,i}->-pp4,{l,j}->pp3,{l,k}->-pp2,{l,l}->pp1,
Band[{1,1}]->1
},
{m,m}
]
]

Tests

Here is a standard 2×2 cyclic sweep through the Jacobi using the simple zeroing strategy!  

We need to quantify the convergence!  We need to quantify how diagonal things need to be? We need to test other startegies for the sweep.

2x2 Jacobi (sweep)

m=4;
A= RandomReal[{-1,1},{m,m}];
A=A+Transpose[A];
AA=A;
Table[
Do[
{ii,jj}={i,j};
JJ=J[{ii,jj},AA[[ii,{ii,jj}]],m];
AA=JJ.AA.Transpose[JJ],
{i,1,m-1},
{j,i+1,m}
];
AA,
{5}]
{MatrixForm[A], MatrixForm[AA]}
Eigenvalues[A]

2x2 Jacobi (random)

Here is a random set of iterations.  Is it better or worse?

m=4;
A= RandomReal[{-1,1},{m,m}];
A=A+Transpose[A];
AA=A;
Table[
Do[
{ii,jj}=RandomSample[{1,2,3,4},2];
JJ=J[{ii,jj},AA[[ii,{ii,jj}]],m];
AA=JJ.AA.Transpose[JJ],
{i,1,m-1},
{j,i+1,m}
];
AA,
{5}]
{MatrixForm[A], MatrixForm[AA]}
Eigenvalues[A]

4x4 Jacobi

Here is an iteration using the 4x4 Jacobi in a purely sequential manner.  It does not seem to wokr very well at all!

m=4;
A= RandomReal[{-1,1},{m,m}];
A=A+Transpose[A];
AA=A;
Do[
Do[
{ii,jj,kk,ll}=RotateLeft[{1,2,3,4},i];
(* Print[{ii,jj,kk,ll}]; *)
JJ=J[{ii,jj,kk,ll},AA[[ii,{ii,jj,kk,ll}]],m];
AA=JJ.AA.Transpose[JJ],
{i,0,m-1}],
{10}]
Map[MatrixForm,{A,AA}]
Map[Eigenvalues,{A,AA}]
Map[MatrixPlot,{A,AA}]

4x4 Jacobi (again)

Here is an iteration using the 4x4 Jacobi by picking out the worst colum.  It does not seem to work very well at all.

FindLeastDiag[A_,m_]:=Module[
{RelativeBadness,i,Worst},
RelativeBadness=Table[ ((Norm[A[[All,i]]])/(Abs[A[[i,i]]])),{i,1,m}];
Worst=Max[RelativeBadness];
{Worst,Position[RelativeBadness,Worst,1][[1,1]]}
]
FindLeastDiag[A,m]
FindLeastDiag[AA,m]

m=4;
A= RandomReal[{-1,1},{m,m}];
A=A+Transpose[A];
AA=A;
Do[
Do[
{err,ii}=FindLeastDiag[AA,m];
{jj,kk,ll}=DeleteCases[{1,2,3,4},ii];
Print[{ii,jj,kk,ll}]; 
JJ=J[{ii,jj,kk,ll},AA[[ii,{ii,jj,kk,ll}]],m];
AA=JJ.AA.Transpose[JJ],
{i,0,1}],
{1}]
Map[MatrixForm,{A,AA}]
Map[Eigenvalues,{A,AA}]
Map[MatrixPlot,Abs[{A,AA}]]

MatrixForm[AA.JJ]
MatrixForm[Transpose[JJ].AA.JJ]
MatrixForm[AA]

Can I make a small angle approximation to the similarity zeroer?  This would be appropriate for a diagonally dominant matrix cause then it would only require small angles. 
First for a 2x2 I can actually make the off diagonals zero!

A=Array[a,{2,2}]/.a[i_,j_]:>If[i>j,a[j,i],a[i,j]];
JJ=J[{1,2},{c,s},2];
ANew=Simplify[(JJ.A.Transpose[JJ]),Assumptions->{(c^2)+(s^2)==1,Element[c,Reals],Element[s,Reals]}];
MatrixForm[ANew]

ANew[[1,1]]

ANew[[2,2]]+ANew[[1,1]]//Simplify
ANew[[2,2]]-ANew[[1,1]]//Simplify

s/.Simplify[
Solve[
Simplify[ANew[[2,1]]/. a[1,2]->β a[1,1]]==0,s]
]
s/.Solve[
Normal[Series[Simplify[ANew[[2,1]]/. a[1,2]->β a[1,1]]==0,{β,0,4}]],
s]
s/.Solve[
Normal[Series[Simplify[ANew[[2,1]]/. a[1,2]->β a[1,1]]==0,{s,0,1}]],
s]

ϵ/.Solve[
Simplify[Normal[Series[Simplify[ANew[[2,1]]/. a[1,2]->β a[1,1]]==0,{ϵ,0, 2}]]],
ϵ]

Simplify[ϵ/.Solve[a[1,2]-2 (ϵ^2) a[1,2]+ϵ (Sqrt[1-(ϵ^2)]) (-a[1,1]+a[2,2])==0,ϵ]]

Series[
Simplify[ϵ/.Solve[a[1,2]-2 (ϵ^2) a[1,2]+ϵ (Sqrt[1-(ϵ^2)]) (-a[1,1]+a[2,2])==0,ϵ]/.a[1,2]->β a[1,1]],
{β,0,2}]

?Simplify

Sweep Patterns for Similarity Transformations