-- test_simeq_newton5.adb  solve nonlinear system of equations
--                         method: newton iteration using Jacobian
--                    use list for higher order terms
--
--  Given problem  A X = Y  where X may have terms x1, x2, x3, x4, ...
--                 and nonlinear such as: x1*x2, x1^2*x3, x3/x2, x3^3/x2^2 ...
--                 A sparse matrix may be used, coefficients given
--                 Y is vector of real, given
--                 independent unknowns are x1, ... , xn
--        a maximal term can have up to power of 3 divided by power of 2
--        in this implementation (if we had symbolic differentiation, no limit)
--        a term could be  x1 , x1*x3 , x2/x5 , x1*x2*x3/(x4*x5) ,  x2^3/x4^2
--
--  Solve by initial guess at values of x1, x2, x3, x4 computing products
--    X_next = X_initial - J_initial^-1 * (A *  X_initial - Y)
--    in general X_next = X_prev - (J_prev^-1 * (A * X_prev - Y))*b
--                                 where 0 < b < 1, often 0.5, for stability
--
--  solved when  abs sum each row A * X_next -Y < epsilon
--
--  There may be no solution
--  It may stall, stop if abs(X_next-X_prev)<epsilon
--  It may diverge, stop, indicate no solution
--                        (or try a different initial guess)
--  It may oscillate, stop, indicate no solution
--                        (or try a different initial guess)
--
--  The matrix equation could be:
--
--   x1    x2    x3   x1*x2 x3*x2^2 x2/x3 x1^3/x3^2 ...
--
-- | A1,1  A1,2  A1,3  A1,4   A1,5  A1,6  A1,7 ...| | x1         | |y1 |
-- | A2,1  A2,2  A2,3  A2,4   A2,5  A2,6  A2,7 ...| | x2         | |y2 |
-- | A3,1  A3,2  A3,3  A3,4   A3,5  A3,6  A3,7 ...| | x3         | |y3 |
-- | A4,1  A4,2  A4,3  A4,4   A4,5  A4,6  A4,7 ...| | x1*x2      | |y4 |
-- | A5,1  A5,2  A5,3  A5,4   A5,5  A5,6  A5,7 ...|*| x3*x2^2    |=|y5 |
-- | A6,1  A6,2  A6,3  A6,4   A6,5  A6,6  A6,7 ...| | x2/x3      | |y6 |
-- | A7,1  A7,2  A7,3  A7,4   A7,5  A7,6  A7,7 ...| | x1^3/x3^2  | |y7 |
-- | ...                                       ...| |            | |   |
--
--  The Jacobian, J is the numerically computed derivatives based on A, X_prev
--  The derivative coefficients may be computed once based on A and the
--  unknown Variables. The numeric value plugs in the X_prev values.
--
--  J1,1 = A1,1 + A1,4*x2 + 3*A1,7*x1^2/x3^2               deriv Row 1 wrt x1
--  J1,2 = A1,2 + A1,4*x1 + 2*A1,5*x3*x2 + A1,6/x3         deriv Row 1 wrt x2
--  J1,3 = A1,3 + A1,5*x2^2 + -A1,6*x2/x3^2 + -2*A1,7*x1^2/x3^3
--                                                         deriv Row 1 wrt x3
--  J2,1 = A2,1 + A2,4*x2 + 3*A2,7*x1^2/x3^2               deriv Row 2 wrt x1
--  J2,2 = A2,2 + A2,4*x1 + 2*A2,5*x3*x2 + A2,6/x3         deriv Row 2 wrt x2
--  J2,3 = A1,3 + A2,5*x2^2 + -A2,6*x2/x3^2 + -2*A2,7*x1^2/x3^3
--                                                         deriv Row 2 wrt x3
--  J3,1 = A3,1 + A3,4*x2 + 3*A3,7*x1^2/x3^2               deriv Row 3 wrt x1
--  J3,2 = A3,2 + A3,4*x1 + 2*A3,5*x3*x2 + A3,6/x3         deriv Row 3 wrt x2
--  J3,3 = A3,3 + A3,5*x2^2 + -A3,6*x2/x3^2 + -2*A3,7*x1^2/x3^3
--                                                         deriv Row 3 wrt x3
--  ...
--  etc.
--
--  Code or a data structure must be available to know the equation
--  of the entries in the Y vector. Symbolic computation of
--  derivatives is assumed to be available, when needed.

with Ada.text_io;
use  Ada.text_io;
with Ada.Numerics.Float_Random;
use  Ada.Numerics.Float_Random;
with real_arrays;
use  real_arrays;
with integer_arrays;
use  integer_arrays;
with simeq_newton5;

procedure test_simeq_newton5 is
  state : real := 30000000.0; -- arbitrary initial value for udrnrt
  S : Generator; -- for Ada.Numerics.Float_Random

  -- udrnrt.adb uniformly distributed random numbers 0.0 .. 1.0
  function udrnrt return real is
    modulus : real := 2147483647.0; -- 2**31-1
    multiplier : real := 16807.0; -- prime to prime power
    t : long_float;
    i : integer;
  begin
    t := state * multiplier;
    i := integer(t / Modulus);
    state := t - (long_float(i) * modulus);
    if state < 0.0 then
     state := state + modulus;
    end if;
    return state / modulus;
  end udrnrt;

  package fltio is new float_io(real);
  use fltio;
  package intio is new integer_io(integer);
  use intio;

  procedure Print_Equation(n : integer; nlin : integer;
                           Var1 : integer_vector; Var2 : integer_vector;
                           Var3 : integer_vector; Vari1 : integer_vector;
                                                  Vari2 : Integer_Vector) is
    Xnames : array(1..5) of String(1..2) := ("X1", "X2", "X3", "X4", "X5");
  begin
    Put_Line("system of equations to be solved, i=1.."&Integer'Image(n));
    for j in 1..n loop
      Put_line("A(i,"&Integer'Image(j)&")*"&Xnames(j)&"+");
    end loop;
    for j in n+1..n+nlin loop
      Put("A(i,"&Integer'Image(J)&")");
      if Var1(J)>0 then
         Put("*"&Xnames(Var1(J)));
      end if;
      if Var2(J)>0 then
         Put("*"&Xnames(Var2(J)));
      end if;
      if Var3(J)>0 then
         Put("*"&Xnames(Var3(J)));
      end if;
      if Vari1(J)>0  then
         Put("/("&Xnames(Vari1(J)));
      end if;
      if Vari2(J)>0  then
         Put("*"&Xnames(Vari2(J)));
      end if;
      if Vari1(J)>0 or Vari2(J)>0 then
         Put(")");
      end if;
      if J /= n+nlin then
        Put_Line("+");
      else
        New_Line;
      end if;
    end loop;
    Put_Line(" = Y(i)");
    New_Line;
  end Print_Equation;

  procedure test_case1 is
    n : integer := 2;    -- number of linear unknowns
    nlin : integer := 2; -- number of non linear terms
    Var1 : integer_vector(1..n+nlin);
    Var2 : integer_vector(1..n+nlin);
    Var3 : integer_vector(1..n+nlin);
    Vari1 : integer_vector(1..n+nlin);
    Vari2 : integer_vector(1..n+nlin);
    A : real_matrix(1..n,1..n+nlin);
    Y : real_vector(1..n);
    X : real_vector(1..n+nlin);
    X_soln : real_vector(1..n+nlin);
    err, aerr : real;           -- error from test solution
  begin
    Put_line("test_case_1");

    -- build test case
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := 1.0;
      end loop;
    end loop;
    A(1,4) := 0.0; -- x + y + x^2 = 4
    A(2,3) := 0.0; -- x + y + y^2 = 7
    -- debug print
    for i in 1..n loop
      for j in 1..n+nlin loop
        Put("A(");
        Put(i,2);
        Put(",");
        Put(j,2);
        Put(")=");
        Put(A(i,j),3,6);
        new_line;
      end loop;
    end loop;
    new_line;

    -- setup definition of linear and nonlinear terms
    Var1  := (1, 2, 1, 2);
    Var2  := (0, 0, 1, 2);
    Var3  := (0, 0, 0, 0);
    Vari1 := (0, 0, 0, 0);
    Vari2 := (0, 0, 0, 0);

    -- choose x1=1 x2=2 , compute x1*x2, etc.
    X_soln(1) := 1.0;
    X_soln(2) := 2.0;
    -- compute dependent terms
    for I in n+1..n+nlin loop
      X_Soln(I) := 1.0;
      if Var1(I) > 0 then
        X_soln(i) := X_soln(Var1(i));
      end if;
      if Var2(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var2(i));
      end if;
      if Var3(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var3(i));
      end if;
      if Vari1(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari1(i));
      end if;
      if Vari2(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari2(i));
      end if;
    end loop;

    -- debug print
    for i in 1..n loop
      Put("X_soln(");
      Put(i,2);
      Put(")=");
      Put(X_soln(i),3,6);
      new_line;
    end loop;
    new_line;

    -- set up Y
    for i in 1..n loop
      Y(i) := 0.0;
      for j in 1..n+nlin loop
        Y(i) := Y(i) + A(i,j)*X_soln(j);
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      Put("Y(");
      Put(i,2);
      Put(")=");
      Put(Y(i),3,6);
      new_line;
    end loop;
    new_line;
    Print_Equation(n, nlin, Var1, Var2, Var3, Vari1, Vari2);

    -- initial condition
    X(1) := 0.5;
    X(2) := 0.5;
    -- higher terms not needed

    -- debug print
    Put_Line("initial guess");
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")=");
      Put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X);

    -- check solution
    Put_line("test case1 solution");
    err := 0.0;
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")= ");
      Put(X(i));
      Put(", soln=");
      Put(X_Soln(i));
      Put(", err=");
      Put(X_Soln(i)-X(i));
      new_line;

      aerr := 0.0;
      for j in 1..n+nlin loop
        aerr := aerr + A(i,j)*X(j);
      end loop;
      aerr := aerr-Y(i);
      err := err + abs(aerr);
    end loop;
    new_line;
    Put("                        residual= ");
    Put(err);
    new_line;
    new_line;
  end test_case1;

  procedure test_case2 is
    n : integer := 3;    -- number of linear unknowns
    nlin : integer := 4; -- number of non linear terms
    Var1 : integer_vector(1..n+nlin);
    Var2 : integer_vector(1..n+nlin);
    Var3 : integer_vector(1..n+nlin);
    Vari1 : integer_vector(1..n+nlin);
    Vari2 : integer_vector(1..n+nlin);
    A : real_matrix(1..n,1..n+nlin);
    Y : real_vector(1..n);
    X : real_vector(1..n+nlin);
    X_soln : real_vector(1..n+nlin);
    err, aerr : real;           -- error from test solution
  begin
    Put_line("test_case_2");

    -- build test case
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := Real(Random(S));
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      for j in 1..n+nlin loop
        Put("A(");
        Put(i,2);
        Put(",");
        Put(j,2);
        Put(")=");
        Put(A(i,j),3,6);
        new_line;
      end loop;
    end loop;
    new_line;

    -- setup definition of linear and nonlinear terms
    Var1  := (1, 2, 3, 1, 3, 2, 1);
    Var2  := (0, 0, 0, 2, 2, 0, 1);
    Var3  := (0, 0, 0, 0, 2, 0, 1);
    Vari1 := (0, 0, 0, 0, 0, 3, 3);
    Vari2 := (0, 0, 0, 0, 0, 0, 3);
    -- choose x1=1.1 x2=1.2 , compute x1*x2, etc.
    X_soln(1) := 1.1;
    X_soln(2) := 1.2;
    X_Soln(3) := 1.4;
    -- compute dependent terms
    for I in n+1..n+nlin loop
      X_Soln(I) := 1.0;
      if Var1(I) > 0 then
        X_soln(i) := X_soln(Var1(i));
      end if;
      if Var2(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var2(i));
      end if;
      if Var3(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var3(i));
      end if;
      if Vari1(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari1(i));
      end if;
      if Vari2(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari2(i));
      end if;
    end loop;

    -- debug print
    for i in 1..n loop
      Put("X_soln(");
      Put(i,2);
      Put(")=");
      Put(X_soln(i),3,6);
      new_line;
    end loop;
    new_line;

    -- set up Y
    for i in 1..n loop
      Y(i) := 0.0;
      for j in 1..n+nlin loop
        Y(i) := Y(i) + A(i,j)*X_soln(j);
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      Put("Y(");
      Put(i,2);
      Put(")=");
      Put(Y(i),3,6);
      new_line;
    end loop;
    new_line;
    Print_Equation(n, nlin, Var1, Var2, Var3, Vari1, Vari2);

    -- initial condition
    X(1) := 1.0;
    X(2) := 1.0;
    X(3) := 1.0;
    -- higher terms not needed

    -- debug print
    Put_Line("initial guess");
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")=");
      Put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X);

    -- check solution
    Put_line("test case2 solution");
    err := 0.0;
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")= ");
      Put(X(i));
      Put(", soln=");
      Put(X_Soln(i));
      Put(", err=");
      Put(X_Soln(i)-X(i));
      new_line;

      aerr := 0.0;
      for j in 1..n+nlin loop
        aerr := aerr + A(i,j)*X(j);
      end loop;
      aerr := aerr-Y(i);
      err := err + abs(aerr);
    end loop;
    new_line;
    Put("                        residual= ");
    Put(err);
    new_line;
    new_line;
  end test_case2;

  procedure test_case3 is
    n : integer := 4;    -- number of linear unknowns
    nlin : integer := 3; -- number of non linear terms
    Var1 : integer_vector(1..n+nlin);
    Var2 : integer_vector(1..n+nlin);
    Var3 : integer_vector(1..n+nlin);
    Vari1 : integer_vector(1..n+nlin);
    Vari2 : integer_vector(1..n+nlin);
    A : real_matrix(1..n,1..n+nlin);
    Y : real_vector(1..n);
    X : real_vector(1..n+nlin);
    X_soln : real_vector(1..n+nlin);
    err, aerr : real;           -- error from test solution
  begin
    Put_line("test_case_3");

    -- build test case
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := udrnrt;
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      for j in 1..n+nlin loop
        Put("A(");
        Put(i,2);
        Put(",");
        Put(j,2);
        Put(")=");
        Put(A(i,j),3,6);
        new_line;
      end loop;
    end loop;
    new_line;

    -- setup definition of linear and nonlinear terms
    Var1  := (1, 2, 3, 4, 1, 1, 4);
    Var2  := (0, 0, 0, 0, 2, 1, 4);
    Var3  := (0, 0, 0, 0, 0, 3, 4);
    Vari1 := (0, 0, 0, 0, 0, 0, 0);
    Vari2 := (0, 0, 0, 0, 0, 0, 0);

    -- choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, etc.
    X_soln(1) := 1.1;
    X_soln(2) := 1.2;
    X_soln(3) := 1.4;
    X_soln(4) := 1.5;
    -- X_soln(5..7) :=  -- computed from "Var's"

    -- compute dependent terms
    for I in n+1..n+nlin loop
      X_Soln(I) := 1.0;
      if Var1(I) > 0 then
        X_soln(i) := X_soln(Var1(i));
      end if;
      if Var2(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var2(i));
      end if;
      if Var3(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var3(i));
      end if;
      if Vari1(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari1(i));
      end if;
      if Vari2(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari2(i));
      end if;
    end loop;

    -- debug print
    for i in 1..n loop
      Put("X_soln(");
      Put(i,2);
      Put(")=");
      Put(X_soln(i),3,6);
      new_line;
    end loop;
    new_line;

    -- set up Y
    for i in 1..n loop
      Y(i) := 0.0;
      for j in 1..n+nlin loop
        Y(i) := Y(i) + A(i,j)*X_soln(j);
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      Put("Y(");
      Put(i,2);
      Put(")=");
      Put(Y(i),3,6);
      new_line;
    end loop;
    new_line;
    Print_Equation(n, nlin, Var1, Var2, Var3, Vari1, Vari2);

    -- initial condition
    X(1) := 1.0;
    X(2) := 1.0;
    X(3) := 1.0;
    X(4) := 1.0;

    -- debug print
    Put_Line("initial guess");
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")=");
      Put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X);

    -- check solution
    Put_line("test case3 solution");
    err := 0.0;
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")= ");
      Put(X(i));
      Put(", soln=");
      Put(X_Soln(i));
      Put(", err=");
      Put(X_Soln(i)-X(i));
      new_line;

      aerr := 0.0;
      for j in 1..n+nlin loop
        aerr := aerr + A(i,j)*X(j);
      end loop;
      aerr := aerr-Y(I);
      err := err + abs(aerr);
    end loop;
    new_line;
    Put("                        residual= ");
    Put(err);
    new_line;
    new_line;
  end test_case3;

  procedure test_case4 is
    n : integer := 4;    -- number of linear unknowns
    nlin : integer := 6; -- number of non linear terms
    Var1 : integer_vector(1..n+nlin);
    Var2 : integer_vector(1..n+nlin);
    Var3 : integer_vector(1..n+nlin);
    Vari1 : integer_vector(1..n+nlin);
    Vari2 : integer_vector(1..n+nlin);
    A : real_matrix(1..n,1..n+nlin);
    Y : real_vector(1..n);
    X : real_vector(1..n+nlin);
    X_soln : real_vector(1..n+nlin);
    err, aerr : real;           -- error from test solution
  begin
    Put_line("test_case_4");

    -- build test case
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := udrnrt;
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      for j in 1..n+nlin loop
        Put("A(");
        Put(i,2);
        Put(",");
        Put(j,2);
        Put(")=");
        Put(A(i,j),3,6);
        new_line;
      end loop;
    end loop;
    new_line;

    -- setup definition of linear and nonlinear terms
    Var1  := (1, 2, 3, 4, 1, 2, 3, 4, 1, 2);
    Var2  := (0, 0, 0, 0, 1, 2, 3, 4, 2, 3);
    Var3  := (0, 0, 0, 0, 1, 2, 3, 4, 3, 4);
    Vari1 := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
    Vari2 := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0);

    -- choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, etc.
    X_soln(1) := 1.1;
    X_soln(2) := 1.2;
    X_soln(3) := 1.4;
    X_soln(4) := 1.5;
    -- X_soln(5..10) :=  -- computed from "Var's"

    -- compute dependent terms
    for I in n+1..n+nlin loop
      X_Soln(I) := 1.0;
      if Var1(I) > 0 then
        X_soln(i) := X_soln(Var1(i));
      end if;
      if Var2(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var2(i));
      end if;
      if Var3(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var3(i));
      end if;
      if Vari1(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari1(i));
      end if;
      if Vari2(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari2(i));
      end if;
    end loop;

    -- debug print
    for i in 1..n loop
      Put("X_soln(");
      Put(i,2);
      Put(")=");
      Put(X_soln(i),3,6);
      new_line;
    end loop;
    new_line;

    -- set up Y
    for i in 1..n loop
      Y(i) := 0.0;
      for j in 1..n+nlin loop
        Y(i) := Y(i) + A(i,j)*X_soln(j);
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      Put("Y(");
      Put(i,2);
      Put(")=");
      Put(Y(i),3,6);
      new_line;
    end loop;
    new_line;

    Print_Equation(n, nlin, Var1, Var2, Var3, Vari1, Vari2);

    -- initial condition
    X(1) := 0.7;
    X(2) := 0.75;
    X(3) := 0.8;
    X(4) := 0.85;

    -- debug print
    Put_Line("initial guess");
    for i in 1..n loop
      Put("X_init(");
      Put(i,2);
      Put(")=");
      Put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X,
                  12, 1.0e-5, 1);

    -- check solution
    Put_line("test case4 solution");
    err := 0.0;
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")= ");
      Put(X(i));
      Put(", soln=");
      Put(X_Soln(i));
      Put(", err=");
      Put(X_Soln(i)-X(i));
      new_line;

      aerr := 0.0;
      for j in 1..n+nlin loop
        aerr := aerr + A(i,j)*X(j);
      end loop;
      aerr := aerr-Y(I);
      err := err + abs(aerr);
    end loop;
    new_line;
    Put("                        sum of errors= ");
    Put(err);
    new_line;
    new_line;
  end test_case4;

  procedure test_case5 is
    n : integer := 5;    -- number of linear unknowns
    nlin : integer := 15; -- number of non linear terms
    Var1 : integer_vector(1..n+nlin);
    Var2 : integer_vector(1..n+nlin);
    Var3 : integer_vector(1..n+nlin);
    Vari1 : integer_vector(1..n+nlin);
    Vari2 : integer_vector(1..n+nlin);
    A : real_matrix(1..n,1..n+nlin);
    Y : real_vector(1..n);
    X : real_vector(1..n+nlin);
    X_soln : real_vector(1..n+nlin);
    err, aerr : real;           -- error from test solution
  begin
    Put_line("test_case_5");

    -- build test case, first and last equation solved
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := 0.0;
      end loop;
      A(i,i) := 1.0;
    end loop;
    for i in 2..n-1 loop
      for j in n+1..n+nlin loop
        A(i,j) := udrnrt;
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      for j in 1..n+nlin loop
        Put("A(");
        Put(i,2);
        Put(",");
        Put(j,2);
        Put(")=");
        Put(A(i,j),3,6);
        new_line;
      end loop;
    end loop;
    new_line;

    -- setup definition of linear and nonlinear terms
    Var1  := (1, 2, 3, 4, 5, 2, 2, 2, 3, 3, 4,
              2, 2, 2, 3, 3, 3, 4, 4, 4);
    Var2  := (0, 0, 0, 0, 0, 2, 3, 4, 3, 4, 4,
              2, 2, 2, 3, 3, 3, 4, 4, 4);
    Var3  := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              2, 3, 4, 2, 3, 4, 2, 3, 4);
    Vari1 := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0, 0);
    Vari2 := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 0, 0);

    -- choose x1=1.1 x2=1.2 x3=1.4 x4 =1.5, compute x1*x2, etc.
    X_soln(1) := 1.1;
    X_soln(2) := 1.2;
    X_soln(3) := 1.4;
    X_soln(4) := 1.5;
    X_soln(5) := 1.7;

    -- compute dependent terms
    for I in n+1..n+nlin loop
      X_Soln(I) := 1.0;
      if Var1(I) > 0 then
        X_soln(i) := X_soln(Var1(i));
      end if;
      if Var2(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var2(i));
      end if;
      if Var3(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var3(i));
      end if;
      if Vari1(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari1(i));
      end if;
      if Vari2(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari2(i));
      end if;
    end loop;

    -- debug print
    for i in 1..n loop
      Put("X_soln(");
      Put(i,2);
      Put(")=");
      Put(X_soln(i),3,6);
      new_line;
    end loop;
    new_line;

    -- set up Y
    for i in 1..n loop
      Y(i) := 0.0;
      for j in 1..n+nlin loop
        Y(i) := Y(i) + A(i,j)*X_soln(j);
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      Put("Y(");
      Put(i,2);
      Put(")=");
      Put(Y(i),3,6);
      new_line;
    end loop;
    new_line;

    Print_Equation(n, nlin, Var1, Var2, Var3, Vari1, Vari2);

    -- initial condition
    for i in 1..n loop
      X(i) := 1.0;
    end loop;
    -- debug print
    Put_Line("initial guess");
    for i in 1..n loop
      Put("X_init(");
      Put(i,2);
      Put(")=");
      Put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X,
                  24, 1.0e-4, 1);

    -- check solution
    Put_line("test case5 solution, computed, expected not unique");
    err := 0.0;
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")= ");
      Put(X(i));
      Put(", soln=");
      Put(X_Soln(i));
      Put(", err=");
      Put(X_Soln(i)-X(i));
      new_line;

      aerr := 0.0;
      for j in 1..n+nlin loop
        aerr := aerr + A(i,j)*X(j);
      end loop;
      aerr := aerr-Y(I);
      err := err + abs(aerr);
    end loop;
    new_line;
    Put("                        sum of errors= ");
    Put(err);
    new_line;
    new_line;

    Put_Line("check when given solution as initial condition");
    -- initial condition
    for i in 1..n loop
      X(i) := X_Soln(i);
    end loop;
    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X,
                  4, 1.0e-4, 1);
    new_line;
  end test_case5;

  procedure test_case6 is
    n : integer := 5;    -- number of linear unknowns
    nlin : integer := 1; -- number of non linear terms
    Var1 : integer_vector(1..n+nlin);
    Var2 : integer_vector(1..n+nlin);
    Var3 : integer_vector(1..n+nlin);
    Vari1 : integer_vector(1..n+nlin);
    Vari2 : integer_vector(1..n+nlin);
    A : real_matrix(1..n,1..n+nlin);
    Y : real_vector(1..n);
    X : real_vector(1..n+nlin);
    X_soln : real_vector(1..n+nlin);
    err, aerr : real;           -- error from test solution
  begin
    Put_line("test_case_6");

    -- build test case
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := udrnrt;
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      for j in 1..n+nlin loop
        Put("A(");
        Put(i,2);
        Put(",");
        Put(j,2);
        Put(")=");
        Put(A(i,j),3,6);
        new_line;
      end loop;
    end loop;
    new_line;

    -- setup definition of linear and nonlinear terms
    Var1  := (1, 2, 3, 4, 5, 1);
    Var2  := (0, 0, 0, 0, 0, 2);
    Var3  := (0, 0, 0, 0, 0, 3);
    Vari1 := (0, 0, 0, 0, 0, 4);
    Vari2 := (0, 0, 0, 0, 0, 5);

    -- choose x1=1.1 x2=1.2 x3=1.4 x4=1.5, x5=1.9
    X_soln(1) := 1.1;
    X_soln(2) := 1.2;
    X_soln(3) := 1.4;
    X_soln(4) := 1.5;
    X_soln(5) := 1.9;
    -- X_soln(6) :=  -- computed from "Var's"

    -- compute dependent terms
    for I in n+1..n+nlin loop
      X_Soln(I) := 1.0;
      if Var1(I) > 0 then
        X_soln(i) := X_soln(Var1(i));
      end if;
      if Var2(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var2(i));
      end if;
      if Var3(i) > 0  then
         X_soln(i) := X_soln(i) * X_soln(Var3(i));
      end if;
      if Vari1(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari1(i));
      end if;
      if Vari2(i) > 0  then
         X_soln(i) := X_soln(i) / X_soln(Vari2(i));
      end if;
    end loop;

    -- debug print
    for i in 1..n loop
      Put("X_soln(");
      Put(i,2);
      Put(")=");
      Put(X_soln(i),3,6);
      new_line;
    end loop;
    new_line;

    -- set up Y
    for i in 1..n loop
      Y(i) := 0.0;
      for j in 1..n+nlin loop
        Y(i) := Y(i) + A(i,j)*X_soln(j);
      end loop;
    end loop;
    -- debug print
    for i in 1..n loop
      Put("Y(");
      Put(i,2);
      Put(")=");
      Put(Y(i),3,6);
      new_line;
    end loop;
    new_line;
    Print_Equation(n, nlin, Var1, Var2, Var3, Vari1, Vari2);

    -- initial condition
    X(1) := 0.7;
    X(2) := 0.75;
    X(3) := 0.5;
    X(4) := 0.9;
    X(5) := 0.6;

    -- debug print
    Put_Line("initial guess");
    for i in 1..n loop
      Put("X_init(");
      Put(i,2);
      Put(")=");
      Put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    simeq_newton5(n, nlin, A, Y, Var1, Var2, Var3, Vari1, Vari2, X,
                  24, 1.0e-3, 1);

    -- check solution
    Put_line("test case6 solution");
    err := 0.0;
    for i in 1..n loop
      Put("X(");
      Put(i,2);
      Put(")= ");
      Put(X(i));
      Put(", soln=");
      Put(X_Soln(i));
      Put(", err=");
      Put(X_Soln(i)-X(i));
      new_line;

      aerr := 0.0;
      for j in 1..n+nlin loop
        aerr := aerr + A(i,j)*X(j);
      end loop;
      aerr := aerr-Y(I);
      err := err + abs(aerr);
    end loop;
    new_line;
    Put("                        sum of errors= ");
    Put(err);
    new_line;
    new_line;
  end test_case6;

begin
  Put_line("test_simeq_newton5.adb running");
  Put_Line("a random number="&Float'Image(Random(S)));

  test_case1;
  test_case2;
  test_case3;
  test_case4; -- close but fails
  test_case5; -- tough
  test_case6; -- stress max nonlinear

  Put_line("end test_simeq_newton5");
end test_simeq_newton5;