'***************************************************************
'*  Calculate the greatest eigenvalue of a real square matrix  *
'*     and the associated eigenvector by the power method      *
'* ----------------------------------------------------------- *
'* Ref.: "Algèbre - Algorithmes et programmes en Pascal By     *
'*        Jean-Louis Jardrin, Dunod Editor, Paris, 1988"       *
'*        [BIBLI 10].                                          *
'* ----------------------------------------------------------- *
'* SAMPLE RUN:                                                 *
'*                                                             *
'* ----------------------------------------------              *
'*  Calculate the greatest eigenvalue of a real                *
'*  square matrix and the associated eigenvector               *
'*  by the power method.                                       *
'* ----------------------------------------------              *
'*                                                             *
'*  Size of matrix (maximum 25): 4                             *
'*                                                             *
'*   Line  1                                                   *
'*    Element  1: 1                                            *
'*    Element  2: 3                                            *
'*    Element  3: 0                                            *
'*    Element  4: 0                                            *
'*                                                             *
'*   Line  2                                                   *
'*    Element  1: 4                                            *
'*    Element  2: 2                                            *
'*    Element  3: 0                                            *
'*    Element  4: 0                                            *
'*                                                             *
'*   Line  3                                                   *
'*    Element  1: 1                                            *
'*    Element  2: -1                                           *
'*    Element  3: 5                                            *
'*    Element  4: -3                                           *
'*                                                             *
'*   Line  4                                                   *
'*    Element  1: 2                                            *
'*    Element  2: 0                                            *
'*    Element  3: 4                                            *
'*    Element  4: -2                                           *
'*                                                             *
'*  Precision: 1e-10                                           *
'*  Epsilon  : 1e-10                                           *
'*  Maximum number of iterations: 27                           *
'*                                                             *
'*    Eigenvalue: 5.000000000086734                            *
'*                                                             *
'*    Eigenvector:                                             *
'*    .7499999999913266                                        *
'*    1                                                        *
'*   -.5208333332778221                                        *
'*   -8.333333327999046D-02                                    *          
'*                                                             *
'*                 English Basic version by J-P Moreau, Paris. *
'*                            (www.jpmoreau.fr)                *
'***************************************************************
'  Exact values are: gamma = 5
'                    eigenvector = (1/48)(36,48,-25,-4)
'---------------------------------------------------------------                               
'Program Test_PWM
DEFDBL A-H, O-Z
DEFINT I-N

NMAX = 25

DIM A(NMAX, NMAX), X(NMAX)

  CLS
  PRINT " ----------------------------------------------"
  PRINT "  Calculate the greatest eigenvalue of a real  "
  PRINT "  square matrix and the associated eigenvector "
  PRINT "  by the power method.                         "
  PRINT " ----------------------------------------------"

  GOSUB 500  'call Read_data(n,A,dta,eps,m)

  GOSUB 1000 'call PWM(eps,dta,m,n,A,it,gamma,X)

  IF it + 1 = 0 THEN
    PRINT
    PRINT " No convergence !"
  ELSEIF it + 1 = 1 THEN
    PRINT
    PRINT "  Method does not apply."
  ELSE
    PRINT
    PRINT "  Eigenvalue: "; gamma
    PRINT
    PRINT "  Eigenvector:"
    FOR i = 1 TO n
	  PRINT "  "; X(i)
    NEXT i
  END IF
  PRINT
  PRINT " Number of iterations: "; l
  PRINT

END

500 'Subroutine Read_data(n,A,dta,eps,m)
  PRINT " Size of matrix (maximum "; NMAX; "): "; : INPUT "", n
  FOR i = 1 TO n
    PRINT
    PRINT "  Line "; i
    FOR j = 1 TO n
      PRINT "   Element "; j; : INPUT "", A(i, j)
    NEXT j
  NEXT i
  PRINT
  INPUT " Precision: ", dta
  INPUT " Epsilon  : ", eps
  INPUT " Maximum number of iterations: ", m
RETURN


'************************************************************
'* calculate greatest eigenvalue and associated eigenvector *
'* by the power method                                      *
'* -------------------------------------------------------- *
'* INPUTS:                                                  *
'*         eps   : smallest number in double precision      *
'*         dta   : required precision                       *
'*         m     : maximum number of iterations             *
'*         n     : size of real square matrix A(n,n)        *
'*         A     : real square matrix A(n,n)                *
'* OUTPUTS:                                                 *
'*         it    : error indicator: -1=no convergence,      *
'*                 0=method cannot be applied,              *
'*                 1=convergence ok.                        *
'*         gamma : greatest eigenvalue (in absolute value)  *
'*                 of input matrix A(n,n)                   *
'*         X1    : associated eigenvector                   *
'************************************************************
1000 'Subroutine PWM(eps,dta,m,n,A,it,gamma,X)
DIM X0(NMAX)
  FOR i = 1 TO n
    X0(i) = 1# / SQR(i)
  NEXT i
  it = -1: l = 1
1010 gamma = 0#
    FOR i = 1 TO n
      X(i) = 0#
      FOR j = 1 TO n
	X(i) = X(i) + A(i, j) * X0(j)
      NEXT j
      IF ABS(X(i)) > ABS(gamma) THEN gamma = X(i)
	NEXT i
    IF ABS(gamma) < eps THEN
      it = 0
    ELSE
      FOR i = 1 TO n
	X(i) = X(i) / gamma
      NEXT i
      phi = 0#
      FOR i = 1 TO n
	s = ABS(X(i) - X0(i))
	IF s > phi THEN phi = s
      NEXT i
      IF phi < dta THEN
	it = 1
      ELSE
	FOR i = 1 TO n
	  X0(i) = X(i)
	NEXT i
	l = l + 1
      END IF
    END IF
  IF it = -1 AND l <= m THEN GOTO 1010
RETURN

'end of file tpwm.bas