# Jacobi eigenvalue and eigenvectors algorithm

Suppose we have a symmetric matrix with dimensions n x n. I need to find the eigenvalues and eigenvectors of this matrix using the Jacobi method. I wrote this code:

j[A0_List] :=
Module[{A = A0, c, s, tau, t, p, q, n = Length@A0, m, R},
V = IdentityMatrix[n];
m = Abs[A - DiagonalMatrix@Table[A[[i, i]], {i, n}]]; (*got rid of the diagonal*)
{p, q} = First@Position[m, Max@m]; (*position of the absolute greater element off the diagonal*)
tau = (A[[q, q]] - A[[p, p]])/(2*A[[p, q]]);
Which[tau < 10^19 && tau >= 0, t = 1/(tau + Sqrt[1 + tau^2]),
tau >= 10^19, t = 1/Abs[2 tau], tau < 0,
t = -1/(-tau + Sqrt[1 + tau^2])];
c = 1/Sqrt[1 + t^2]; s = t c; (*c is cosine and s is sine*)
R = {{c, s}, {-s, c}}; (*rotation matrix*)
A[[All, {p, q}]] = A[[All, {p, q}]].R;  (*implementing the transformation*)
A[[{p, q}, All]] = Transpose[R].A[[{p, q}, All]];
V[[All, {p, q}]] = V[[All, {p, q}]].R;  (*here the eigenvectors are supposed to compute*)
A];
jacobi[A0_List?SymmetricMatrix] :=
Module[{A = N[A0], list, n = Length@A0, zero},
zero = ConstantArray[0, {n, n}];
list = FixedPointList[Chop[j[#], 10^(-5)] &, A,
SameTest -> (Chop[j[#], 10^(-5)] -
DiagonalMatrix@Table[Chop[j[#], 10^(-5)][[i, i]], {i, n}] ==
zero &)]; (*iteration unitl the last matix has the off diagonal elements zero*)
Last@list (*print the matrix that has the eigenvalues on the diagonal*)]


Can someone please help me understand how can I now find the eigenvectors (something else than (A- λI)v = 0), because mine don't come out right. I would also appreciate any comment on my code. Thank you.

-
FixedPointList[] seems to be a poor fit for this. Have you tried starting from a procedural implementation? –  Ｊ. Ｍ. Apr 4 at 18:11
If you have not already done so, might have a look at Demmel's "Applied Numerical Linear Algebra" section 5.4.3. Also it looks like you might be rotating the 2x2 overlap, that is, A[[{p,q},{p,q}]], twice. Do you really want the second rotation of A? The one-sided version in the above reference does not do that. –  Daniel Lichtblau Apr 4 at 18:28
I try to avoid procedural programming because of the longer computation time. –  Gregor Apr 4 at 19:16
Maybe, but you want something correct to compare against, no? And it is not true that procedural methods give longer computational times; it's entirely dependent on the cleverness of the programmer. –  Ｊ. Ｍ. Apr 6 at 13:26