{ ------------------------ UNIT brown.pas ------------------------------

************************************************************************
*                                                                      *
* Brown's method for nonlinear systems                                 *
* -------------------------------------                                *
*                                                                      *
* Programing language: Pascal                                          *
* Compiler:            Borland TPW 1.5                                 *
* Computer:            IBM Compatible PC, 400 MHz                      *
* Literature:          Brown, K. M.:                                   *
*                      A quadratically convergent Newton-like method   *
*                      based upon Gaussian elimination                 *
*                      SIAM J. Numer. Anal. Vol. 6 , (1969), p. 560    *
* Source:              QuickBASIC program                              *
* Author:              Johannes Karfusehr (FORTRAN)                    *
* Adaption:            Juergen Dietel, Computer Center, RWTH Aachen    *
*                      for C version                                   *
*                      J-P Moreau, Paris for this Pascal version       *
* Date:                9.16.2002                                       *
*                      (www.jpmoreau.fr)                               *
* -------------------------------------------------------------------- *
* Ref.: "Numerical Algorithms with C By G. Engeln-Mueller and F. Uhlig,*
*        Springer-Verlag, 1996" [BIBLI 11].                            *
***********************************************************************}
UNIT Brown;

INTERFACE

Uses Basis;

Const
      EIGHT= 8.0;
      HALF = 0.5;
      ONE  = 1.0;
      PI   = 3.1415926535;
      SIZE = 24;
      TEN  = 10.0;
      TWO  = 2.0;
      ZERO = 0.0;

Type  REAL = Double;
      MAT  = Array[0..SIZE,0..SIZE] of REAL;
      IMAT = Array[0..SIZE,0..SIZE+1] of INTEGER;
      VEC  = Array[0..SIZE] of REAL;

Var
      ihf: IMAT;
      dquot, rslin: VEC;
      hf: MAT;

Procedure brown1  { Brown's method for nonlinear systems of equations }
         (
          n:      integer;   { number of equations ...................}
          x0:     VEC;       { Starting value for iteration ..........}
          eps:    REAL;      { error bound ...........................}
          prot:   Integer;   { Protokol switch .......................}
          maxit:  integer;   { maximal number of steps ...............}
      VAR x1:     VEC;       { solution ..............................}
      VAR itanz:  integer;   { actual steps performed ................}
      VAR rc:     Integer    { error code ............................}
         );


IMPLEMENTATION


Function func0(k:Integer; x:VEC): REAL;
{***********************************************************************
*                         Test example 0                               *
***********************************************************************}
Begin
  Case k of
    0:  func0 := sqr(x[0]) - x[1] - ONE;
    1:  func0 := sqr(x[0] - TWO) + sqr(x[1] - HALF) - ONE
  End
End;

Function func1(k:Integer; x:VEC): REAL;
{***********************************************************************
*                         Test example 1                               *
***********************************************************************}
Begin
  Case k of
    0:  func1 := x[0] - x[1] + sin(x[0]) + cos(x[1]) - PI/EIGHT;
    1:  func1 := sqr(x[0]) - 4*sqr(x[1]) - Ln(sqr(x[2]));
    2:  func1 := x[2]-1 + 2*x[0] - 4*x[1]
  End
End;

Function func(k:Integer; x:VEC): REAL;
{***********************************************************************
*                         Test example 3                               *
***********************************************************************}
Begin
  Case k of
    0: func := 10*x[0] + x[1] + x[2] + x[3] -20 + sin(x[0])*sin(x[0]) + cos(x[1])*cos(x[1]);
    1: func := x[0] + 20*x[1] + x[2] + x[3] -48 + 1/(x[0]*x[0]*x[0]*x[0]*x[0]*x[0]);
    2: func := (x[0]+x[1])*(x[0]+x[1]) + 30*x[2] + x[3] -97 + Ln(x[0]) + Ln(x[1]+x[2]);
    3: func := x[0] + x[1] + x[2] + 40*x[3] -166 + x[0]*x[0]
  End
End;

{ -------------------------------------------------------------------- } 
Procedure subst (
                  n:integer;
                  k:integer;
                  ihf:IMAT;
                  hf:MAT;
                  rslin:VEC;
              VAR x1:VEC
                 );
{***********************************************************************
* 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                                      *
*                                                                      *
* Constants in use :  REAL, ZERO                                       *
* ==================                                                   *
*                                                                      *
***********************************************************************}
Var
       sum: REAL;        { aux variable for finding x1[kmax]         }
       km,               { Loop counter                              }
       kmax,             { original row index                        }
       jsub,             { original column index                     }
       j:Integer;        { Loop counter                              }
Begin

  for km := k Downto 1 do
  begin
    kmax := ihf[km - 1,n];
    sum:=ZERO;
    for j := km to n-1 do
    begin
      jsub :=  ihf[km,j];
      sum  := sum + hf[km - 1,jsub] * x1[jsub]
    end;
    x1[kmax] := sum + rslin[km - 1]
  end
End;


{ -------------------------------------------------------------------- }
Procedure iter4
                (
                 n: integer;
                 epsm: REAL;
                 xalt: VEC;
             VAR x1: VEC;
             VAR sing: boolean;
             VAR rc: Integer
                );
{***********************************************************************
* Compute one approximation via  Brown's method.                       *
*                                                                      *
* Input parameters:                                                    *
* =================                                                    *
* fkt   Function, which evaluates the right hand side. It has          *
*       components from  0 to n - 1.                                   *
* n     number of equations                                            *
* epsm  machine constant                                               *
* xalt  previous iterate                                               *
*                                                                      *
* Output parameters:                                                   *
* ==================                                                   *
* x1    [0..n-1] vector, the new iterate for the zero                  *
* sing  error code, matrix is singular                                 *
*                                                                      *
* Return code rc:                                                      *
* ==============                                                       *
* = 0: all ok                                                          *
* = 1: error in evaluating func                                        *
*                                                                      *
* global names used:                                                   *
* ==================                                                   *
* ihf, rslin, dquot, hf, nlgls, REAL, boolean, subst, ABS, FALSE,      *
* TRUE, ZERO                                                           *
***********************************************************************}
Label 10, 20;
Var
       i,          { Loop counters                                    }
       j,
       k,          { index of current function                        }
       anzahl,     { counter for those difference quotients that are 0}
       temp,       { aux variable                                     }
       kmax,       { Index of largest difference quotient             }
       jsub:Integer; { aux variable for an index in ihf               }
       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:REAL; { value of function k at                           }
                   { (x1[0],..., x1[temp]+h,...,x1[n-1])              }
Begin

  { --------------------- initialize variables ----------------------- }
  temp:=0; kmax:=0;
  for j := 0 to n-1 do
  begin
    ihf[0,j] := j;
    for i:=0 to n-1 do x1[n]:=xalt[n]
  end;

  { ----------- linearize the kth coordinate function ---------------- }
  for k := 0 to n-1 do
  begin

    anzahl := 0;
    faktor := 0.001;

    for j := 0 to 2 do
    begin
      if k > 0 then
        subst(n, k, ihf, hf, rslin, x1);
      f:=func(k, x1);

      { ---- find ith diskretization size and difference quotient ---- }
      for i := k to n-1 do
      begin
        temp := ihf[k,i];
        hold := x1[temp];
        h    := faktor * hold;
        if ABS(h) <= epsm then
          h := 0.001;
        x1[temp] := hold + h;
        if k > 0 then
          subst(n, k, ihf, hf, rslin, x1);
        fplus:=func(k, x1);
        x1[temp] := hold;

        dquot[temp] := (fplus - f) / h;

        if ABS(dquot[temp]) <= epsm then
          Inc(anzahl)
        else if ABS(f / dquot[temp]) >= 1e20 then
          Inc(anzahl)
      end;

      if anzahl < n - k then
      begin
        sing := FALSE;
        goto 10;
      end
      else
      begin
        sing  :=  TRUE;
        faktor := faktor * TEN;
        anzahl :=  0
      end
    end;


10: if sing then goto 20
    else if k < n - 1 then
    begin
      kmax := ihf[k,k];

      { --- find largest magnitude difference quotient ------------ }
      dermax := ABS(dquot[kmax]);
      for i := k + 1 to n-1 do
      begin
        jsub := ihf[k,i];
        if ABS(dquot[jsub]) < dermax then
          ihf[k + 1,i] := jsub
        else
        begin
          ihf[k + 1,i] := kmax;
          kmax         := jsub
        end
      end;
      if ABS(dquot[kmax]) <= epsm then
        sing := TRUE;

      ihf[k,n] := kmax;


      if sing then goto 20
      else
      begin
        { ---------- solve kth equation for  xmax ------------------ }
        sum := ZERO;
        for j := k + 1 to n-1 do
        begin
          jsub       := ihf[k + 1,j];
          hf[k,jsub] := -dquot[jsub] / dquot[kmax];
          sum        := sum + dquot[jsub] * x1[jsub]
        end;
        rslin[k] := (sum - f) / dquot[kmax] + x1[kmax]
      end
    end
    else
    begin
      { ----- solve (n-1)th coordinate function via discrete ------- 
        ----- Newton method for one variable ----------------------- }
      if ABS(dquot[temp]) <= epsm then
        sing := TRUE
      else
      begin
        kmax     := temp;
        rslin[k] := -f / dquot[kmax] + x1[kmax]
      end
    end;
20: end;

  if Not sing then
  begin
    x1[kmax] := rslin[n - 1];                 { compute approximation }
    if n > 1 then                            { by back substitution  }
      subst(n, n - 1, ihf, hf, rslin, x1)
  end;
  rc:=0
End;

{ ------------------------------------------------------------------- }
Procedure brown1  { Brown's method for nonlinear systems of equations }
         (
          n:      Integer;   { number of equations ...................}
          x0:     VEC;       { Starting value for iteration ..........}
          eps:    REAL;      { error bound ...........................}
          prot:   Integer;   { Protokol switch .......................}
          maxit:  integer;   { maximal number of steps ...............}
      VAR x1:     VEC;       { solution ..............................}
      VAR itanz:  integer;   { actual steps performed ................}
      VAR rc:     Integer    { 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 mgnitude less than MACH_EPS at the new x   *
* 3. The desired accuracy has been reached.                            *
*                                                                      *
* If specified in prot, intermediate results are tabulated, see prot   *
*                                                                      *
* Input parameters:                                                    *
* =================                                                    *
* fkt    Function, whichdescribes the right hand side of the system;   *
*        its components are numbered 0 to n-1.                         *
* 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.                                         *
* prot   Protocol flag. If TRUE we tabulate the differenz of successive*
*        iterates, the last iterate and the functional value there     *
*        after each iteration in the standard output file. If FALSE,   *
*        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 rc:                                                     *
* ===============                                                      *
* error code.                                                          *
* = 0: successsful iteration                                           *
* = 1: desired accuracy not achieved after maxit iterations            *
* = 2: system matrix singular                                          *
* = 3: lack of memory                                                  *
* = 4: wrong input : fkt = NULL or n < 1 or maxit < 1                  *
* = 5: error in evaluating  fkt()                                      *
*                                                                      *
* global names used:                                                   *
* ==================                                                   *
* ihf, rslin, dquot, hf, nlgls, REAL, boolean, MACH_EPS, FALSE,        *
* iter4, subst, IMATRIX, MATRIX, writeln, ABS, EIGHT                   *
***********************************************************************}
Label 30,35,40,50;
Var
    i,             { Loop variable                                     }
    m,             { iterations counter                                }
    krit:Integer;  { 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                     }
    sing:boolean;  { Flag, indicating whether the linearized system    
                     is singular                                       }
    relf,          { aux variables                                     }
    fwert,         {                                                   }
    delta0,        { maximum norm of previous vector                   }
    delta1,        { maximum norm of current vector                    }
    epsm:REAL;     { machine constant                                  }
    xalt:VEC;      { [0..n-1] vector used to store old iterate         }
Begin
  relf:=ZERO; 
  { ---------------- catch wrong input ------------------------------- }
  if (n < 1) or (maxit < 1) then
  begin
    rc := 4;
    goto 50  {exit}
  end;

  if eps < MACH_EPS then           { desired accuracy rediculously }
    eps := EIGHT * MACH_EPS;       { small ? correct }

  epsm   := 100*MACH_EPS;          { initialize variables }
  sing   := FALSE;
  krit   := 0;
  delta0 := 0.01;
  for i:=0 to n-1 do xalt[n]:=x0[n];

  if prot<>0 then writeln;

  for m := 1 to maxit do           { start iteration }
  begin

    iter4(n, epsm, xalt, x1, sing, rc);
    if rc<>0 then
    begin
      rc:=5;
      goto 50 {exit}
    end;

    { ----------- if desired document each step ---------------------- }
    if prot<>0 then
    begin
      writeln(m:3,' Iteration step');
      if Not sing then
      begin
        writeln('      Difference         Component       Approximation             Function value');
        for i := 0 to n-1 do
        begin
          fwert:=func(i,x1);
          writeln(x1[i] - xalt[i], i, x1[i], fwert)
        end
      end
      else
        writeln(' Jacobi matrix singular!')
    end;

    if Not sing then
    begin

      { ----------------- test break-off criteria -------------------- }
      for i := 0 to n-1 do       { test relative change in new iterate }
      begin
        relf := (x1[i] - xalt[i]) / (xalt[i] + eps);
        if ABS(relf) >= eps then goto 30
      end;
30:   if ABS(relf) < eps then
      begin
        krit := 1;                         { step size too small }
        goto 40
      end;

      for i := 0 to n-1 do                 { check function value }
      begin
        fwert:=func(i, x1);
        if ABS(fwert) > epsm then
          goto 35
      end;

35:   if ABS(fwert) <= epsm then
      begin
        krit := 2;                       { function value small enough }
        goto 40
      end;

      delta1 := ABS(x1[0] - xalt[0]);            { compare with desired }
      for i := 1 to n-1 do                      { accuracy             }
        if delta1 < ABS(x1[i] - xalt[i]) then
          delta1 := ABS(x1[i] - xalt[i]);
      if delta1 <= 0.001 then
        if delta0 <= delta1 then
        begin
          krit := 3;                        { desired accuracy reached }
          goto 40
        end;
      delta0 := delta1;

      if m < maxit then
        for i:=0 to n-1 do xalt[i]:=x1[i]
    end
    else
      goto 40;

  end; {of m loop}


40: itanz := m;


  if sing then
    rc:= 2
  else if krit = 0 then
    rc:= 1
  else
    rc:=0;
50:End;

END.
{ -------------------------- END brown.pas -------------------------- }