type 'a matrix = 'a array array;;

let make rows cols = Array.make_matrix rows cols 0.;;
let init rows cols f = 
  Array.init rows (fun row -> Array.init cols (fun col -> f row col));
;;

let rows = Array.length;;
let cols mat = Array.length mat.(0);;

let float_sum mat = 
  let sum      = ref 0. in 
  let add i    = sum := !sum +. abs_float i in 
  let iter row = add (Array.fold_left ( +. ) 0. row) in 
  let _        = Array.iter iter mat in !sum;
;;

let iter f matrix = Array.iter (Array.iter f) matrix;;
let map f matrix  = Array.map (Array.map f) matrix;;
let dmap f matrix =
  for i = 0 to rows matrix - 1 do
    let row = matrix.(i) in
      for j = 0 to cols matrix - 1 do
        row.(j) <- f row.(j);
      done;
  done;
  matrix;
;;

let add a b =
  assert (rows a = rows b);
  assert (cols a = cols b);
  let c = init (rows a) (cols a) (fun i j -> a.(i).(j) +. b.(i).(j)) in c;
;;

let scalar_div a s = dmap (fun v -> v /. s) a;;

let print matrix file  =
  for i = 0 to rows matrix - 1 do
    Array.iter (Printf.fprintf file "%f ") matrix.(i);
    Printf.fprintf file "\n";
  done;
;;