Working on parallelism

Basis/ERI.ml

@@ -314,11 +314,9 @@ let of_basis_parallel basis =
   |> Stream.iter (fun l ->
      Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
       set_chem eri_array i_c j_c k_c l_c value) l);
-  if not Parallel.master then
-    exit 0;
 
   Printf.printf "Computed ERIs in parallel in %f seconds\n%!" (Unix.gettimeofday () -. t0);
-  eri_array
+  Parallel.broadcast (lazy eri_array)
 
 
 

CI/CI.ml

@@ -53,14 +53,73 @@ let make det_space =
   let ndet = Ds.size det_space in
   let det  = Ds.determinants det_space in
   let mo_basis = Ds.mo_basis det_space in
+  (*
   let h_matrix = lazy (
-    Array.init ndet (fun i -> let ki = det.(i) in
+    Util.list_range 0 (ndet-1)
+    |> List.map (fun i -> let ki = det.(i) in
                   Array.init ndet (fun j -> let kj = det.(j) in
                                     h_ij mo_basis ki kj
-                                      ))
-    |> Mat.of_array
+                                  )
+                  |> Vec.of_array)
+    |> Mat.of_col_vecs_list
   )
   in
+  *)
+  (*
+  let h_matrix = lazy (
+    let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in
+    List.init ntasks (fun i -> 
+      let m = 
+        if i = (ntasks-1) then
+          (ndet - i*ntasks) 
+        else
+          ntasks
+      in
+      List.init m (fun j -> i*ntasks + j)
+      )
+    |> Stream.of_list
+    |> Farm.run ~ordered:true
+      ~f:(fun l ->
+          Printf.eprintf "Start\n%!";
+      List.map (fun i ->
+          let ki = det.(i) in
+          Printf.eprintf "%d / %d\n%!" i ndet;
+          Array.init ndet (fun j -> let kj = det.(j) in
+                            h_ij mo_basis ki kj)
+          ) l)
+    |> Util.stream_to_list
+    |> List.concat
+    |> List.map Vec.of_array
+    |> Mat.of_col_vecs_list
+  ) in
+  *)
+  let h_matrix = lazy (
+    let h =
+      if Parallel.master then
+        Array.make_matrix ndet ndet 0.
+        |> Mat.of_array
+      else
+        Array.make_matrix 1 1 0.
+        |> Mat.of_array
+    in
+    List.init ndet (fun i -> i)
+    |> Stream.of_list
+    |> Farm.run ~ordered:false
+      ~f:(fun i ->
+          let ki = det.(i) in
+          Printf.eprintf "%d / %d\n%!" i ndet;
+          List.init ndet (fun j ->
+            let kj = det.(j) in
+            let x = h_ij mo_basis ki kj in
+            if x <> 0. then Some (i,j,x) else None
+          )
+          |> Util.list_some
+          )
+    |> Util.stream_to_list
+    |> List.iter (fun l ->  if Parallel.master then
+      List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l);
+    Parallel.broadcast (lazy h)
+  ) in
   let s2_matrix = lazy (
     Array.init ndet (fun i -> let ki = det.(i) in
                       Array.init ndet (fun j -> let kj = det.(j) in
@@ -70,8 +129,9 @@ let make det_space =
   )
   in
   let eigensystem = lazy (
-    Lazy.force h_matrix
-    |> Util.diagonalize_symm 
+    let h_matrix = Lazy.force h_matrix in
+    Parallel.broadcast @@
+      lazy (Util.diagonalize_symm h_matrix)
   )
   in
   { det_space ; h_matrix ; s2_matrix ; eigensystem }

INSTALL.md

@@ -50,6 +50,14 @@ cd odoc-ltxhtml
 make install
 ```
 
+# Getopt
+
+Parsing of command line arguments (similar to GNU GetOpt)
+
+```bash
+opam install getopt
+```
+
 # Alcotest
 
 Lightweight and colourful test framework

MOBasis/MOBasis.ml

@@ -172,6 +172,7 @@ let of_rhf hf =
     make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
   in
 
+(*
   let () =
     let e = ref 0. in
     let t = KinInt.matrix (Lazy.force result.kin_ints) in
@@ -196,6 +197,7 @@ let of_rhf hf =
     Printf.printf "Energy two-e = %20.15f\n" !e2;
     Printf.printf "Energy       = %20.15f\n" (Si.nuclear_repulsion simulation +. !e +. !e2)
   in
+  *)
   result
 
 

SCF/RHF.ml

@@ -140,7 +140,9 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
       in
 
       let gap =
-        eigenvalues.{nocc+1} -. eigenvalues.{nocc};
+        if nocc < Vec.dim eigenvalues then
+          eigenvalues.{nocc+1} -. eigenvalues.{nocc} 
+        else 0.
       in
 
       let () = 

Utils/Command_line.ml

@@ -0,0 +1,199 @@
+type short_opt     = char
+type long_opt      = string 
+type optional      = Mandatory | Optional 
+type documentation = string
+type argument      = With_arg of string | Without_arg | With_opt_arg of string
+
+type description   = {
+  short: short_opt ;
+  long : long_opt  ;
+  opt  : optional  ;
+  doc  : documentation ;
+  arg  : argument ;
+  }
+
+let anon_args  = ref []
+and header_doc = ref ""
+and description_doc = ref ""
+and footer_doc = ref ""
+and specs      = ref []
+
+let set_header_doc s = header_doc := s
+let set_description_doc s = description_doc := s
+let set_footer_doc s = footer_doc := s
+
+(* Hash table containing all the options *)
+let dict = Hashtbl.create 67 
+
+let get_bool x = Hashtbl.mem dict x
+
+let show_help () = get_bool "help"
+
+let get x =
+    try Some (Hashtbl.find dict x)
+    with Not_found -> None
+
+let anonymous name opt doc =
+  { short=' ' ; long=name; opt; doc; arg=Without_arg; }
+
+let output_text t =
+  Format.printf "@[<v 0>";
+  begin
+    match Str.split (Str.regexp "\n") t with
+    | x :: [] -> Format.printf "@[<hov 0>";
+                    Str.split (Str.regexp " ") x
+                    |> List.iter (fun y -> Format.printf "@[%s@]@ " y) ;
+                    Format.printf "@]"
+    | t       -> List.iter (fun x -> 
+                    Format.printf "@[<hov 0>";
+                    Str.split (Str.regexp " ") x
+                    |> List.iter (fun y -> Format.printf "@[%s@]@ " y) ;
+                    Format.printf "@]@;") t
+  end;
+  Format.printf "@]"
+;;
+  
+
+let output_short x =
+  match x.short, x.opt, x.arg  with
+  | ' ', Mandatory, _ -> Format.printf "@[%s@]"   x.long
+  | ' ', Optional , _ -> Format.printf "@[[%s]@]" x.long
+  |  _ , Mandatory, Without_arg -> Format.printf "@[-%c@]"   x.short
+  |  _ , Optional , Without_arg -> Format.printf "@[[-%c]@]" x.short
+  |  _ , Mandatory, With_arg arg -> Format.printf "@[-%c %s@]"   x.short arg
+  |  _ , Optional , With_arg arg -> Format.printf "@[[-%c %s]@]" x.short arg
+  |  _ , Mandatory, With_opt_arg arg -> Format.printf "@[-%c [%s]@]"   x.short arg
+  |  _ , Optional , With_opt_arg arg -> Format.printf "@[[-%c [%s]]@]" x.short arg
+
+
+let output_long max_width x =
+  let arg =
+    match x.short, x.arg with
+    | ' ' , _                -> x.long
+    |  _  , Without_arg      -> x.long
+    |  _  , With_arg     arg -> Printf.sprintf "%s=%s" x.long arg
+    |  _  , With_opt_arg arg -> Printf.sprintf "%s[=%s]" x.long arg
+  in
+  let long =
+     let l = String.length arg in
+     arg^(String.make (max_width-l) ' ')
+  in
+  Format.printf "@[<v 0>";
+  begin
+    match x.short with
+    | ' '   -> Format.printf "@[%s       @]" long
+    | short -> Format.printf "@[-%c  --%s @]" short long
+  end;
+  Format.printf "@]";
+  output_text  x.doc
+
+
+let help () =
+
+    (* Print the header *)
+    output_text !header_doc;
+    Format.printf "@.@.";
+
+    (* Find the anonymous arguments *)
+    let anon =
+      List.filter (fun x -> x.short = ' ') !specs
+    in
+
+    (* Find the options *)
+    let options =
+      List.filter (fun x -> x.short <> ' ') !specs
+      |> List.sort (fun x y -> Char.compare x.short y.short)
+    in
+
+    (* Find column lengths *)
+    let max_width =
+      List.map (fun x -> 
+        ( match x.arg with
+          | Without_arg      -> String.length x.long
+          | With_arg arg     -> String.length x.long + String.length arg
+          | With_opt_arg arg -> String.length x.long + String.length arg + 2
+        )
+        + ( if x.opt = Optional then 2 else 0)
+      ) !specs
+      |> List.fold_left max 0
+    in
+
+
+    (* Print usage *)
+    Format.printf "@[<v>@[<v 2>Usage:@,@,@[<hov 4>@[%s@]" Sys.argv.(0);
+    List.iter (fun x -> Format.printf "@ "; output_short x) options;
+    Format.printf "@ @[[--]@]";
+    List.iter (fun x -> Format.printf "@ "; output_short x;) anon;
+    Format.printf "@]@,@]@,";
+
+
+    (* Print arguments and doc *)
+    Format.printf "@[<v 2>Arguments:@,";
+    Format.printf "@[<v 0>" ;
+    List.iter (fun x -> Format.printf "@ "; output_long max_width x) anon;
+    Format.printf "@]@,@]@,";
+
+
+    (* Print options and doc *)
+    Format.printf "@[<v 2>Options:@,";
+
+    Format.printf "@[<v 0>" ;
+    List.iter (fun x -> Format.printf "@ "; output_long max_width x) options;
+    Format.printf "@]@,@]@,";
+
+
+    (* Print footer *)
+    if !description_doc <> "" then
+    begin
+      Format.printf "@[<v 2>Description:@,@,";
+      output_text !description_doc;
+      Format.printf "@,"
+    end;
+
+    (* Print footer *)
+    output_text !footer_doc;
+    Format.printf "@."
+
+
+
+let set_specs specs_in =
+    specs := { short='h' ;
+               long ="help" ;
+               doc  ="Prints the help message." ;
+               arg  =Without_arg ;
+               opt  =Optional ;
+             } :: specs_in;
+
+    let cmd_specs =
+      List.filter (fun x -> x.short != ' ') !specs
+      |> List.map (fun { short ; long ; opt ; doc ; arg } ->
+                  match arg with
+                  |  With_arg _ ->
+                      (short, long, None, Some (fun x -> Hashtbl.replace dict long x) )
+                  |  Without_arg ->
+                      (short, long, Some (fun () -> Hashtbl.replace dict long ""), None)
+                  |  With_opt_arg _ ->
+                      (short, long, Some (fun () -> Hashtbl.replace dict long ""),
+                      Some (fun x -> Hashtbl.replace dict long x) )
+                  )
+    in
+
+    Getopt.parse_cmdline cmd_specs (fun x -> anon_args := !anon_args @ [x]);
+
+    if show_help () then
+        (help () ; exit 0);
+
+    (* Check that all mandatory arguments are set *)
+    List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
+    |> List.iter (fun x -> 
+        match get x.long with
+        | Some _ -> ()
+        | None -> failwith ("Error: --"^x.long^" option is missing.")
+        )
+;;
+
+
+let anon_args () = !anon_args
+
+
+

Utils/Util.ml

@@ -98,7 +98,7 @@ let fact = function
 
 let binom n k =
   (*TODO : slow function *)
-  assert (n > k);
+  assert (n >= k);
   let rec aux n k =
     if k = 0 || k = n then
       1
@@ -188,8 +188,15 @@ let stream_range first last =
     else None
   )
 
+let stream_to_list stream = 
+  let rec aux accu =
+    try aux (Stream.next stream :: accu) with
+    Stream.Failure -> List.rev accu
+  in aux []
+
 
 let list_range first last = 
+  if last < first then [] else
   let rec aux accu = function
   | 0 -> first :: accu
   | i -> aux ( (first+i)::accu ) (i-1)

Utils/Util.mli

@@ -73,6 +73,8 @@ val list_range : int -> int -> int list
 val stream_range : int -> int -> int Stream.t
 (** [stream_range first last] returns a stream <first ; first+1 ; ... ; last-1 ; last>. *)
 
+val stream_to_list : 'a Stream.t -> 'a list
+(** Read a stream and put items in a list. *)
 
 (** {2 Linear algebra } *)