' ******************************************************************
' Eigenvalues and eigenvectors of a positive semi-definite
' symmetric matrix by singular value decomposition
' ******************************************************************

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

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

' ------------------------------------------------------------------
' External subroutine
' ------------------------------------------------------------------

DECLARE SUB SV_Decomp(A() AS DOUBLE, S() AS DOUBLE, V() AS DOUBLE)
' Singular value decomposition

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

SUB EigenSym(A()      AS DOUBLE,_
             Lambda() AS DOUBLE,_
             V()      AS DOUBLE)
' ------------------------------------------------------------------
' Eigenvalues and eigenvectors of a positive semi-definite
' symmetric matrix by singular value decomposition
' ------------------------------------------------------------------
' Input parameters  : A       = matrix
' ------------------------------------------------------------------
' Output parameters : Lambda = eigenvalues in decreasing order
'                     V      = matrix of eigenvectors (stored by
'                              columns)
' ------------------------------------------------------------------
' Possible results  : MatOK
'                     MatErrDim
'                     MatNonConv
' ------------------------------------------------------------------
' NB : 1. The eigenvectors are normalized, with their first
'         component > 0
'      2. This procedure destroys the original matrix A
' ------------------------------------------------------------------

  DIM AS INTEGER I, J, K, Lb, Ub
  DIM AS DOUBLE  R

  Lb = LBOUND(A, 1)
  Ub = UBOUND(A, 1)

  SV_Decomp A(), Lambda(), V()

  IF ErrCode <> 0 THEN EXIT SUB

  ' 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
END SUB

