DECLARE SUB Data1 (n%, eps#, d1#, d2#, m%)
DECLARE SUB EPMRI (eps#, d1#, d2#, m%, n%, it%, R#())
DECLARE SUB DECCRM (eps#, n%, it%)
DECLARE SUB VAMR (eps#, dta#, m%, n%, it%, R#())
DECLARE SUB IIM (eps#, dta#, m%, n%, it%, gamma#, X1#())
'********************************************************
'* Eigenvalues and eigenvectors of a real square matrix *
'* by Rutishauser's method and inverse iteration method *
'* ---------------------------------------------------- *
'* Reference:                                           *
'*                                                      *
'*   "ALGEBRE Algorithmes et programmes en Pascal       *
'*    de Jean-Louis Jardrin - Dunod BO-PRE 1988."       *
'*    [BIBLI 10].                                       *
'*                                                      *
'*                  Basic Release By J-P Moreau, Paris. *
'*                           (www.jpmoreau.fr)          *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'* Input file: elpro.dat                                *
'*                                                      *
'* 5                                                    *
'*   1   2   3  -7    12                                *
'*   2   4   7   3    -1                                *
'*   3   7  10   8     4                                *
'*  -7   3   8  -0.75 -9                                *
'*  12  -1   4  -9    10                                *
'*                                                      *
'* Output to screen:                                    *
'*                                                      *
'*    *** EIGENVALUES AND EIGENVECTORS ***              *
'*          OF A REAL SQUARE MATRIX                     *
'*          BY RUTISHAUSER''S METHOD                    *
'*        AND INVERSE ITERATION METHOD                  *
'*                                                      *
'* Matrix A:                                            *
'*   1.0000   2.0000   3.0000  -7.0000  12.0000         *
'*   2.0000   4.0000   7.0000   3.0000  -1.0000         *
'*   3.0000   7.0000  10.0000   8.0000   4.0000         *
'*  -7.0000   3.0000   8.0000  -0.7500  -9.0000         *
'*  12.0000  -1.0000   4.0000  -9.0000  10.0000         *
'*                                                      *
'* Eigenvalue 1 :  23.75595                             *
'* Eigenvector:                                         *
'*   0.70582  -0.01268   0.13149  -0.52750   1.00000    *
'*                                                      *
'* Eigenvalue 2 : -10.48655                             *
'* Eigenvector:                                         *
'*   0.46717  -0.00795  -0.50783   1.00000   0.26443    *
'*                                                      *
'* Eigenvalue 3 :   0.46335                             *
'* Eigenvector:                                         *
'*   0.20881   1.00000  -0.47803  -0.27500  -0.21691    *
'*                                                      *
'* Eigenvalue 4 :  -7.77458                             *
'* Eigenvector:                                         *
'*   1.00000  -0.32610   0.20962  -0.14769  -0.81542    *
'*                                                      *
'* Eigenvalue 5 :   0.46335                             *
'* Eigenvector:                                         *
'*   0.20881   1.00000  -0.47803  -0.27500  -0.21691    *
'*                                                      *
'********************************************************
'PROGRAM Elpro
DEFINT I-N
DEFDBL A-H, O-Z

OPTION BASE 1

    DIM SHARED NMAX
    NMAX = 20

'   double d1,d2,eps
    DIM SHARED A(NMAX, NMAX), VX(NMAX, NMAX), U(NMAX, NMAX), V(NMAX, NMAX)
    DIM R(NMAX)

    CLS
    PRINT
    PRINT "    *** EIGENVALUES AND EIGENVECTORS ***"
    PRINT "          OF A REAL SQUARE MATRIX"
    PRINT "          BY RUTISHAUSER''S METHOD"
    PRINT "        AND INVERSE ITERATION METHOD"
    PRINT

'   read data from input text file
    Data1 n, eps, d1, d2, m

    EPMRI eps, d1, d2, m, n, it, R()

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

    IF it = -1 THEN
      PRINT "  No convergence!"
      PRINT
    ELSE
      FOR j = 1 TO n
	PRINT " Eigenvalue"; j; ":";
	PRINT USING F$; R(j)
	PRINT " Eigenvector:"
	FOR i = 1 TO n
	  PRINT USING F$; VX(i, j);
	NEXT i
	PRINT
	PRINT
      NEXT j
    END IF
    PRINT

END 'of main program

' Read data from input text file
  SUB Data1 (n, eps, d1, d2, m)
    
    OPEN "elpro.dat" FOR INPUT AS #1

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

    INPUT #1, n
    PRINT
    PRINT " Matrix A:"
    PRINT
    FOR i = 1 TO n
      FOR j = 1 TO n
	INPUT #1, A(i, j)
	PRINT USING F$; A(i, j);
      NEXT j
      PRINT
    NEXT i
    PRINT
    CLOSE #1
    eps = .0000000001#
    d1 = .00000001#
    d2 = .00000001#
    m = 600
  
  END SUB

  '***************************************************************
  '* Subroutine DECCRM determines the lower triangular matrix and*
  '* the upper triangukar matrix of Crout's decomposition of a   *
  '* given square real matrix, A.                                *
  '* ----------------------------------------------------------- *
  '* INPUTS:                                                     *
  '*         eps: required precision (double)                    *
  '*          n : size of matrix A (integer)                     *
  '*          A : input matrix (n x n)                           *
  '* OUTPUTS:                                                    *
  '*          it: flag, =0 if method does not apply              *
  '*                    =1 if method is ok.                      *
  '*           U: lower triangular matrix.                       *
  '*           V: upper triangular matrix.                       *
  '***************************************************************
  ' NOTE: A, U and V are shared with main program
  ' ---------------------------------------------
  SUB DECCRM (eps, n, it)

    IF ABS(A(1, 1)) < eps THEN
	  it = 0
    ELSE
      FOR i = 1 TO n
	    U(i, 1) = A(i, 1)
      NEXT i
      V(1, 1) = 1#
      FOR j = 2 TO n
	    V(1, j) = A(1, j) / U(1, 1)
      NEXT j
      it = 1: k = 2
      WHILE it <> 0 AND k <= n
	FOR i = 1 TO n
	  IF (i < k) THEN
	    U(i, k) = 0#
	  ELSE
	    s = 0#
	    FOR j = 1 TO k - 1
	      s = s + U(i, j) * V(j, k)
	    NEXT j
	    U(i, k) = A(i, k) - s
	  END IF
	NEXT i
	IF ABS(U(k, k)) < eps THEN
	  it = 0
	ELSE
	  FOR j = 1 TO n
	    IF j < k THEN
	      V(k, j) = 0#
	    ELSEIF j = k THEN
	      V(k, j) = 1#
	    ELSE
	      s = 0#
	      FOR i = 1 TO k - 1
		s = s + U(k, i) * V(i, j)
	      NEXT i
	      V(k, j) = A(k, j) / U(k, k)
	    END IF
	  NEXT j
	  k = k + 1
	END IF
      WEND
    END IF
  END SUB

  '******************************************************
  '* INPUTS:                                            *
  '* EPS : precision (Double)                           *
  '* D1  : precision d1 (Double)                        *
  '* D2  : precision d2 (Double)                        *
  '* M   : maximum number of iterations (integer)       *
  '* N   : order of matrix A (integer)                  *
  '* A   : input matrix to study (of MAT type)          *
  '* -------------------------------------------------- *
  '* OUTPUTS:                                           *
  '* IT  : -1 if no convergence is obtained (integer)   *
  '* R   : table of eigenvalues (of VEC type)           *
  '* VX  : table of eigenvectors (of MAT type)          *
  '******************************************************
  SUB EPMRI (eps, d1, d2, m, n, it, R())
    DIM X(NMAX), A1(NMAX, NMAX)

    FOR i = 1 TO n
      FOR j = 1 TO n
	A1(i, j) = A(i, j)
      NEXT j
    NEXT i

    VAMR eps, d2, m, n, it, R()

' restore A after VAMR
    FOR i = 1 TO n
      FOR j = 1 TO n
	A(i, j) = A1(i, j)
      NEXT j
    NEXT i

    j = 1
    WHILE it = 1 AND j <= n

      IIM eps, d1, m, n, it, R(j), X()
  
' restore A after IIM
      FOR i = 1 TO n
	FOR k = 1 TO n
	  A(i, k) = A1(i, k)
	NEXT k
	VX(i, j) = X(i)
      NEXT i
      j = j + 1
    WEND
	
  END SUB

  '************************************************************
  '* Procedure IIM calculates a real eigenvalue and the asso- *
  '* ciated eigenvector of a real square matrix the inverse   *
  '* iteration method.                                        *
  '* -------------------------------------------------------- *
  '* INPUTS:                                                  *
  '*         eps : absolute precision (double)                *
  '*         dta : relative precision (double)                *
  '*          m  : maximum number of iterations (integer)     *
  '*          n  : size of matrix A                           *
  '*          A  : input real square matrix (n x n)           *
  '* OUTPUTS:                                                 *
  '*          it : flag, =-1 if convergence is not obtained   *
  '*                     =1 if convergence is ok.             *
  '*        Gamma: starting value for the eigenvalue as input *
  '*               approximation of the eigenvalue with preci-*
  '*               sion dta in output.                        *
  '*          X1 : contains in output the associated eigen-   *
  '*               vector.                                    *
  '*                                                          *
  '************************************************************
  SUB IIM (eps, dta, m, n, it, gamma, X1())
    DIM W(NMAX), X0(NMAX)
    DIM LP(NMAX)

    FOR i = 1 TO n
      A(i, i) = A(i, i) - gamma
    NEXT i
    FOR k = 1 TO n - 1
      p0 = A(k, k): l0 = k
      FOR i = k + 1 TO n
	IF ABS(A(i, k)) > ABS(p0) THEN
	  p0 = A(i, k): l0 = i
	END IF
      NEXT i
      LP(k) = l0
      IF ABS(p0) < eps THEN
	p0 = eps: A(l0, k) = eps
      END IF
      IF l0 <> k THEN
	FOR j = k TO n
	  t0 = A(k, j): A(k, j) = A(l0, j): A(l0, j) = t0
	NEXT j
      END IF
      FOR i = k + 1 TO n
	A(i, k) = A(i, k) / p0
	FOR j = k + 1 TO n
	  A(i, j) = A(i, j) - A(i, k) * A(k, j)
	NEXT j
      NEXT i
    NEXT k

    IF ABS(A(n, n)) < eps THEN A(n, n) = eps
    FOR i = 1 TO n
      X0(i) = 1# / SQR(1# * i)
    NEXT i
    it = -1: l = 1
    WHILE it = -1 AND l <= m
      FOR i = 1 TO n
	W(i) = X0(i)
      NEXT i
      FOR k = 1 TO n - 1
	l0 = LP(k)
	IF l0 <> k THEN
	  t0 = W(k): W(k) = W(l0): W(l0) = t0
	END IF
	FOR i = k + 1 TO n
	  W(i) = W(i) - A(i, k) * W(k)
	NEXT i
      NEXT k
      X1(n) = W(n) / A(n, n)
      FOR i = n - 1 TO 1 STEP -1
	s = 0#
	FOR j = i + 1 TO n
	  s = s + A(i, j) * X1(j)
	NEXT j
	X1(i) = (W(i) - s) / A(i, i)
      NEXT i
      p0 = 0#
      FOR i = 1 TO n
	IF ABS(X1(i)) > ABS(p0) THEN p0 = X1(i)
      NEXT i
      FOR i = 1 TO n
	X1(i) = X1(i) / p0
      NEXT i
      phi = 0#
      FOR i = 1 TO n
	s = ABS(X1(i) - X0(i))
	IF s > phi THEN phi = s
      NEXT i

      IF phi < dta THEN
	gamma = gamma + 1# / p0
	it = 1
      ELSE
	FOR i = 1 TO n
	  X0(i) = X1(i)
	NEXT i
	l = l + 1
      END IF
    WEND
  END SUB  'IIM                                                            

  '*********************************************************
  '* Calculate the eigenvalues of a real square matrix by  *
  '* Rutishauser's Method.                                 *
  '* ----------------------------------------------------- *
  '* INPUTS:                                               *
  '*        eps: absolute precision (double)               *
  '*        dta: relative precision (double)               *
  '*         m : maximum number of iterations (integer)    *
  '*         n : size of matrix A (integer)                *
  '*         A : input real square matrix (n x n)          *
  '* OUTPUTS:                                              *
  '*         it: flag, =-1 if convergence is not obtained  *
  '*                   =1 if convergence is ok.            *
  '*         R : contains in output the n eigenvalues of A *
  '*                                                       *
  '*********************************************************
  SUB VAMR (eps, dta, m, n, it, R())

    t0 = 0#
    l = 1
    WHILE l <= m AND it <> 1
      FOR i = 1 TO n
	R(i) = A(i, i)
      NEXT i

      DECCRM eps, n, it

      IF it = 0 THEN
	FOR i = 1 TO n
	  A(i, i) = A(i, i) + 1#
	NEXT i
	t0 = t0 + 1#
      ELSE
	FOR i = 1 TO n
	  FOR j = 1 TO n
	    s = 0#
	    FOR k = 1 TO n
	      s = s + V(i, k) * U(k, j)
	    NEXT k
	    A(i, j) = s
	  NEXT j
	NEXT i
	phi = 0#
	FOR i = 1 TO n
	  s = ABS(A(i, i) - R(i))
	  IF s > phi THEN phi = s
	NEXT i
	IF phi < dta THEN
	  FOR i = 1 TO n
	    R(i) = A(i, i) - t0
	  NEXT i
	ELSE
	  l = l + 1
	  it = -1
	END IF
      END IF
    WEND
  END SUB