Jacobi Eigenvalue Algorithm - Algorithm

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(SRn×n; out eRn; out ERn×n) var i, k, l, m, stateN s, c, t, p, yR indNn changedLn function maxind(kN) ∈ N ! index of largest off-diagonal element in row k m := k+1 for i := k+2 to n do ifSki│ > │Skmthen m := i endif endfor return m endfunc procedure update(kN; tR) ! 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 (yek) then changedk := true; state := state+1 endif endproc procedure rotate(k,l,i,jN) ! perform rotation of Sij, Skl ┐ ┌ ┐┌ ┐ │Skl│ │cs││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 ifSk indk│ > │Sm indmthen m := k endif endfor k := m; l := indm; p := Skl ! calculate c = cos φ, s = sin φ y := (elek)/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│ │cs││Eki│ │ │ := │ ││ │ │Eli│ │s c││Eli│ └ ┘ └ ┘└ endfor ! rows k, l have changed, update rows indk, indl indk := maxind(k); indl := maxind(l) loop endproc

Read more about this topic:  Jacobi Eigenvalue Algorithm