-- test_simeq_newton4.adb  solve nonlinear system of equations
--                         method: newton iteration using Jacobian
--
--  Given problem  A X = Y  where X may have terms x1, x2, x3, x4, ...
--                 and higher order such as: x1*x2, x1*x1*x3, x4*x4*x4, ...
--                 A sparse matrix may be used, coefficients given
--                 Y is vector of real, given
--                 independent unknowns are x1, x2, x3, x4, ...
--
--  Solve by initial guess at values of x1, x2, x3, x4 computing products
--    X_next = X_initial - J_initial^-1 * (A *  X_initial - Y)
--
--  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    x4   x1*x2 x1*x1*x3  x4*x4*x4
--
-- | 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 |*| x4       |=| y4 |
-- | A5,1  A5,2  A5,3  A5,4  A5,5   A5,6    A5,7 | | x1*x2    | | y5 |
-- | A6,1  A6,2  A6,3  A6,4  A6,5   A6,6    A6,7 | | x1*x1*x3 | | y6 |
-- | A7,1  A7,2  A7,3  A7,4  A7,5   A7,6    A7,7 | | x4*x4*x4 | | 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,5*x2 + 2*A1,6*x3           deriv Row 1 wrt x1
--  J1,2 = A1,2 + A1,5*x1                       deriv Row 1 wrt x2
--  J1,3 = A1,3 + A1,6*x1*x1                    deriv Row 1 wrt x3
--  J1,4 = A1,4 + 3*A1,7*x4*x4                  deriv Row 1 wrt x4
--  J1,5 = A1,5                                 deriv Row 1 wrt x1, wrt x2
--  J1,6 = A1,6                                 deriv Row 1 wrt x1, wrt x3
--  J1,7 = A1,7                                 deriv Row 1 wrt x1, wrt x4
--
--  J2,1 = A2,2 + A2,5*x2 + 2*A2,6*x3           deriv Row 1 wrt x1
--  J2,2 = A2,2 + A2,5*x1                       deriv Row 1 wrt x2
--  J2,3 = A2,3 + A2,6*x1*x1                    deriv Row 1 wrt x3
--  J2,4 = A2,4 + 3*A2,7*x4*x4                  deriv Row 1 wrt x4
--  J2,5 = A2,5                                 deriv Row 1 wrt x1, wrt x2
--  J2,6 = A2,6                                 deriv Row 1 wrt x1, wrt x3
--  J2,7 = A2,7                                 deriv Row 1 wrt x1, wrt x4
--
--  etc.
--
--  J1,4 ... J4,7  should last ones be needed? may not be independent
--

with Ada.text_io;
use  Ada.text_io;
with real_arrays;
use  real_arrays;
with integer_arrays;
use  integer_arrays;
with simeq_newton4;

procedure test_simeq_newton4 is
  state : real := 30000000.0; -- arbitrary initial value

  -- 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_Terms(n:integer; nlin:integer; var1:integer_vector;
                        var2:integer_vector; var3:integer_vector;
                        var4:integer_vector) is
    names : array(1..6) of string(1..2) := ("x1","x2","x3","x4","x5","x6");
  begin
    put_Line("terms:");
    for i in 1..n loop
      put_Line(names(var1(i)));
    end loop;
    put_Line("power terms:");
    for i in n+1..n+nlin loop
      put(names(var1(i))&" * "&names(var2(i))); -- must be at least two
      if(var3(i)/=0) then
         put(" * "&names(var3(i)));
      end if;
      if(var4(i)/=0) then
         put(" * "&names(var4(i)));
      end if;
      new_Line;
    end loop;
  end print_terms;

  procedure test_case0 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);
    var4 : 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_0");

    -- build test case
    -- try "C" test case
    A(1,1):=0.041631;
    A(1,2):=0.176643;
    A(1,3):=0.364602;
    A(1,4):=0.091331;
    A(2,1):=0.092298;
    A(2,2):=0.487217;
    A(2,3):=0.526750;
    A(2,4):=0.454433;
    -- 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, 1, 2);
    var4 := (0, 0, 1, 2);
    -- choose x1=1.1 x2=2.3 , compute non-linear terms
    X_soln(1) := 1.1;
    X_soln(2) := 2.3;
    -- compute dependent terms
    for i in n+1..n+nlin loop
      X_soln(i) := X_soln(var1(i));  -- var1 must not be empty
      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(var4(i) /= 0) then
         X_soln(i) := X_soln(i)* X_soln(var4(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;

    -- initial condition
    X(1) := 1.0;
    X(2) := 1.0;
    -- debug print
    for i in 1..n loop
      put("X(");
      put(i,2);
      put(")=");
      put(X(i),3,6);
      new_line;
    end loop;
    new_line;
    Print_Terms(N, Nlin, var1, var2, var3, var4);

    simeq_newton4(n, nlin, A, Y, var1, var2, var3, var4, X, 12, 1.0e-6, 3);

    -- check solution
    put_line("test case0 solution");
    err := 0.0;
    for i in 1..n loop
      put("X(");
      put(i,2);
      put(")= ");
      put(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 residuals= ");
    put(err);
    new_line;
    put_Line("end test case0");
    new_line;
  end test_case0;

  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);
    var4 : 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) := 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, 1, 2);
    var2 := (0, 0, 1, 2);
    var3 := (0, 0, 1, 2);
    var4 := (0, 0, 0, 2);
    -- choose x1=1.1 x2=2.3 , compute non-linear terms
    X_soln(1) := 1.1;
    X_soln(2) := 2.3;
    -- compute dependent terms
    for i in n+1..n+nlin loop
      X_soln(i) := X_soln(var1(i));  -- var1 must not be empty
      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(var4(i) /= 0) then
         X_soln(i) := X_soln(i)* X_soln(var4(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;

    -- initial condition
    X(1) := 1.0;
    X(2) := 1.0;
    -- debug print
    for i in 1..n loop
      put("X(");
      put(i,2);
      put(")=");
      put(X(i),3,6);
      new_line;
    end loop;
    new_line;
    Print_Terms(N, Nlin, var1, var2, var3, var4);

    simeq_newton4(n, nlin, A, Y, var1, var2, var3, var4, X, 12, 1.0e-6, 3);

    -- 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));
      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 residuals= ");
    put(err);
    new_line;
    put_Line("end test case1");
    new_line;
  end test_case1;

  procedure test_case2 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);
    var4 : 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) := 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);
    var4 := (0, 0, 0, 0, 0, 3, 4);

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

    -- compute dependent terms
    for i in n+1..n+nlin loop
      X_soln(i) := X_soln(var1(i));  -- var1 must not be empty
      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(var4(i) /= 0) then
         X_soln(i) := X_soln(i)* X_soln(var4(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;

    -- initial condition
    X(1) := 1.0;
    X(2) := 1.0;
    X(3) := 1.0;
    X(4) := 1.0;
    X(5) := X(1)*X(2);
    X(6) := X(1)*X(1)*X(3);
    X(7) := X(4)*X(4)*X(4);
    -- debug print
    for i in 1..n loop
      put("X(");
      put(i,2);
      put(")=");
      put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    Print_Terms(N, Nlin, var1, var2, var3, var4);

    simeq_newton4(n, nlin, A, Y, var1, var2, var3, var4, X, 12, 1.0e-6, 3);

    -- 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));
      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 residuals= ");
    put(err);
    new_line;
    put_Line("end test case2");
    new_line;
  end test_case2;

  procedure test_case3 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);
    var4 : 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, 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);
    var4 := (0, 0, 0, 0, 0, 0, 3, 4, 3, 4);

    -- 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) := X_soln(var1(i));  -- var1 must not be empty
      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(var4(i) /= 0) then
         X_soln(i) := X_soln(i)* X_soln(var4(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;

    -- initial condition
    X(1) := 0.7;
    X(2) := 0.75;
    X(3) := 0.8;
    X(4) := 0.85;
    X(5) := X(1)*X(1)*X(1);
    X(6) := X(2)*X(2)*X(2);
    X(7) := X(3)*X(3)*X(3);
    X(8) := X(4)*X(4)*X(4);
    X(9) := X(1)*X(2)*X(3);
    X(10):= X(2)*X(3)*X(4);
    -- debug print
    for i in 1..n loop
      put("X_init(");
      put(i,2);
      put(")=");
      put(X(i),3,6);
      new_line;
    end loop;
    new_line;

    Print_Terms(N, Nlin, var1, var2, var3, var4);

    simeq_newton4(n, nlin, A, Y, var1, var2, var3, var4, X, 17, 1.0e-5, 5);

    -- 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));
      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 residuals= ");
    put(err);
    new_line;
    put_Line("end test case3");
    new_line;
  end test_case3;

  procedure Test_Case_Big is
    n : integer := 4;    -- number of linear unknowns
    nlin : integer := 65; -- number of non linear terms
    var1 : integer_vector(1..n+nlin);
    var2 : integer_vector(1..n+nlin);
    var3 : integer_vector(1..n+nlin);
    var4 : 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);
    Xchk : real_vector(1..n);
    err, aerr : real;           -- error from test solution
  begin
    put_line("test_case_big");

    -- build test case
    for i in 1..n loop
      for j in 1..n+nlin loop
        A(i,j) := udrnrt;
      end loop;
    end loop;

    -- setup definition of linear and nonlinear terms
    var1 := (1, 2, 3, 4, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4,
             1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 4,
             1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
             2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4);
    var2 := (0, 0, 0, 0, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4,
             1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 2, 2, 2, 3, 3, 4, 3, 3, 4, 4,
             1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 4,
             2, 2, 2, 2, 2, 2, 3, 3, 3, 4, 3, 3, 3, 4, 4);
    var3 := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
             1, 2, 3, 4, 2, 3, 4, 3, 4, 4, 2, 3, 4, 3, 4, 4, 3, 4, 4, 4,
             1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 2, 2, 2, 3, 3, 4, 3, 3, 4, 4,
             2, 2, 2, 3, 3, 4, 3, 3, 4, 4, 3, 3, 4, 4, 4);
    var4 := (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
             1, 2, 3, 4, 2, 3, 4, 3, 4, 4, 2, 3, 4, 3, 4, 4, 3, 4, 4, 4,
             2, 3, 4, 3, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4, 4);

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

    -- compute dependent terms
    for i in n+1..n+nlin loop
      X_soln(i) := X_soln(var1(i));  -- var1 must not be empty
      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(var4(i) /= 0) then
         X_soln(i) := X_soln(i)* X_soln(var4(i));
      end if;
    end loop;

    -- 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;

    -- initial condition
    for i in 1..n+nlin loop
      X(i) := X_soln(i);
    end loop;

    Print_Terms(N, Nlin, var1, var2, var3, var4);


    simeq_newton4(n, nlin, A, Y, var1, var2, var3, var4, X,  10, 1.0e-6, 1);

    -- check solution should be solved
    put_line("test case_big solution  init = solution");
    err := 0.0;
    for i in 1..n loop
      put("X(");
      put(i,2);
      put(")= ");
      put(X(i));
      Put(", err=");
      Put(X(i)-Xchk(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 residuals= ");
    put(err);
    new_line;

    for i in 1..n+nlin loop
      X(i) := 1.0;
    end loop;

    simeq_newton4(n, nlin, A, Y, var1, var2, var3, var4, X,  15, 1.0e-6, 1);

    -- check solution
    put_line("test case_big solution, init 1.0");
    err := 0.0;
    for i in 1..n loop
      put("X(");
      put(i,2);
      put(")= ");
      put(X(i));
      Put(", err=");
      Put(X(i)-Xchk(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 residuals= ");
    put(err);
    new_line;

    put_Line("end test case_big");
    new_line;
  end Test_Case_Big;

begin
  put_line("test_simeq_newton4.adb running");
  test_case0;
  test_case1;
  test_case2;
  test_case3; -- fails
  Test_Case_Big;

  put_line("end test_simeq_newton4.adb");
end test_simeq_newton4;