struct
  type ('a, 'b, 'c) matrix =
      { mat_i : 'a array ;
        mat_j : 'b array ;
        mat_mat : 'c array array ;
      }

  let latex oc ?(cols=10)
      string_of_i string_of_j string_of_v m =
    let p = Printf.fprintf in
    let i_elements = m.mat_i in
    let j_elements = m.mat_j in
    let ilen = Array.length i_elements in
    let jlen = Array.length j_elements in
    let f_begin i =
      p oc "\\begin{tabular}{|l|";
      for k = i to min (i + cols - 1) (ilen - 1) do p oc "r|" done;
      p oc "}\n\\hline\n";
      for k = i to min (i + cols - 1) (ilen - 1) do
        p oc "& %s" (string_of_i i_elements.(k))
      done;
      p oc "\\\\\n\\hline\n"
    in
    let f_end () =
      p oc "\\end{tabular}\n\n"
    in
    let rec iter i =
      if i < ilen then
        (
         f_begin i;
         for j = 0 to jlen - 1 do
           p oc "%s " (string_of_j j_elements.(j)) ;
           for k = i to min (i + cols - 1) (ilen - 1) do
             p oc "& %s" (string_of_v m.mat_mat.(k).(j))
           done;
           p oc "\\\\\n\\hline\n";
         done;
         f_end ();
         iter (i + cols)
        )
      else
        ()
    in
    iter 0

  let store_matrix f dm =
    let oc = open_out_bin f in
    output_value oc dm;
    close_out oc

  let load_matrix f =
    let ic = open_in_bin f in
    let v = input_value ic in
    close_in ic;
    v

  exception Uncompatible_matrices
  let concat_matrices_i ?check mat1 mat2 =
    if Array.length mat1.mat_j <> Array.length mat2.mat_j then
      raise Uncompatible_matrices;
    (
     match check with
       None -> ()
     | Some f -> Array.iteri (fun i v -> f v mat2.mat_j.(i)) mat1.mat_j
    );
    let len_i1 = Array.length mat1.mat_i in
    let len_i2 = Array.length mat2.mat_i in
    let len_i = len_i1 + len_i2 in
    let len_j = Array.length mat1.mat_j in
    let init_v =
      if len_j < 1 then
        None
      else
        if len_i1 < 1 then
          if len_i2 < 1 then
            None
          else
            Some mat2.mat_mat.(0).(0)
        else
          Some mat1.mat_mat.(0).(0)
    in
    let mat =
      match init_v with
        None -> [| |]
      | Some v ->
          let mat = Array.make_matrix len_i len_j v in
          for i = 0 to len_i1 - 1 do
            for j = 0 to len_j - 1 do
              mat.(i).(j) <- mat1.mat_mat.(i).(j)
            done
          done;
          for i = 0 to len_i2 - 1 do
            for j = 0 to len_j - 1 do
              mat.(len_i1 + i).(j) <- mat2.mat_mat.(i).(j)
            done
          done;
          mat
    in
    {
      mat_i = Array.append mat1.mat_i mat2.mat_i ;
      mat_j = mat1.mat_j ;
      mat_mat = mat;
    }

  let concat_matrices_j ?check mat1 mat2 =
    if Array.length mat1.mat_i <> Array.length mat2.mat_i then
      raise Uncompatible_matrices;
    (
     match check with
       None -> ()
     | Some f -> Array.iteri (fun i v -> f v mat2.mat_i.(i)) mat1.mat_i
    );
    let len_j1 = Array.length mat1.mat_j in
    let len_j2 = Array.length mat2.mat_j in
    let len_j = len_j1 + len_j2 in
    let len_i = Array.length mat1.mat_i in
    let init_v =
      if len_i < 1 then
        None
      else
        if len_j1 < 1 then
          if len_j2 < 1 then
            None
          else
            Some mat2.mat_mat.(0).(0)
        else
          Some mat1.mat_mat.(0).(0)
    in
    let mat =
      match init_v with
        None -> [| |]
      | Some v ->
          let mat = Array.make_matrix len_i len_j v in
          for j = 0 to len_j1 - 1 do
            for i = 0 to len_i - 1 do
              mat.(i).(j) <- mat1.mat_mat.(i).(j)
            done
          done;
          for j = 0 to len_j2 - 1 do
            for i = 0 to len_i - 1 do
              mat.(i).(len_j1 + j) <- mat2.mat_mat.(i).(j)
            done
          done;
          mat
    in
    {
      mat_i = mat1.mat_i ;
      mat_j = Array.append mat1.mat_j mat2.mat_j;
      mat_mat = mat;
    }
end