let basic_solve a b =
let m = Array.length a in
let n = Array.length a.(0) in
let j i =
try
for j = 0 to n-1 do
if not (nearly_zero a.(i).(j)) then raise (Found j)
done ; -1
with Found j -> j
in
let x = Array.make n 0. in
for i = m-1 downto 0 do
let j = j i in
if j = -1 then
( if abs_float b.(i) > 0.1 then
( Printf.printf "Error at %d: 0. <> %f\n%!" i b.(i) ;
raise Unsolvable ) else
Printf.printf "Admitted error %f at %d.\n%!" b.(i) i )
else
( x.(j) <- b.(i) /. a.(i).(j) ;
for ii = 0 to i do
b.(ii) <- b.(ii) -. x.(j) *. a.(ii).(j)
done )
done ;
x