DECLARE SUB CDIV (AR#, AI#, BR#, BI#, ZR#, ZI#)
DECLARE SUB CPRO (alpha#, c#(), c1#())
DECLARE SUB CONJ (c#(), c1#())
DECLARE SUB CMUL (c1#(), c2#(), c3#())
DECLARE SUB CDIF (c1#(), c2#(), c3#())
DECLARE SUB CADD (c1#(), c2#(), c3#())
DECLARE FUNCTION Squ# (A#)
DECLARE FUNCTION CABS# (c#())
'***********************************************************************
'*     Calculate eigenvalues and eigenvectors of a Square Hermitian    *
'*     Matrix By Jacobi's Method                                       *
'* ------------------------------------------------------------------- *
'* SAMPLE RUN:                                                         *
'* File Matc4.dat contains:                                            *
'* 4                                                                   *
'* 12   0   -1   -1   2   0   3   3                                    *
'* -1   1   12    0   1  -1  -2   0                                    *
'*  2   0    1    1   8   0  -1  -1                                    *
'*  3  -3   -2    0  -1   1   8   0                                    *
'*                                                                     *
'*  Precision: 1D-10                                                   *
'*  Max. number of iterations: 25                                      *
'*                                                                     *
'*  Eigenvalue  1:  4.000000                                           *
'*  Eigenvector:                                                       *
'*  -0.500000 -  0.500000 I                                            *
'*   0.000000 +  0.000000 I                                            *
'*   0.500000 +  0.500000 I                                            *
'*   1.000000 +  0.000000 I                                            *
'*                                                                     *
'*  Eigenvalue  2:  8.000000                                           *
'*  Eigenvector:                                                       *
'*  -0.000000 +  0.000000 I                                            *
'*  -0.500000 +  0.500000 I                                            *
'*   1.000000 +  0.000000 I                                            *
'*  -0.500000 +  0.500000 I                                            *
'*                                                                     *
'*  Eigenvalue  3: 12.000000                                           *
'*  Eigenvector:                                                       *
'*   0.500000 +  0.500000 I                                            *
'*   1.000000 +  0.000000 I                                            *
'*   0.500000 +  0.500000 I                                            *
'*   0.000000 +  0.000000 I                                            *
'*                                                                     *
'*  Eigenvalue  4: 16.000000                                           *
'*  Eigenvector:                                                       *
'*   1.000000 +  0.000000 I                                            *
'*  -0.500000 +  0.500000 I                                            *
'*  -0.000000 +  0.000000 I                                            *
'*   0.500000 -  0.500000 I                                            *
'*                                                                     *
'* ------------------------------------------------------------------- *
'* Reference: "ALGEBRE Algorithmes et programmes en Pascal             *
'*             By Jean-Louis Jardrin - Dunod BO-PRE 1988" [BIBLI 10]   *
'*                                                                     *
'*                            Basic Release 1.1 By J-P Moreau, Paris.  *
'*                                      (www.jpmoreau.fr)              *
'*                                                                     *
'* Release 1.1 (Nov. 30th, 2007): Added function NORMAL.               *
'***********************************************************************
'Program TEPHJ

DEFDBL A-H, O-Z
DEFINT I-N

ISIZE = 25

'integer I,it,J,M,N
'double dta

DIM A(ISIZE, ISIZE, 2), VX(ISIZE, ISIZE, 2) 'Complex matrices
DIM R(ISIZE)

GOSUB 500  'call LME(N,A)    (read complex matrix from text file)

CLS
PRINT
PRINT "  ** Compute Eigenvalues and Eigenvectors **"
PRINT "         of a Square Hermitian Matrix"
PRINT "             By Jacobi's Method"
PRINT

INPUT "  Precision: ", dta
INPUT "  Max. number of iterations: ", M
PRINT

GOSUB 1000  'call EPHJ(dta, M, N, A, it, R, VX)

GOSUB 2000  'call NORMAL(N, R, VX)

F$ = "##.######"

IF it = -1 THEN
  PRINT "  No convergence."
ELSE                             'display results
  FOR J = 1 TO N
    PRINT "  Eigenvalue "; J; ": ";
    PRINT USING F$; R(J)
    PRINT "  Eigenvector:"
    FOR I = 1 TO N
      PRINT "  ";
      PRINT USING F$; VX(I, J, 1);
      IF VX(I, J, 2) >= 0# THEN
	PRINT " + ";
      ELSE
	PRINT " - ";
      END IF
      PRINT USING F$; ABS(VX(I, J, 2));
      PRINT " I"
    NEXT I
    PRINT
  NEXT J
  PRINT
END IF

END  'of main program


500 'Sub LME(N,A)
'Read from input file complex hermitian matrix
  OPEN "matc4.dat" FOR INPUT AS #1
  INPUT #1, N
  FOR I = 1 TO N
    FOR J = 1 TO N
      INPUT #1, A(I, J, 1), A(I, J, 2)
    NEXT J
  NEXT I
  CLOSE #1
RETURN


1000 'Sub EPHJ(dta, M, N, A, it, R, VX)
'************************************************************************
'*   Compute all eigenvalues/eigenvectors of a square hermitian matrix  *
'*   using Jacobi's method.                                             *
'* -------------------------------------------------------------------- *
'* Inputs:                                                              *
'*         dta: required precision (double)                             *
'*         M  : maximum number of iterations (integer)                  *
'*         N  : size of matrix                                          *
'*         A  : pointer to hermitian (complex) matrix                   *
'* Outputs:                                                             *
'*         it : flag for convergence (integer) =-1, no convergence      *
'*         R  : vector(1%%N) of eigenvalues (double)                    *
'*         VX : matrix storing complex eigenvectors in columns          *
'************************************************************************

' complex numbers
DIM c0(2), c1(2), c2(2), c3(2), u0(2), tmp(2), u1(2), z0(2), Z1(2)

  z0(1) = 0#: z0(2) = 0#: Z1(1) = 1#: Z1(2) = 0#
  FOR I = 1 TO N
    FOR J = 1 TO N
      IF I = J THEN
	VX(I, J, 1) = Z1(1)
	VX(I, J, 2) = Z1(2)
      ELSE
	VX(I, J, 1) = z0(1)
	VX(I, J, 2) = z0(2)
      END IF
    NEXT J
  NEXT I
  it = -1: L = 1
  WHILE L <= M AND it <> 1
    s = 0#
    FOR I = 1 TO N
      FOR J = I + 1 TO N
	tmp(1) = A(I, J, 1): tmp(2) = A(I, J, 2)
	t0 = CABS(tmp())
	IF t0 > s THEN
	  s = t0: K0 = I: L0 = J
	END IF
      NEXT J
    NEXT I
    IF s = 0# THEN
      it = 1
    ELSE
      tmp(1) = A(K0, L0, 1): tmp(2) = A(K0, L0, 2)
      delta = Squ(A(L0, L0, 1) - A(K0, K0, 1)) + 4# * Squ(CABS(tmp()))
      t0 = A(L0, L0, 1) - A(K0, K0, 1) + SQR(delta)
      t1 = A(L0, L0, 1) - A(K0, K0, 1) - SQR(delta)
      IF ABS(t0) >= ABS(t1) THEN
	w0 = t0
      ELSE
	w0 = t1
      END IF
      s0 = ABS(w0) / SQR(Squ(w0) + 4# * Squ(CABS(tmp())))
      t0 = 2# * s0 / w0:
      CPRO t0, tmp(), c0()
      CONJ c0(), c1()
      FOR I = 1 TO K0 - 1
	u0(1) = A(I, K0, 1): u0(2) = A(I, K0, 2)
	CMUL c0(), u0(), c2(): tmp(1) = A(I, L0, 1): tmp(2) = A(I, L0, 2)
	CPRO s0, tmp(), c3(): CADD c2(), c3(), tmp()
	A(I, K0, 1) = tmp(1): A(I, K0, 2) = tmp(2)
	tmp(1) = A(I, L0, 1): tmp(2) = A(I, L0, 2)
	CMUL c1(), tmp(), c2(): CPRO s0, u0(), c3(): CDIF c2(), c3(), tmp()
	A(I, L0, 1) = tmp(1): A(I, L0, 2) = tmp(2)
      NEXT I
      FOR K = K0 + 1 TO L0 - 1
	u0(1) = A(K0, K, 1): u0(2) = A(K0, K, 2)
	CMUL c1(), u0(), c2(): tmp(1) = A(K, L0, 1): tmp(2) = A(K, L0, 2)
	CONJ tmp(), u1(): CPRO s0, u1(), c3()
	CADD c2(), c3(), tmp(): A(K0, K, 1) = tmp(1): A(K0, K, 2) = tmp(2)
	tmp(1) = A(K, L0, 1): tmp(2) = A(K, L0, 2)
	CMUL c1(), tmp(), c2(): CONJ u0(), u1(): CPRO s0, u1(), c3()
	CDIF c2(), c3(), tmp()
	A(K, L0, 1) = tmp(1): A(K, L0, 2) = tmp(2)
      NEXT K
      FOR J = L0 + 1 TO N
	u0(1) = A(K0, J, 1): u0(2) = A(K0, J, 2)
	CMUL c1(), u0(), c2(): tmp(1) = A(L0, J, 1): tmp(2) = A(L0, J, 2)
	CPRO s0, tmp(), c3(): CADD c2(), c3(), tmp()
	A(K0, J, 1) = tmp(1): A(K0, J, 2) = tmp(2)
	tmp(1) = A(L0, J, 1): tmp(2) = A(L0, J, 2)
	CMUL c0(), tmp(), c2(): CPRO s0, u0(), c3(): CDIF c2(), c3(), tmp()
	A(L0, J, 1) = tmp(1): A(L0, J, 2) = tmp(2)
      NEXT J
      t0 = A(K0, K0, 1)
      tmp(1) = A(K0, L0, 1): tmp(2) = A(K0, L0, 2)
      t1 = 4# * Squ(s0 * CABS(tmp())) / w0
      A(K0, K0, 1) = Squ(CABS(c0())) * t0 + t1 + Squ(s0) * A(L0, L0, 1)
      A(L0, L0, 1) = Squ(s0) * t0 - t1 + Squ(CABS(c0())) * A(L0, L0, 1)
      A(K0, L0, 1) = z0(1): A(K0, L0, 2) = z0(2)
      FOR I = 1 TO N
	u0(1) = VX(I, K0, 1): u0(2) = VX(I, K0, 2)
	CMUL c0(), u0(), c2(): tmp(1) = VX(I, L0, 1): tmp(2) = VX(I, L0, 2)
	CPRO s0, tmp(), c3(): CADD c2(), c3(), tmp()
	VX(I, K0, 1) = tmp(1): VX(I, K0, 2) = tmp(2)
	tmp(1) = VX(I, L0, 1): tmp(2) = VX(I, L0, 2)
	CMUL c1(), tmp(), c2(): CPRO s0, u0(), c3(): CDIF c2(), c3(), tmp()
	VX(I, L0, 1) = tmp(1): VX(I, L0, 2) = tmp(2)
      NEXT I
      t0 = 0#
      FOR I = 1 TO N - 1
	FOR J = I + 1 TO N
	  tmp(1) = A(I, J, 1): tmp(2) = A(I, J, 2)
	  t0 = t0 + Squ(CABS(tmp()))
	NEXT J
      NEXT I
      s = 2# * t0
      IF (s < dta) THEN
	it = 1
      ELSE
	L = L + 1
      END IF
    END IF
  WEND
  IF it = 1 THEN
    FOR I = 1 TO N
      R(I) = A(I, I, 1)
    NEXT I
  END IF
RETURN

2000 'Sub NORMAL(N, R, VX)
'--------------------------------------------------------------
'   Sort eigenvalues by fabsolute values in ascending order and
'   normalize.
'--------------------------------------------------------------
DIM Tr(ISIZE), Ti(ISIZE)
'double ZM,Z1,Z2,VR
DIM V(2)                'Complex number

'   SORT SOLUTIONS IN ASCENDING ORDER
    FOR J = 2 TO N
      VR = R(J)
      FOR K = 1 TO N
	Tr(K) = VX(K, J, 1)
	Ti(K) = VX(K, J, 2)
      NEXT K

      FOR I = J - 1 TO 1 STEP -1
	IF ABS(R(I)) <= ABS(VR) THEN GOTO 5
	R(I + 1) = R(I)
	FOR K = 1 TO N
	  VX(K, I + 1, 1) = VX(K, I, 1)
	  VX(K, I + 1, 2) = VX(K, I, 2)
	NEXT K
      NEXT I
      I = 0
5     R(I + 1) = VR
      FOR K = 1 TO N
	VX(K, I + 1, 1) = Tr(K)
	VX(K, I + 1, 2) = Ti(K)
      NEXT K
    NEXT J

'   NORMALIVXE WITH RESPECT TO BIGGEST ELEMENT
    FOR J = N TO 1 STEP -1
      ZM = 0#
      FOR I = 1 TO N
	V(1) = VX(I, J, 1)
	V(2) = VX(I, J, 2)
	Z1 = CABS(V())
	IF Z1 >= ZM THEN
	  IM = I
	  ZM = Z1
	END IF
      NEXT I
      Z1 = VX(IM, J, 1)
      Z2 = VX(IM, J, 2)
      FOR I = 1 TO N
	CDIV VX(I, J, 1), VX(I, J, 2), Z1, Z2, V(1), V(2)
	VX(I, J, 1) = V(1)
	VX(I, J, 2) = V(2)
      NEXT I
    NEXT J
RETURN

FUNCTION CABS (c())
'Absolute value of a complex number
  CABS = SQR(Squ(c(1)) + Squ(c(2)))
END FUNCTION

SUB CADD (c1(), c2(), c3())
'Add two complex numbers
  c3(1) = c1(1) + c2(1): c3(2) = c1(2) + c2(2)
END SUB

SUB CDIF (c1(), c2(), c3())
'Substract two complex numbers
  c3(1) = c1(1) - c2(1): c3(2) = c1(2) - c2(2)
END SUB

SUB CDIV (AR, AI, BR, BI, ZR, ZI)
'      EFFECT DIVISION Z=(ZR+I*ZI)=(AR+I*AI)/(BR+I*BI)
' -----DO NOT USE IF BR=BI=0
' double YR,YI,W
  YR = BR
  YI = BI
  IF ABS(YR) > ABS(YI) THEN
    W = YI / YR
    YR = W * YI + YR
    ZR = (AR + W * AI) / YR
    ZI = (AI - W * AR) / YR
  ELSE
    W = YR / YI
    YI = W * YR + YI
    ZR = (W * AR + AI) / YI
    ZI = (W * AI - AR) / YI
  END IF
END SUB

SUB CMUL (c1(), c2(), c3())
'Multiply two complex numbers
  c3(1) = c1(1) * c2(1) - c1(2) * c2(2)
  c3(2) = c1(1) * c2(2) + c1(2) * c2(1)
END SUB

SUB CONJ (c(), c1())
'Return conjugate of a complex number
  c1(1) = c(1): c1(2) = -c(2)
END SUB

SUB CPRO (alpha, c(), c1())
'Multiply a complex number by a real number
  c1(1) = alpha * c(1): c1(2) = alpha * c(2)
END SUB

FUNCTION Squ (A)
  Squ = A * A
END FUNCTION