diff --git a/CI/CI.ml b/CI/CI.ml
index d2d258b..0f425b2 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -338,6 +338,73 @@ let create_matrix_spin f det_space =
 
 
 
+(* Create a matrix using the fact that the determinant space is made of
+   the outer product of spindeterminants.  *)
+let create_matrix_spin_computed f det_space = 
+  lazy (
+    let ndet = Ds.size det_space in
+    let a, b = 
+      match Ds.determinants det_space with
+      | Ds.Spin (a,b) -> (a,b)
+      | _ -> assert false
+    in
+    let n_beta = Array.length b in
+
+    let h i_alfa  =
+      let deg_a = Spindeterminant.degree i_alfa in
+      fun j_alfa ->
+        match deg_a j_alfa with
+        | 0 | 1 | 2 ->
+          (fun i_beta ->
+             let deg_b = Spindeterminant.degree i_beta in
+             let ki = Determinant.of_spindeterminants i_alfa i_beta in
+             fun j_beta ->
+               match deg_b j_beta with
+               | 0 | 1 | 2 -> (
+                   let kj = Determinant.of_spindeterminants j_alfa j_beta in
+                   f ki kj)
+               | _ -> 0.
+          )
+        | _ -> (fun _ _ -> 0.)
+    in
+
+    let i_prev = ref (-10)
+    and result = ref (fun _ -> 0.)
+    in
+
+    let g i =
+      if i <> !i_prev then
+        begin
+          i_prev := i;
+          let i_a = (i-1)/n_beta in
+          let i_alfa = i_a + 1 in
+          let h1 = 
+            h a.(i_alfa-1) 
+          in
+          let i_beta = i - i_a*n_beta  in
+          let bi = b.(i_beta-1) in
+          let h123_prev = ref (fun _ -> 0.) in
+          let j_alfa_prev = ref (-10) in
+          result := fun j ->
+              let j_a = (j-1)/n_beta in
+              let j_alfa = j_a + 1 in
+              let h123 = 
+                if j_alfa <> !j_alfa_prev then
+                  begin
+                    j_alfa_prev := j_alfa ;
+                    h123_prev := (h1 a.(j_alfa-1) bi)
+                  end;
+                !h123_prev
+              in 
+              let j_beta = j - j_a*n_beta  in
+              h123 b.(j_beta-1)
+        end;
+      !result
+    in
+
+    Matrix.of_fun ndet ndet g
+  )
+
 
 
 let make ?(n_states=1) ?(algo=`Direct) det_space = 
@@ -362,7 +429,11 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
     let f = 
       match Ds.determinants det_space with
       | Ds.Arbitrary _ -> create_matrix_arbitrary
-      | Ds.Spin      _ -> create_matrix_spin
+      | Ds.Spin      _ -> 
+        if algo = `Direct then
+          create_matrix_spin_computed
+        else 
+          create_matrix_spin
     in
     f (fun ki kj ->
          if ki <> kj then
@@ -393,9 +464,31 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
         ))
       in
       let matrix_prod psi = 
-        (*
-        Matrix.mm ~transa:`T m_H psi
-        *)
+        let result = 
+          Matrix.mm ~transa:`T m_H psi 
+        in
+        Parallel.broadcast (lazy result)
+      in
+      let eigenvectors, eigenvalues = 
+        let result = 
+          Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
+        in
+        Parallel.broadcast (lazy result)
+      in
+      let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
+      eigenvectors, eigenvalues
+    in
+
+    let eigensystem_direct () = 
+      let m_H =
+        Lazy.force m_H
+      in
+      let diagonal =
+        Parallel.broadcast (lazy (
+          Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
+        ))
+      in
+      let matrix_prod psi = 
         let result = 
           Matrix.parallel_mm ~transa:`T psi m_H
           |> Matrix.transpose
@@ -412,33 +505,6 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
       eigenvectors, eigenvalues
     in
 
-    let eigensystem_direct () = 
-      eigensystem_incore ()
-    in
-(*
-      let diagonal =
-        let stream = Ds.determinant_stream det_space in
-        Vec.init (Ds.size det_space) (fun _ -> 
-          let ki = Stream.next stream in
-          h_ij mo_basis ki ki -. e_shift)
-      in
-
-
-      let matrix_prod psi = 
-        (*TODO*)
-      in
-
-      let eigenvectors, eigenvalues = 
-        let result = 
-          Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
-        in
-        Parallel.broadcast (lazy result)
-      in
-      let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
-      eigenvectors, eigenvalues
-    in
-*)
-
     match algo with
     | `Direct -> eigensystem_direct ()
     | `InCore -> eigensystem_incore ()
diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml
index 2f2a827..64e4083 100644
--- a/Utils/Davidson.ml
+++ b/Utils/Davidson.ml
@@ -100,15 +100,23 @@ let make
 
     (* Compute the residual as proposed new vectors *)
     let u_proposed = 
-      Mat.init_cols n m (fun i k -> (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. 
-                              (max (diagonal.{i} -. lambda.{k}) 0.01) )
+      Mat.init_cols n m (fun i k ->
+        (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. 
+          (max (diagonal.{i} -. lambda.{k}) 0.01) )
+    in
+    let maxu = lange u_proposed ~norm:`M in
+    let thr = maxu *. 0.001 in
+    let u_proposed = 
+      Mat.map (fun x -> if abs_float x < thr then 0. else x) u_proposed
       |> Mat.to_col_vecs_list
     in
 
+
     let residual_norms = List.map nrm2 u_proposed in
     let residual_norm  = List.fold_left (fun accu i -> max accu i) 0. residual_norms in
 
-    Printf.printf "%3d %16.10f  %16.8e%!\n" iter lambda.{1} residual_norm;
+    if Parallel.master then
+      Printf.printf "%3d %16.10f  %16.8e%!\n" iter lambda.{1} residual_norm;
 
     if residual_norm > threshold then
       let u_next, w_next, iter =
diff --git a/Utils/Matrix.ml b/Utils/Matrix.ml
index 2b1b243..53f6c85 100644
--- a/Utils/Matrix.ml
+++ b/Utils/Matrix.ml
@@ -7,39 +7,56 @@ type sparse_matrix =
     v: Vector.t array;
   }
 
+type computed = 
+  {
+    m: int;
+    n: int;
+    f: int -> int -> float;
+  }
+
 type t =
 | Dense  of Mat.t
 | Sparse of sparse_matrix
+| Computed of computed
 
 let epsilon = Constants.epsilon
 
+let is_computed = function
+  | Computed _ -> true
+  | _ -> false
+
 let is_sparse = function
   | Sparse _ -> true
-  | Dense  _ -> false
-
+  | _ -> false
 
 let is_dense = function
-  | Sparse _ -> false
   | Dense  _ -> true
+  | _ -> false
 
 
 let dim1 = function
   | Dense  m -> Mat.dim1 m
-  | Sparse {m ; n ; v} -> m
+  | Sparse {m ; _} -> m
+  | Computed {m ; _} -> m
 
 
 let dim2 = function
   | Dense  m -> Mat.dim2 m
-  | Sparse {m ; n ; v} -> n
+  | Sparse {n ; _} -> n
+  | Computed {n ; _} -> n
 
 
+let check_bounds m n i j =
+  if (i <= 0 || i > m || j <= 0 || j > n) then
+     raise (Invalid_argument "Index out of bounds")
+
 let get = function
   | Dense  m -> (fun i j -> m.{i,j})
-  | Sparse {m ; n ; v } -> (fun i j -> Vector.get v.(j-1) i)
+  | Sparse { m ; n ; v } -> (fun i j -> Vector.get v.(j-1) i)
+  | Computed { m ; n ; f } -> (fun i j -> check_bounds m n i j ; f i j)
 
 
 let sparse_of_dense ?(threshold=epsilon) = function
-  | Sparse _ -> invalid_arg "Expected a dense matrix"
   | Dense m' ->
     let m = Mat.dim1 m'
     and n = Mat.dim2 m'
@@ -47,23 +64,41 @@ let sparse_of_dense ?(threshold=epsilon) = function
       Mat.to_col_vecs m'
       |> Array.map (fun v -> Vector.sparse_of_vec ~threshold v)
     in Sparse {m ; n ; v}
+  | _ -> invalid_arg "Expected a dense matrix"
 
 
 let dense_of_sparse = function
-  | Dense  _ -> invalid_arg "Expected a sparse matrix"
   | Sparse {m ; n ; v} -> 
     let m' =
       Array.map (fun v -> Vector.to_vec v) v
       |> Mat.of_col_vecs
     in Dense m'
+  | _ -> invalid_arg "Expected a sparse matrix"
 
 
+let sparse_of_computed ?(threshold=epsilon) = function
+  | Computed {m ; n ; f} ->
+    Sparse { m ; n ; v=Array.init n (fun j ->
+        Util.list_range 1 m
+        |> List.map (fun i ->
+            let x = f i (j+1) in
+            if abs_float x > threshold then Some (i, x)
+            else None)
+        |> Util.list_some
+        |> Vector.sparse_of_assoc_list m
+      ) }
+  | _ -> invalid_arg "Expected a computed matrix"
+
+let dense_of_computed x = dense_of_sparse @@ sparse_of_computed x
+
 let dense_of_mat m = Dense m
 
+let of_fun m n f = Computed {m ; n ; f}
 
 let rec to_vector_array ?(threshold=epsilon) = function
   | Sparse {m ; n ; v} -> v
   | Dense m -> to_vector_array (sparse_of_dense ~threshold (Dense m))
+  | Computed m -> to_vector_array (sparse_of_computed ~threshold (Computed m))
 
 
 let identity n = 
@@ -92,6 +127,7 @@ let rec to_mat = function
   | Sparse m -> 
     dense_of_sparse (Sparse m)
     |> to_mat 
+  | Computed m -> sparse_of_computed (Computed m) |> dense_of_sparse |> to_mat
 
 let transpose = function
   | Dense  m -> Dense (Mat.transpose_copy m)
@@ -109,6 +145,9 @@ let transpose = function
       in
       Sparse {m=n ; n=m ; v=v'}
     end
+  | Computed {m ; n ; f} ->
+     let f' i j = f j i in
+     Computed { m=n ; n=m ; f=f' }
 
 
 let outer_product ?(threshold=epsilon) v1 v2 = 
@@ -152,7 +191,7 @@ let outer_product ?(threshold=epsilon) v1 v2 =
 
 
 
-let mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b = 
+let rec mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b = 
 
   let f, f' = 
     match transa, transb with
@@ -259,11 +298,60 @@ let mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b =
     Sparse {m=m' ; n=n ; v=v''}
   in
 
+  let mmcc transa transb a b =
+    let {m ; n ; f} =
+      if transb = `T then
+        match transpose (Computed b) with
+        | Computed x -> x
+        | _ -> assert false
+      else
+        b
+    in
+    let m', n', f' =
+      if transa = `T then
+        match transpose (Computed a) with
+        | Computed {m ; n ; f} -> m, n, f
+        | _ -> assert false
+      else
+        let {m ; n ; f} = a in
+        m, n, f
+    in
+    if n' <> m then
+      invalid_arg "Inconsistent dimensions";
+    let g i j =
+      let result = ref 0. in
+      for k=1 to m do
+        let  a = f k j in
+        if a <> 0. then
+          result := !result +. (f' i k) *. a ;
+      done;
+      !result
+    in
+    Computed {m=m' ; n=n ; f=g}
+  in
+      
   match a, b with
   | (Dense  a), (Dense  b) -> Dense (gemm ~transa ~transb a b)
   | (Sparse a), (Dense  b) -> spmm transa transb a b
   | (Dense  a), (Sparse b) -> mmsp transa transb a b
   | (Sparse a), (Sparse b) -> mmspmm transa transb a b
+  | (Computed a), (Computed b) -> mmcc transa transb a b
+  | (Computed a), (Dense    _) ->
+      let b = { m = dim1 b ; n = dim2 b ; f = get b } in
+      mmcc transa transb a b
+      |> dense_of_computed
+  | (Computed a), (Sparse   _) ->
+      let b = { m = dim1 b ; n = dim2 b ; f = get b } in
+      mmcc transa transb a b
+      |> sparse_of_computed ~threshold
+  | _, (Computed _) ->
+     begin
+       match transa, transb with
+       | `N, `N -> mm ~transa:`T ~transb:`T ~threshold b a
+       | `N, `T -> mm ~transa:`N ~transb:`T ~threshold b a
+       | `T, `N -> mm ~transa:`T ~transb:`N ~threshold b a
+       | `T, `T -> mm ~transa:`N ~transb:`N ~threshold b a
+     end |> transpose
 
 
 let mv ?(sparse=false) ?(trans=`N) ?(threshold=epsilon) a b =
@@ -300,12 +388,28 @@ let mv ?(sparse=false) ?(trans=`N) ?(threshold=epsilon) a b =
         |> Vector.dot b )
   in
 
+  let cmv a b =
+    match trans with
+    | `N -> Vec.init a.m (fun i -> 
+         let accu = ref 0. in
+         for j=1 to a.n do
+           accu := !accu +. a.f i j *. Vector.get b j
+         done;
+         !accu)
+    | `T -> Vec.init a.m (fun i -> 
+         let accu = ref 0. in
+         for j=1 to a.n do
+           accu := !accu +. a.f j i *. Vector.get b j
+         done;
+         !accu)
+  in
+
   let dense_result =
     match a, Vector.is_dense b with
     | Dense  a, true   -> gemv ~trans a (Vector.to_vec b)
     | Dense  a, false  -> mv a b
-    | Sparse a, true   -> spmv a b
-    | Sparse a, false  -> spmv a b
+    | Sparse a, _ -> spmv a b
+    | Computed a, _ -> cmv a b
   in
 
   if sparse then
@@ -326,6 +430,7 @@ let rec op2 dense_op sparse_op a b =
                                 { m=a.m ; n=a.n ;
                                   v = Array.map2 sparse_op a.v b.v
                                 }
+  | _ -> failwith "Not implemented"
 
 let add = op2 (fun a b -> Mat.add a b) (fun a b -> Vector.add a b)
 let sub = op2 (fun a b -> Mat.sub a b) (fun a b -> Vector.sub a b)
@@ -336,12 +441,14 @@ let scale f = function
                   { a with
                     v = if f = 1.0 then a.v
                       else Array.map (fun v -> Vector.scale f v) a.v }
+  | _ -> failwith "Not implemented"
 
 let frobenius_norm = function
   | Dense  a -> lange ~norm:`F a
   | Sparse a -> 
     Array.fold_left (fun accu v -> accu +. Vector.dot v v) 0. a.v
     |> sqrt
+  | _ -> failwith "Not implemented"
 
 
 
@@ -361,6 +468,13 @@ let split_cols nrows = function
       |> List.map Array.of_list
       |> List.map (fun v -> Sparse { m=a.m ; n= Array.length v ; v })
     end
+  | Computed a -> 
+    begin
+      Util.list_range 0 (a.n-1)
+      |> Util.list_pack nrows
+      |> List.map Array.of_list
+      |> List.map (fun v -> Computed { m=a.m ; n= Array.length v ; f = (fun i j -> a.f i (j+v.(0)) ) })
+    end
   
 
 let join_cols l =
@@ -376,8 +490,9 @@ let join_cols l =
 
   and aux = function
     | [] -> Sparse { m=0 ; n=0 ; v=[| |] }
-    | (Dense  a) :: rest -> aux_dense  [] ((Dense  a) :: rest)
-    | (Sparse a) :: rest -> aux_sparse 0 0 [] ((Sparse a) :: rest)
+    | (Dense    a) :: rest -> aux_dense  [] ((Dense  a) :: rest)
+    | (Sparse   a) :: rest -> aux_sparse 0 0 [] ((Sparse a) :: rest)
+    | (Computed a) :: rest -> aux_sparse 0 0 [] (List.map sparse_of_computed ( (Computed a) :: rest ))
 
   in aux (List.rev l)
 
@@ -440,11 +555,21 @@ let rec ax_eq_b ?(trans=`N) a b =
 
 let parallel_mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b = 
 
-  let n = 4 in
+  let n = 
+    match transa with
+    | `N -> dim2 a
+    | `T -> dim1 a
+  in
+  let n = n / (Parallel.size * 4) in
   split_cols n b
   |> Stream.of_list 
   |> Farm.run ~ordered:true ~f:(fun b ->
-            mm ~transa ~transb ~threshold a b
+            match a, b with
+              | Computed _, Computed _ -> 
+                  mm ~transa ~transb ~threshold a b
+                  |> sparse_of_computed ~threshold
+              | _ ->
+                  mm ~transa ~transb ~threshold a b
             )
   |> Util.stream_to_list
   |> join_cols
@@ -453,8 +578,9 @@ let parallel_mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b =
 (* ------------ Printers ------------ *)
 
 let rec pp_matrix ppf = function
-  | Dense  m -> Util.pp_matrix ppf m
-  | Sparse m -> pp_matrix ppf @@ dense_of_sparse (Sparse m)
+  | Dense    m -> Util.pp_matrix ppf m
+  | Sparse   m -> pp_matrix ppf @@ dense_of_sparse (Sparse m)
+  | Computed m -> pp_matrix ppf @@ dense_of_computed (Computed m)
 
 
 
@@ -593,25 +719,40 @@ let test_case () =
     and m5   = dense_of_mat  x2 |> transpose
     and m5_s = sparse_of_mat x2 |> transpose
     in
+    let c1 = of_fun (Mat.dim1 x1) (Mat.dim2 x1) (fun i j -> x1.{i,j}) in
+    let c2 = of_fun (Mat.dim1 x2) (Mat.dim2 x2) (fun i j -> x2.{i,j}) in
+    let c3 = of_fun (Mat.dim1 x3) (Mat.dim2 x3) (fun i j -> x3.{i,j}) in
+    let c4 = of_fun (dim1 m4) (dim2 m4) (fun i j -> get m4 i j ) in
+    let c5 = of_fun (dim1 m5) (dim2 m5) (fun i j -> get m5 i j ) in
+    Alcotest.(check (float 1.e-10)) "dense  dense   0" 0. (norm_diff m3 c3);
     Alcotest.(check (float 1.e-10)) "dense  dense   1" 0. (norm_diff (mm m1 m2) m3);
-    Alcotest.(check (float 1.e-10)) "dense  dense   2" 0. (norm_diff (mm ~transa:`T m4 m2) m3);
-    Alcotest.(check (float 1.e-10)) "dense  dense   3" 0. (norm_diff (mm ~transb:`T m1 m5) m3);
-    Alcotest.(check (float 1.e-10)) "dense  dense   4" 0. (norm_diff (mm ~transa:`T ~transb:`T m2 m1) (transpose m3));
+    Alcotest.(check (float 1.e-10)) "dense  dense   2" 0. (norm_diff (mm c1 c2) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   3" 0. (norm_diff (mm c1 m2) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   4" 0. (norm_diff (mm c1 m2_s) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   5" 0. (norm_diff (mm m1 c2) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   6" 0. (norm_diff (mm m1_s c2) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   7" 0. (norm_diff (mm ~transa:`T m4 m2) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   8" 0. (norm_diff (mm ~transa:`T c4 m2) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   9" 0. (norm_diff (mm ~transb:`T m1 m5) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   10" 0. (norm_diff (mm ~transb:`T m1 c5) m3);
+    Alcotest.(check (float 1.e-10)) "dense  dense   11" 0. (norm_diff (mm ~transa:`T ~transb:`T m2 m1) (transpose m3));
                                                    
-    Alcotest.(check (float 1.e-10)) "dense  sparse  5" 0. (norm_diff (mm m1 m2_s) m3_s);
-    Alcotest.(check (float 1.e-10)) "dense  sparse  6" 0. (norm_diff (mm ~transa:`T m4 m2_s) m3_s);
-    Alcotest.(check (float 1.e-10)) "dense  sparse  7" 0. (norm_diff (mm ~transb:`T m1 m5_s) m3_s);
-    Alcotest.(check (float 1.e-10)) "dense  sparse  8" 0. (norm_diff (transpose (mm m2 m1_s ~transa:`T ~transb:`T)) m3_s);
+    Alcotest.(check (float 1.e-10)) "dense  sparse  12" 0. (norm_diff (mm m1 m2_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "dense  sparse  13" 0. (norm_diff (mm ~transa:`T m4 m2_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "dense  sparse  14" 0. (norm_diff (mm ~transa:`T c4 m2_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "dense  sparse  15" 0. (norm_diff (mm ~transb:`T m1 m5_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "dense  sparse  16" 0. (norm_diff (transpose (mm m2 m1_s ~transa:`T ~transb:`T)) m3_s);
 
-    Alcotest.(check (float 1.e-10)) "sparse dense   9" 0. (norm_diff (mm m1_s m2) m3_s);
-    Alcotest.(check (float 1.e-10)) "sparse dense  10" 0. (norm_diff (mm ~transa:`T m4_s m2) m3_s);
-    Alcotest.(check (float 1.e-10)) "sparse dense  11" 0. (norm_diff (mm ~transb:`T m1_s m5) m3_s);
-    Alcotest.(check (float 1.e-10)) "sparse dense  12" 0. (norm_diff (transpose (mm m2_s m1 ~transa:`T ~transb:`T)) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse dense   17" 0. (norm_diff (mm m1_s m2) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse dense   18" 0. (norm_diff (mm ~transa:`T m4_s m2) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse dense   19" 0. (norm_diff (mm ~transb:`T m1_s m5) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse dense   20" 0. (norm_diff (mm ~transb:`T m1_s c5) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse dense   21" 0. (norm_diff (transpose (mm m2_s m1 ~transa:`T ~transb:`T)) m3_s);
 
-    Alcotest.(check (float 1.e-10)) "sparse sparse 13" 0. (norm_diff (mm m1_s m2_s) m3_s);
-    Alcotest.(check (float 1.e-10)) "sparse sparse 14" 0. (norm_diff (mm ~transa:`T m4_s m2_s) m3_s);
-    Alcotest.(check (float 1.e-10)) "sparse sparse 15" 0. (norm_diff (mm ~transb:`T m1_s m5_s) m3_s);
-    Alcotest.(check (float 1.e-10)) "sparse sparse 16" 0. (norm_diff (transpose (mm m2_s m1_s ~transa:`T ~transb:`T)) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse sparse  22" 0. (norm_diff (mm m1_s m2_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse sparse  23" 0. (norm_diff (mm ~transa:`T m4_s m2_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse sparse  24" 0. (norm_diff (mm ~transb:`T m1_s m5_s) m3_s);
+    Alcotest.(check (float 1.e-10)) "sparse sparse  25" 0. (norm_diff (transpose (mm m2_s m1_s ~transa:`T ~transb:`T)) m3_s);
   in
 
   let test_solve () = 
diff --git a/Utils/Matrix.mli b/Utils/Matrix.mli
index d81ab7b..cbf5a0c 100644
--- a/Utils/Matrix.mli
+++ b/Utils/Matrix.mli
@@ -30,9 +30,19 @@ val to_mat : t -> Mat.t
 val to_vector_array : ?threshold:float -> t -> Vector.t array
 (** Convert the matrix into an array of column vectors. *)
 
+val of_fun : int -> int -> (int -> int -> float) -> t
+(** [of_fun m n f] creates a computed matrix of dimension m times n, where the function 
+    [f i j] is called to evaluate element [i j] *)
+
 val sparse_of_dense : ?threshold:float -> t -> t
 (** Creates a sparse matrix from a dense matrix. Default threshold is {!Constants.epsilon}. *)
 
+val sparse_of_computed : ?threshold:float -> t -> t
+(** Creates a sparse matrix from a computed matrix. *)
+
+val dense_of_computed : t -> t
+(** Creates a dense matrix from a computed matrix. *)
+
 val dense_of_sparse : t -> t
 (** Creates a dense matrix from a sparse matrix. *)