-- simeq.adb  solve equations AX=Y

with real_arrays; use real_arrays;

function simeq ( A : real_matrix ; Y : real_vector ) return real_vector is
  --      purpose : solve the linear system of equations with real
  --                coefficients   [a] * |x| = |y|
  --
  --      input  : the real matrix  a
  --               the real vector  y
  --
  --      output : the real vector  x
  --
  --      method : GAUSS-jordan elimination using maximum element
  --               for pivot.
  --
  --      exception : matrix_data_error if 'a' is singular
  --                  array_index_error if 'a' is not square or
  --                                       length of 'y' /= rows of 'a'
  --
  --      usage  :     x := linear_equations ( a , y ) ;
  --
  --      written by : jon squire , 28 may 1983 (yes, its that old!)
  --      ORIGONAL DEC 1959 for IBM 650, TRANSLATED TO OTHER LANGUAGES   */
  --      e.g. FORTRAN converted to Ada converted to C                   */


  n : integer := a'length(1) ;          -- number of equations
  x : real_vector ( 1 .. n ) ;          -- result being computed
  b : real_matrix ( 1 .. n , 1 .. n + 1 ) ;  -- working matrix
  row : array ( 1 .. n ) of integer ;   -- row interchange indicies
  hold , i_pivot : integer ;            -- pivot indicies
  pivot : real ;                        -- pivot element value
  abs_pivot : real ;                    -- abs of pivot element
  norm1 : real := 0.0 ;                 -- 1 norm of matrix
  matrix_data_error : exception ;
begin
  if a'length ( 1 ) /= a'length ( 2 ) or a'length ( 1 ) /= y'length then
    raise array_index_error ;
  end if ;
  --                               build working data structure
  for i in 1 .. n loop
    for j in 1 .. n loop
      b ( i , j ) := a ( i - 1 + a'first( 1 ) , j - 1 + a'first ( 2 )) ;
      if abs b ( i , j ) > norm1 then
        norm1 := abs b ( i , j ) ;
      end if;
    end loop ;
    b ( i , n + 1 ) := y ( i - 1 + y'first ) ;
  end loop ;
  --                               set up row  interchange vectors
  for k in 1 .. n loop
    row ( k ) := k ;
  end loop ;
  --                               begin main reduction loop
  for k in 1 .. n loop
    --                             find largest element for pivot
    pivot := b ( row( k ) , k) ;
    abs_pivot := abs ( pivot ) ;
    i_pivot := k ;
    for i in k .. n loop
      if abs ( b( row( i ) , k)) > abs_pivot then
        i_pivot := i ;
        pivot := b ( row( i ) , k) ;
        abs_pivot := abs ( pivot ) ;
      end if ;
    end loop ;
    --                             check for near singular
    if abs_pivot < real'epsilon * norm1 then
      raise matrix_data_error;
    end if;
    --                             have pivot, interchange row pointers
    hold := row ( k ) ;
    row ( k ) := row ( i_pivot ) ;
    row ( i_pivot ) := hold ;
    --                           reduce about pivot
    for j in k + 1 .. n + 1 loop
      b ( row( k ) , j) := b ( row( k ) , j) / b ( row( k ) , k) ;
    end loop ;
    --                           inner reduction loop
    for i in 1 .. n loop
      if i /= k then
        for j in k + 1 .. n + 1 loop
          b ( row( i ) , j) := b ( row( i ) , j) - b ( row( i ) , k) * b (
             row( k ) , j) ;
        end loop ;
      end if ;
    end loop ;
    --                             finished inner reduction
  end loop ;
  --                               build  x  for return, unscrambling rows
  for i in 1 .. n loop
    x ( i ) := b ( row( i ) , n + 1) ;
  end loop ;
  return x ;
end simeq ;