of the pairs . Also the updating of the index array can be done in O(n) average-case complexity: First, the maximum entry in the updated rows k and l can be found in O(n) steps. In the other rows i, only the entries in columns k and l change. Looping over these rows, if is neither k nor l, it suffices to compare the old maximum at to the new entries and update if necessary. If should be equal to k or l and the corresponding entry decreased during the update, the maximum over row i has to be found from scratch in O(n) complexity. However, this will happen on average only once per rotation. Thus, each rotation has O(n) and one sweep O(n3) average-case complexity, which is equivalent to one matrix multiplication. Additionally the must be initialized before the process starts, which can be done in n2 steps. Typically the Jacobi method converges within numerical precision after a small number of sweeps. Note that multiple eigenvalues reduce the number of iterations since . 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, d, r ∈ 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;