Algorithm
The following algorithm is a description of the Jacobi method in math-like notation. It calculates a vector e which contains the eigenvalues and a matrix E which contains the corresponding eigenvectors, i.e. is an eigenvalue and the column an orthonormal eigenvector for, i = 1, …, n.
procedure jacobi(S ∈ Rn×n; out e ∈ Rn; out E ∈ Rn×n) var i, k, l, m, state ∈ N s, c, t, p, y ∈ R ind ∈ Nn changed ∈ Ln function maxind(k ∈ N) ∈ N ! index of largest off-diagonal element in row k m := k+1 for i := k+2 to n do if │Ski│ > │Skm│ then m := i endif endfor return m endfunc procedure update(k ∈ N; t ∈ R) ! update ek and its status y := ek; ek := y+t if changedk and (y=ek) then changedk := false; state := state−1 elsif (not changedk) and (y≠ek) then changedk := true; state := state+1 endif endproc procedure rotate(k,l,i,j ∈ N) ! perform rotation of Sij, Skl ┌ ┐ ┌ ┐┌ ┐ │Skl│ │c −s││Skl│ │ │ := │ ││ │ │Sij│ │s c││Sij│ └ ┘ └ ┘└ ┘ endproc ! init e, E, and arrays ind, changed E := I; state := n for k := 1 to n do indk := maxind(k); ek := Skk; changedk := true endfor while state≠0 do ! next rotation m := 1 ! find index (k,l) of pivot p for k := 2 to n−1 do if │Sk indk│ > │Sm indm│ then m := k endif endfor k := m; l := indm; p := Skl ! calculate c = cos φ, s = sin φ y := (el−ek)/2; t := │y│+√(p2+y2) s := √(p2+t2); c := t/s; s := p/s; t := p2/t if y<0 then s := −s; t := −t endif Skl := 0.0; update(k,−t); update(l,t) ! rotate rows and columns k and l for i := 1 to k−1 do rotate(i,k,i,l) endfor for i := k+1 to l−1 do rotate(k,i,i,l) endfor for i := l+1 to n do rotate(k,i,l,i) endfor ! rotate eigenvectors for i := 1 to n do ┌ ┐ ┌ ┐┌ ┐ │Eki│ │c −s││Eki│ │ │ := │ ││ │ │Eli│ │s c││Eli│ └ ┘ └ ┘└ ┘ endfor ! rows k, l have changed, update rows indk, indl indk := maxind(k); indl := maxind(l) loop endprocRead more about this topic: Jacobi Eigenvalue Algorithm