'************************************************************************
'* This program tests the procedure brown1 which is designed to solve   *
'* a nonlinear system of equations                                      *
'*             f0  (x0,x1,...,xn-1)  =  0                               *
'*             f1  (x0,x1,...,xn-1)  =  0                               *
'*             ...                                                      *
'*             fn-1(x0,x1,...,xn-1)  =  0                               *
'* using Brown's method.                                                *
'*                                                                      *
'* Scope of program:                                                    *
'* ================                                                     *
'* The program reads the input from screen and writes output onto file  *
'* browntst.lst.                                                        *
'*                                                                      *
'* After reading, the input is printed out for control purposes. Then   *
'* Brown's method is performed and, barring premature errors, the end   *
'* results are printed.                                                 *
'*                                                                      *
'* The system to solve is defined in Subroutine 100.                    *
'* -------------------------------------------------------------------- *
'* SAMPLE RUN:                                                          *
'*                                                                      *
'* Desired precision : 1e-8                                             *
'* Protocole (0 or 1): 0                                                *
'* Maximum number of iterations : 25                                    *
'* Components of starting vector:                                       *
'* 1                                                                    *
'* 1                                                                    *
'* 1                                                                    *
'* 1                                                                    *
'*                                                                      *
'* The output file browntst.lst contains:                               *
'*                                                                      *
'* -------------------------------------------------------------------- *
'* Brown's method for nonlinear systems of equations                    *
'* -------------------------------------------------------------------- *
'* Example nø 3                                                         *
'* System to be solved:                                                 *
'* 10*x0 + x1 + x2 + x3 - 20 + sin(x0)^2 + cos(x1)^2 = 0                *
'* x0 + 20*x1 + x2 + x3 - 48 + 1/(x[0]^6)            = 0                *
'* (x0+x1)^2 + 30*x2 + x3 - 97 + Ln(x0) + Ln(x1+x2)  = 0                *
'*  x0 + x1 + x2 + 40*x3 -166 + x0*x0                = 0                *
'* Starting vector:                                                     *
'*  1   1   1   1                                                       *
'* Error bound =  .00000001                                             *
'* Maximal number of iterations =  25                                   *
'*                                                                      *
'* Intermediate results are not kept.                                   *
'* Solution vector:                                                     *
'*  1.040648  1.972398  2.745049  3.978974                              *
'* Number of iterations:  4                                             *
'* -------------------------------------------------------------------- *
'* Ref.: "Numerical algorithms with C by Gisela Engeln-Muellges and     *
'*        Franck Uhlig, Springer-Verlag, 1996" [BIBLI 11].              *
'*                                                                      *
'*                                  Basic Version By J-P Moreau, Paris. *
'*                                          (www.jpmoreau.fr)           *
'************************************************************************
'Program Test_Brown
DEFINT I-N
DEFDBL A-H, O-Z

OPTION BASE 0

ISIZE = 24
XMACHEPS = 2E-16

DIM ihf(ISIZE, ISIZE + 1)
DIM dquot(ISIZE), rslin(ISIZE)
DIM hf(ISIZE, ISIZE)

'integer n,            ! size of system
'        maxit,        ! maximal number of iterations
'        ierror,       ! error code of subroutine 2000
'        itanz,        ! number of iterations performed
'        i,            ! Loop variable
'        iprot         ! protocole flag (0=no intermediate results)
'double  eps,          ! desired accuracy
DIM x0(ISIZE)          ' [0..n-1] starting vector
DIM x1(ISIZE)          ' [0..n-1] approximate solution

' -------------------- read input ----------------------------------
  CLS
  PRINT
  n = 4                ' size of system
  INPUT " Desired precision : ", eps
  PRINT
  INPUT " Protocole (0 or 1): ", iprot
  PRINT
  INPUT " Maximal number of iterations : ", maxit
  PRINT
  PRINT " Componants of starting vector: "
  FOR i = 0 TO n - 1
    INPUT "", x0(i)
  NEXT i
  PRINT
' ------------ print input for checking purposes -------------------
  OPEN "browntst.lst" FOR OUTPUT AS #1
  PRINT #1, " ----------------------------------------------------"
  PRINT #1, "  Brown's method for nonlinear systems of equations "
  PRINT #1, " ----------------------------------------------------"
  PRINT #1, " Example nø 3"
  PRINT #1, " System to be solved:"
  PRINT #1, " 10*x0 + x1 + x2 + x3 - 20 + sin(x0)^2 + cos(x1)^2 = 0"
  PRINT #1, " x0 + 20*x1 + x2 + x3 - 48 + 1/(x[0]^6)            = 0"
  PRINT #1, " (x0+x1)^2 + 30*x2 + x3 - 97 + Ln(x0) + Ln(x1+x2)  = 0"
  PRINT #1, " x0 + x1 + x2 + 40*x3 -166 + x0*x0                 = 0"
  PRINT #1, " Starting vector:"
  FOR i = 0 TO n - 1
    PRINT #1, " "; x0(i);
  NEXT i
  PRINT #1,
  PRINT #1, " Error bound = "; eps
  PRINT #1, " Maximal number of iterations = "; maxit
  PRINT #1,
  IF iprot <> 0 THEN
    PRINT #1, " Intermediate results are saved."
  ELSE
    PRINT #1, " Intermediate results are not kept."
  END IF
' ------------ solve nonlinear system ------------------------------
  GOSUB 2000 'call Brown1(n, x0, eps, iprot, maxit, x1, itanz, irc)

  ierror = irc
  IF ierror = 0 THEN
	GOTO 10
  ELSEIF ierror = 1 THEN
	PRINT " brown: too many steps"
  ELSEIF ierror = 2 THEN
	PRINT " brown: linearized system singular"
  ELSEIF ierror = 3 THEN
	PRINT " brown: lack of memory"
  ELSEIF ierror = 4 THEN
	PRINT " brown: wrong input parameter: n < 1 or maxit < 1"
  ELSEIF ierror = 5 THEN
	PRINT " brown: error calling function"
  END IF

' ---------------------- print solution -----------------------------
10 PRINT #1, " Solution vector:"
  FOR i = 0 TO n - 1
    PRINT #1, USING "###.######"; x1(i);
  NEXT i
  PRINT #1,
  PRINT #1, " Number of iterations: "; itanz
  PRINT #1, " ----------------------------------------------------"
  CLOSE #1
  
  PRINT " Results in file browntst.lst."
  PRINT

END 'of main program

100 'Subroutine func3(k,x1,y)
'************************************************************************
'*                         Test example 3                               *
'************************************************************************
  IF k = 0 THEN
    y = 10# * x1(0) + x1(1) + x1(2) + x1(3) - 20# + SIN(x1(0)) ^ 2 + COS(x1(1)) ^ 2
  ELSEIF k = 1 THEN
    y = x1(0) + 20# * x1(1) + x1(2) + x1(3) - 48# + 1# / (x1(0) ^ 6)
  ELSEIF k = 2 THEN
    y = (x1(0) + x1(1)) ^ 2 + 30# * x1(2) + x1(3) - 97# + LOG(x1(0)) + LOG(x1(1) + x1(2))
  ELSEIF k = 3 THEN
    y = x1(0) + x1(1) + x1(2) + 40# * x1(3) - 166# + x1(0) ^ 2
  END IF
RETURN

' -----------------------------------------------------------------------
500 'Subroutine subst (n, k, ihf, hf, rslin, x1)
'************************************************************************
'* Solve a linear system of equations                                   *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* n      size of system                                                *
'* k      Index for coordinates: 0,...,n-1                              *
'* ihf    [0..n-1,0..n] matrix, register of row and column interchanges *
'* hf     [0..n-1,0..n-1] matrix, the system matrix                     *
'* rslin  [0..n-1] right nahd side vector                               *
'*                                                                      *
'* Output parameter:                                                    *
'* =================                                                    *
'* x1     [0..n-1] solution vector                                      *
'*                                                                      *
'************************************************************************
'double sum              ! aux variable for finding x1[kmax]
'integer km,             ! Loop counter
'       kmax,            ! original row index
'       jsub,            ! original column index
'       j                ! Loop counter
  FOR km = k TO 1 STEP -1
    kmax = ihf(km - 1, n)
    sum = 0#
    FOR j = km TO n - 1
      jsub = ihf(km, j)
      sum = sum + hf(km - 1, jsub) * x1(jsub)
    NEXT j
    x1(kmax) = sum + rslin(km - 1)
  NEXT km
RETURN

' -----------------------------------------------------------------------
1000 'Subroutine iter4(n, epsm, xalt, x1, ising, irc)
'************************************************************************
'* Compute one approximation via  Brown's method.                       *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* n     number of equations                                            *
'* epsm  machine constant                                               *
'* xalt  previous iterate                                               *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* x1    [0..n-1] vector, the new iterate for the zero                  *
'* ising error code, (1=matrix is singular)                             *
'*                                                                      *
'* Return code irc:                                                     *
'* ===============                                                      *
'* = 0: all ok                                                          *
'* = 1: error in evaluating function (not used here)                    *
'*                                                                      *
'************************************************************************
'Labels:1010, 1020
'integer i,         ! Loop counters
'        j,
'        k,         ! index of current function
'        ianzahl,   ! counter for those difference quotients that are 0
'        itemp,     ! aux variable
'        kmax,      ! Index of largest difference quotient
'        jsub       ! aux variable for an index in ihf
'double  hold,      ! aux variable for x1[temp]
'        h,         ! Step size for difference quotient in direction of x1[temp]
'        faktor,    ! quotient  h / hold
'        dermax,    ! maximal difference quotient magnitude
'        sum,       ! sum in rslin[k]
'        f,         ! value of function k at x1
'        fplus      ! value of function k at x1[0],..., x1[temp]+h,...,x1[n-1]
' --------------------- initialize variables --------------------------
  itemp = 0: kmax = 0
  FOR j = 0 TO n - 1
    ihf(0, j) = j
    x1(j) = xalt(j)
  NEXT j

' ----------- linearize the kth coordinate function -------------------
  FOR k = 0 TO n - 1

    ianzahl = 0
    faktor = .001#

    FOR j = 0 TO 2
      IF k > 0 THEN
	GOSUB 500 'call subst(n, k, ihf, hf, rslin, x1)
      END IF
      
      GOSUB 100  'call func3(k, x1, f)
      f = y

'     ---- find ith diskretization size and difference quotient ----
      FOR i = k TO n - 1
	itemp = ihf(k, i)
	hold = x1(itemp)
	h = faktor * hold
	IF ABS(h) <= epsm THEN
	  h = .001#
	END IF
	x1(itemp) = hold + h
	IF k > 0 THEN
	  GOSUB 500 'call subst(n, k, ihf, hf, rslin, x1)
	END IF
	
	GOSUB 100 'call func3(k, x1, fplus)
	fplus = y

	x1(itemp) = hold

	dquot(itemp) = (fplus - f) / h

	IF ABS(dquot(itemp)) <= epsm THEN
	  ianzahl = ianzahl + 1
	ELSEIF ABS(f / dquot(itemp)) >= 1E+20 THEN
	  ianzahl = ianzahl + 1
	END IF
      NEXT i

      IF ianzahl < n - k THEN
	ising = 0
	GOTO 1010
      ELSE
	ising = 1
	faktor = faktor * 10#
	ianzahl = 0
      END IF
    NEXT j

1010 IF ising = 1 THEN
      GOTO 1020
    ELSEIF k < n - 1 THEN
      kmax = ihf(k, k)

'     --- find largest magnitude difference quotient ------------
      dermax = ABS(dquot(kmax))
      FOR i = k + 1 TO n - 1
	jsub = ihf(k, i)
	IF ABS(dquot(jsub)) < dermax THEN
	  ihf(k + 1, i) = jsub
	ELSE
	  ihf(k + 1, i) = kmax
	  kmax = jsub
	END IF
      NEXT i
      IF ABS(dquot(kmax)) <= epsm THEN
	ising = 1
      END IF

      ihf(k, n) = kmax

      IF ising = 1 THEN
	GOTO 1020
      ELSE
'       ---------- solve kth equation for  xmax ------------------
	sum = 0#
	FOR j = k + 1 TO n - 1
	  jsub = ihf(k + 1, j)
	  hf(k, jsub) = -dquot(jsub) / dquot(kmax)
	  sum = sum + dquot(jsub) * x1(jsub)
	NEXT j
	rslin(k) = (sum - f) / dquot(kmax) + x1(kmax)
      END IF
    ELSE
'     ----- solve (n-1)th coordinate function via discrete -------
'     ----- Newton method for one variable -----------------------
      IF ABS(dquot(itemp)) <= epsm THEN
	ising = 1
	  ELSE
	kmax = itemp
	rslin(k) = -f / dquot(kmax) + x1(kmax)
      END IF
    END IF
1020 NEXT k

  IF ising = 0 THEN
    x1(kmax) = rslin(n - 1)                  ' compute approximation
    IF n > 1 THEN                            ' by back substitution
      GOSUB 500 'call subst(n, n - 1, ihf, hf, rslin, x1)
    END IF
  END IF
  irc = 0
RETURN

' ---------------------------------------------------------------------
' Brown's method for nonlinear systems of equations
2000 'Subroutine brown1(n, x0, eps, iprot, maxit, x1, itanz, irc)
'integer   n              !number of equations
'double    x0(0:SIZE)     !Starting value for iteration
'double    eps            !error bound
'integer   prot,          !Protocole switch
'          maxit          !maximal number of steps
'double    x1(0:SIZE)     !solution vector
'integer   itanz,         !actual steps performed
'          irc            !error code
'************************************************************************
'* find a zero of a nonlinear system of n equations in n unknown using  *
'* Brown's method                                                       *
'*                                                                      *
'* One step of Brown's method is performed by calling iter4().          *
'* The iteration is continued until the admissable number of iterations *
'* is reached or one of the following break-off criteria is satisfied : *
'* 1. the relative change between two successive iterates is less than  *
'*    eps.                                                              *
'* 2. The function value has magnitude less than XMACHEPS at the new x  *
'* 3. The desired accuracy has been reached.                            *
'*                                                                      *
'* If specified in iprot, intermediate results are tabulated, see iprot *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* n      size of system                                                *
'* x0     [0..n-1] starting vector for the iteration                    *
'* eps    desired accuracy; if specified as <= 0, we set                *
'*        eps = 0.8 * MACH_EPS.                                         *
'* iprot  Protocole flag. If =1, we tabulate the differenz of succes-   *
'*        sive iterates, the last iterate and the functional value      *
'*        there after each iteration in the standard output file.       *
'*        If =0, there is no output.                                    *
'* maxit  Maximal number of iterations                                  *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* x1     [0..n-1] vector, the approximate solution for the zero        *
'* itanz  number of iterations executed                                 *
'*                                                                      *
'* Return value irc:                                                    *
'* ================                                                     *
'* error code.                                                          *
'* = 0: successsful iteration                                           *
'* = 1: desired accuracy not achieved after maxit iterations            *
'* = 2: system matrix singular                                          *
'* = 3: lack of memory (not used here)                                  *
'* = 4: wrong input : n < 1 or maxit < 1                                *
'* = 5: error in evaluating function (not used here)                    *
'*                                                                      *
'************************************************************************
'Labels: 2030,2035,2040,2050
'integer i,           ! Loop variable
'        m,           ! iterations counter
'        krit         ! list of causes for stop of program:
'                     ! = 0: too many steps                       
'                     ! = 1: step size too small                  
'                     ! = 2: function value sufficiently small    
'                     ! = 3: desired accuracy reached
'integer ising        ! Flag, indicating whether the linearized system
'                     ! is singular (ising=1, else 0)
'double  relf, fwert, ! aux variables
'        delta0,      ! maximum norm of previous vector
'        delta1,      ! maximum norm of current vector
'        epsm,        ! machine constant
'        xalt(0:SIZE) ! [0..n-1] vector used to store old iterate

  relf = 0#
' ---------------- catch wrong input ---------------------------------
  IF n < 1 OR maxit < 1 THEN
    irc = 4
    GOTO 2050  'return
  END IF

  IF eps < XMACHEPS THEN           ' desired accuracy rediculously
    eps = 8# * XMACHEPS            ' small ? correct
  END IF

  epsm = 100# * XMACHEPS           ' initialize variables
  ising = 0
  krit = 0
  delta0 = .01#
  FOR ii = 0 TO n - 1
    xalt(ii) = x0(ii)
  NEXT ii

  IF iprot <> 0 THEN
    PRINT
  END IF

  FOR m = 1 TO maxit                ' start iteration

    GOSUB 1000  'call iter4(n, epsm, xalt, x1, ising, irc)
    IF irc <> 0 THEN
      irc = 5
      GOTO 2050  'return
    END IF

'    ----------- if desired document each step ------------------------
    IF iprot <> 0 THEN
      PRINT " Iteration step: "; m
      IF ising = 0 THEN
	PRINT " Difference         Component       Approximation"
	PRINT " function value"
	FOR i = 0 TO n - 1
	  k = i
	  GOSUB 100 'call func3(i,x1,fwert)
	  fwert = y
	  PRINT x1(i) - xalt(i), i, x1(i)
	  PRINT fwert
	  INPUT "", R$
	NEXT i
      ELSE
	PRINT " Jacobi matrix singular!"
      END IF
    END IF

    IF ising = 0 THEN

'     ----------------- test break-off criteria ----------------------
      FOR i = 0 TO n - 1            ' test relative change in new iterate
	relf = (x1(i) - xalt(i)) / (xalt(i) + eps)
	IF ABS(relf) >= eps THEN
	  GOTO 2030
	END IF
      NEXT i
2030  IF ABS(relf) < eps THEN
	krit = 1                           ' step size too small
	GOTO 2040
      END IF

      FOR i = 0 TO n - 1                   ' check function value
	k = i
	GOSUB 100 'call func3(i, x1, fwert)
	fwert = y
	IF ABS(fwert) > epsm THEN
	  GOTO 2035
	END IF
      NEXT i

2035  IF ABS(fwert) <= epsm THEN
	krit = 2                           ' function value small enough
	GOTO 2040
      END IF

      delta1 = ABS(x1(0) - xalt(0))             ' compare with desired
      FOR i = 1 TO n - 1                        ' accuracy
	IF delta1 < ABS(x1(i) - xalt(i)) THEN
	  delta1 = ABS(x1(i) - xalt(i))
	END IF
      NEXT i
      IF delta1 <= .001# THEN
	IF delta0 <= delta1 THEN
	  krit = 3                          ' desired accuracy reached
	  GOTO 2040
	END IF
      END IF
      
      delta0 = delta1

      IF m < maxit THEN
	FOR ii = 0 TO n - 1
	  xalt(ii) = x1(ii)
	NEXT ii
      END IF
    ELSE
      GOTO 2040
    END IF
  NEXT m

2040 itanz = m

  IF ising = 1 THEN
    irc = 2
  ELSEIF krit = 0 THEN
    irc = 1
  ELSE
    irc = 0
  END IF
2050 RETURN

' end of file browntst.bas