' ******************************************************************
' This program computes the eigenvalues and eigenvectors of a
' symmetric matrix by singular value decomposition or by the
' iterative method of Jacobi. The methods are demonstrated with
' Hilbert matrices. These matrices are ill-conditioned (i.e. the
' ratio of the lowest to the highest eigenvalue is very low).
'
' The Jacobi method applies a series of rotations to the original
' matrix in order to transform it into a diagonal matrix. The
' diagonal terms of this matrix are the eigenvalues. The product
' of the rotation matrices gives the eigenvectors. The original
' matrix is destroyed during the process.
'
' The constant Tol defines the tolerance with which an off-diagonal
' element of the transformed matrix is considered zero (expressed as
' a fraction of the sum of squared diagonal terms). The constant
' MaxIter defines the maximal number of iterations allowed. These
' two constants are linked, i.e. decreasing Tol may need increasing
' MaxIter to avoid non-convergence of the Jacobi procedure.
'
' The SVD method uses the decomposition A = U*S*V' where S is the
' diagonal matrix of eigenvalues and V is the matrix of eigenvectors
'
' These two methods destroy the original matrix A.
' ******************************************************************

#INCLUDE "fbmath.bi"

CONST MaxIter = 1000  ' Maximum number of iterations (Jacobi)
CONST Tol     = 1D-8  ' Required precision (Jacobi)

SUB Hilbert(N AS INTEGER, A() AS DOUBLE)
' ------------------------------------------------------------------
' Generates an Hilbert matrix of order N
' ------------------------------------------------------------------
'       (                                             )
'       ( 1      1/2     1/3     1/4     ... 1/N      )
'       ( 1/2    1/3     1/4     1/5     ... 1/(N+1)  )
'   A = ( 1/3    1/4     1/5     1/6     ... 1/(N+2)  )
'       ( ........................................... )
'       ( 1/N    1/(N+1) 1/(N+2) 1/(N+3) ... 1/(2N-1) )
'       (                                             )
' ------------------------------------------------------------------

  DIM AS INTEGER I, J

  A(1, 1) = 1
  FOR J = 2 TO N                 ' First row of matrix
    A(1, J) = 1 / J
  NEXT J

  FOR I = 2 TO N
    A(I, N) = 1 / (N + I - 1)    ' Last column of matrix
    FOR J = 1 TO N - 1
      A(I, J) = A(I - 1, J + 1)  ' Fill matrix
    NEXT J
  NEXT I
END SUB

SUB WriteResults(Lambda() AS DOUBLE, V() AS DOUBLE)
' ------------------------------------------------------------------
' Outputs results to screen
' ------------------------------------------------------------------

  DIM AS INTEGER I, J, N

  N = UBOUND(Lambda)

  PRINT : PRINT "Eigenvalues :" : PRINT

  FOR I = 1 TO N
    PRINT Lambda(I)
  NEXT I

  IF N < 8 THEN
    PRINT : PRINT "Eigenvectors (columns) :" : PRINT
    FOR I = 1 TO N
      FOR J = 1 TO N
        PRINT USING "#####.######"; V(I, J);
      NEXT J
      PRINT
    NEXT I
  END IF
END SUB

DIM AS INTEGER N  ' Size of matrix

DO
  PRINT
  INPUT "Order of Hilbert matrix (1 to end) : ", N

  IF N < 2 THEN STOP

  ' Allocate vectors and matrices
  REDIM AS DOUBLE A(1 TO N, 1 TO N)
  REDIM AS DOUBLE V(1 TO N, 1 TO N)
  REDIM AS DOUBLE Lambda(1 TO N)

  ' Generate Hilbert matrix
  Hilbert N, A()

  ' Compute eigenvalues and eigenvectors by Jacobi method
  ' Jacobi A(), MaxIter, Tol, Lambda(), V()

  ' Compute eigenvalues and eigenvectors by SVD method
  EigenSym A(), Lambda(), V()

  IF MathErr = MatOk THEN
    WriteResults Lambda(), V()
  ELSE
    PRINT "Too many iterations!"
  END IF
LOOP
