' ******************************************************************
' Eigenvalues and eigenvectors of a symmetric matrix
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE MatOk 0  
' No error

#DEFINE MatNonConv -1
' Non-convergence

' ------------------------------------------------------------------
' Mathematical constant
' ------------------------------------------------------------------

#DEFINE Sqrt2Div2 0.7071067811865476#  
' Sqrt(2)/2

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' ******************************************************************

SUB Jacobi(A()      AS DOUBLE,_
           MaxIter  AS INTEGER,_ 
           Tol      AS DOUBLE,_ 
           Lambda() AS DOUBLE,_ 
           V()      AS DOUBLE)
' ------------------------------------------------------------------
' Eigenvalues and eigenvectors of a symmetric matrix by the
' iterative method of Jacobi
' ------------------------------------------------------------------
' Input parameters  : A       = matrix
'                     MaxIter = maximum number of iterations
'                     Tol     = required precision
' ------------------------------------------------------------------
' Output parameters : Lambda = eigenvalues in decreasing order
'                     V      = matrix of eigenvectors (stored by
'                              columns)
' ------------------------------------------------------------------
' Possible results  : MatOK
'                     MatNonConv
' ------------------------------------------------------------------
' NB : 1. The eigenvectors are normalized, with their first
'         component > 0
'      2. This procedure destroys the original matrix A
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, J, K, Im, Jm, Iter
  DIM AS DOUBLE  B, C, C2, Na, Nd, P, Q, R, S, S2, T
  
  Lb = LBOUND(A, 1)
  Ub = UBOUND(A, 1)

  Iter = 0
  Na = 0
  Nd = 0
  R = 0

  FOR I = Lb TO Ub
    V(I, I) = 1
    Nd = Nd + A(I, I) * A(I, I)
    IF I <> Ub THEN
      FOR J = I + 1 TO Ub
        R = R + A(I, J) * A(I, J)
        V(I, J) = 0
        V(J, I) = 0
      NEXT J
    END IF
  NEXT I

  Na = Nd + 2 * R

  DO
    R = 0
    FOR I = Lb TO Ub - 1
      FOR J = I + 1 TO Ub
        T = ABS(A(I, J))
        IF T > R THEN
          R = T
          Im = I
          Jm = J
        END IF
      NEXT J
    NEXT I

    B = A(Im, Im) - A(Jm, Jm)

    IF B = 0 THEN
      C = Sqrt2div2
      S = C * SGN(A(Im, Jm))
    ELSE
      P = 2 * A(Im, Jm) * SGN(B)
      Q = ABS(B)
      R = SQR(P * P + Q * Q)
      C = SQR(.5 * (1 + Q / R))
      S = .5 * P / (R * C)
    END IF

    FOR K = Lb TO Ub
      R = V(K, Im)
      V(K, Im) = C * R + S * V(K, Jm)
      V(K, Jm) = C * V(K, Jm) - S * R
    NEXT K

    IF Im <> Lb THEN
      FOR K = Lb TO Im - 1
        R = A(K, Im)
        A(K, Im) = C * R + S * A(K, Jm)
        A(K, Jm) = C * A(K, Jm) - S * R
      NEXT K
    END IF

    IF Jm <> Im + 1 THEN
      FOR K = Im + 1 TO Jm - 1
        R = A(Im, K)
        A(Im, K) = C * R + S * A(K, Jm)
        A(K, Jm) = C * A(K, Jm) - S * R
      NEXT K
    END IF

    IF Jm <> Ub THEN
      FOR K = Jm + 1 TO Ub
        R = A(Im, K)
        A(Im, K) = C * R + S * A(Jm, K)
        A(Jm, K) = C * A(Jm, K) - S * R
      NEXT K
    END IF

    Nd = Nd + 2 * A(Im, Jm) * A(Im, Jm)

    C2 = C * C
    S2 = S * S
    P = 2 * S * C * A(Im, Jm)
    R = A(Im, Im)
    A(Im, Im) = C2 * R + S2 * A(Jm, Jm) + P
    A(Jm, Jm) = S2 * R + C2 * A(Jm, Jm) - P
    A(Im, Jm) = 0

    Iter = Iter + 1
    IF Iter > MaxIter THEN ErrCode = MatNonConv
  LOOP UNTIL ABS(1 - Na / Nd) < Tol

  ' The diagonal terms of the transformed matrix are the eigenvalues
  FOR I = Lb TO Ub
    Lambda(I) = A(I, I)
  NEXT I

  ' Sort eigenvalues and eigenvectors
  FOR I = Lb TO Ub - 1
    K = I
    R = Lambda(I)
    FOR J = I + 1 TO Ub
      IF Lambda(J) > R THEN
        K = J
        R = Lambda(J)
      END IF
    NEXT J

    ' Exchange eigenvalues
    SWAP Lambda(I), Lambda(K)

    ' Exchange eigenvectors
    FOR J = Lb TO Ub
      SWAP V(J, I), V(J, K)
    NEXT J
  NEXT I

  ' Make sure that the first component of each eigenvector is > 0
  FOR J = Lb TO Ub
    IF V(Lb, J) < 0 THEN
      FOR I = Lb TO Ub
        V(I, J) = -V(I, J)
      NEXT I
    END IF
  NEXT J

  ErrCode = MatOk
END SUB

