Faster determinants in OCaml

configure

@@ -438,18 +438,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
         fail
 fi
 
-BWRAP=$(find_exe bwrap)
-if [[ ${BWRAP} = $(not_found) ]] ; then
-        error "Bubblewrap (bwrap) is not installed."
-        fail
-fi
-
 LIBCAP=$(find_lib -lcap)
 if [[ ${LIBCAP} = $(not_found) ]] ; then
         error "Libcap (libcap) is not installed."
         fail
 fi
 
+BWRAP=$(find_exe bwrap)
+if [[ ${BWRAP} = $(not_found) ]] ; then
+        error "Bubblewrap (bwrap) is not installed."
+        fail
+fi
+
 OPAM=$(find_exe opam)
 if [[ ${OPAM} = $(not_found) ]] ; then
         error "OPAM (ocaml) package manager is not installed."

ocaml/Bitlist.ml

@@ -7,82 +7,61 @@ Type for bits strings
 list of Bits
 *)
 
-type t = Bit.t list
+type t = int64 array
+
+
+let n_int = Array.length
+
+(* Create a zero bit list *)
+let zero n_int =
+  Array.make (N_int_number.to_int n_int) 0L
 
 (* String representation *)
 let to_string b =
-  let rec do_work accu = function
+  let int64_to_string x = 
-    | [] -> accu
+    String.init 64 (fun i ->
-    | head :: tail ->
+      if Int64.logand x @@ Int64.shift_left 1L i <> 0L then
-      let new_accu = (Bit.to_string head) ^ accu
+        '+'
-      in do_work new_accu tail
+      else
+        '-')
   in
-  do_work "" b
+  Array.map int64_to_string b
+  |> Array.to_list
+  |> String.concat ""
 
 
 let of_string ?(zero='0') ?(one='1') s =
-  List.init (String.length s) (String.get s)
+  let n_int = ( (String.length s - 1) lsr 6 ) + 1 in
-  |> List.rev_map ( fun c ->
+  let result = Array.make n_int 0L in
-     if (c = zero) then Bit.Zero
+  String.iteri (fun i c ->
-     else if (c = one) then Bit.One
+    if c = one then
-     else (failwith ("Error in bitstring ") ) )
+      begin
+        let iint = i lsr 6 in (* i / 64 *)
+        let k = i - (iint lsl 6) in 
+        result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k;
+      end) s;
+  result
 
-let of_string_mp s =
+let of_string_mp = of_string ~zero:'-' ~one:'+'
-  List.init (String.length s) (String.get s)
-  |> List.rev_map (function
-    | '-' -> Bit.Zero
-    | '+' -> Bit.One
-    | _   -> failwith ("Error in bitstring ") )
 
 
 (* Create a bit list from an int64 *)
-let of_int64 i =
+let of_int64 i = [| i |]
-
-  let rec do_work accu = function
-  | 0L -> Bit.Zero :: accu |> List.rev
-  | 1L -> Bit.One  :: accu |> List.rev
-  | i  ->
-    let b =
-      match (Int64.logand i 1L ) with
-      | 0L -> Bit.Zero
-      | 1L -> Bit.One
-      | _ -> raise (Failure "i land 1 not in (0,1)")
-    in
-    do_work (b :: accu) (Int64.shift_right_logical i 1)
-  in
-
-  let adjust_length result =
-    let rec do_work accu = function
-    | 64 -> List.rev accu
-    | i when i>64 -> raise (Failure "Error in of_int64 > 64")
-    | i when i<0 -> raise (Failure "Error in of_int64 < 0")
-    | i -> do_work (Bit.Zero :: accu)  (i+1)
-    in
-    do_work (List.rev result) (List.length result)
-  in
-  adjust_length (do_work [] i)
-
 
 (* Create an int64 from a bit list *)
-let to_int64 l =
+let to_int64 = function
-  assert ( (List.length l) <= 64) ;
+  | [| i |] -> i
-  let rec do_work accu = function
+  | _ -> failwith "N_int > 1"
-    | [] -> accu
+
-    | Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail
+
-    | Bit.One::tail  -> do_work Int64.(logor one (shift_left accu 1)) tail
+(* Create a bit list from an array of int64 *)
-  in do_work Int64.zero (List.rev l)
+external of_int64_array : int64 array -> t = "%identity"
+external to_int64_array : t -> int64 array = "%identity"
 
  
 (* Create a bit list from a list of int64 *)
 let of_int64_list l =
-  List.map of_int64 l
+  Array.of_list l |> of_int64_array
-  |> List.concat
-
-(* Create a bit list from an array of int64 *)
-let of_int64_array l =
-  Array.map of_int64 l
-  |> Array.to_list
-  |> List.concat
 
 
 (* Compute n_int *)
@@ -91,100 +70,63 @@ let n_int_of_mo_num mo_num =
   N_int_number.of_int ( (mo_num-1)/bit_kind_size + 1 )
 
 
-(* Create a zero bit list *)
-let zero n_int =
-  let n_int = N_int_number.to_int n_int in
-  let a = Array.init n_int (fun i-> 0L)  in
-  of_int64_list ( Array.to_list a )
-
 
 (* Create an int64 list from a bit list *)
 let to_int64_list l =
-  let rec do_work accu buf counter = function
+  to_int64_array l |> Array.to_list
-    | [] ->
-        begin
-          match buf with
-          | [] -> accu
-          | _  -> (List.rev buf)::accu
-        end
-    | i::tail ->
-        if (counter < 64) then
-          do_work accu (i::buf) (counter+1) tail
-        else
-          do_work ( (List.rev (i::buf))::accu) [] 1 tail
-  in
-  let l = do_work [] [] 1 l
-  in
-  List.rev_map to_int64 l
 
-(* Create an array of int64 from a bit list *)
-let to_int64_array l =
-    to_int64_list l
-    |> Array.of_list
 
 (* Create a bit list from a list of MO indices *)
 let of_mo_number_list n_int l =
-  let n_int = N_int_number.to_int n_int in
+  let result = zero n_int in
-  let length = n_int*64 in
+  List.iter (fun j ->
-  let a = Array.make length (Bit.Zero) in
+        let i = (MO_number.to_int j) - 1 in
-  List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l;
+        let iint = i lsr 6 in (* i / 64 *)
-  Array.to_list a
+        let k = i - (iint lsl 6) in
+        result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k;
+      ) l;
+  result
 
 
 let to_mo_number_list l =
-  let a = Array.of_list l in
+  let rec aux_one x shift accu = function
-  let mo_num = MO_number.get_max () in
+    | -1 -> accu
-  let rec do_work accu = function
+    | i  -> if Int64.logand x (Int64.shift_left 1L i) <> 0L then
-  | 0 -> accu
+              aux_one x shift ( (i+shift) ::accu) (i-1)
-  | i ->
+            else
-      begin
+              aux_one x shift accu (i-1)
-        let new_accu =
-        match a.(i-1) with
-        | Bit.One  -> (MO_number.of_int ~max:mo_num i)::accu
-        | Bit.Zero -> accu
   in
-        do_work new_accu (i-1)
+  Array.mapi (fun i x ->
-      end
+    let shift = (i lsr 6) lsl 6 + 1 in
-  in
+    aux_one x shift [] 63
-  do_work [] (List.length l)
+    ) l
-
+  |> Array.to_list
+  |> List.concat
+  |> List.map MO_number.of_int
 
 
 
 (* logical operations on bit_list *)
-let logical_operator2 op a b =
+let and_operator a b = Array.map2 Int64.logand a b
-  let rec do_work_binary result a b =
+let xor_operator a b = Array.map2 Int64.logxor a b
-  match a, b with
+let  or_operator a b = Array.map2 Int64.logor  a b
-  | [], [] -> result
+let not_operator   b = Array.map  Int64.lognot   b
-  | [], _  | _ , [] -> raise (Failure "Lists should have same length")
+
-  | (ha::ta), (hb::tb) ->
+
-    let newbit = op ha hb
+
-    in  do_work_binary (newbit::result) ta tb
+let pop_sign = 
+  let mask =
+      (Int64.pred (Int64.shift_left 1L 63))
   in
-  List.rev (do_work_binary [] a b)
+  fun x -> Int64.logand mask x
-
-
-let logical_operator1 op b =
-  let rec do_work_unary result b =
-  match b with
-  | [] -> result
-  | (hb::tb) ->
-    let newbit = op hb
-    in  do_work_unary (newbit::result) tb
-  in
-  List.rev (do_work_unary [] b)
-
-
-let and_operator a b = logical_operator2 Bit.and_operator a b
-let xor_operator a b = logical_operator2 Bit.xor_operator a b
-let  or_operator a b = logical_operator2  Bit.or_operator a b
-let not_operator   b = logical_operator1 Bit.not_operator   b
 
 
 let popcnt b =
-  List.fold_left (fun accu -> function
+  Array.fold_left (fun accu x ->
-    | Bit.One -> accu+1
+    if x >= 0L then
-    | Bit.Zero -> accu
+      accu + (Z.popcount @@ Z.of_int64 x)
+    else
+      accu + 1 + (Z.popcount @@ Z.of_int64 (pop_sign x))
       ) 0 b
 
 

ocaml/Bitlist.mli

@@ -1,4 +1,4 @@
-type t = Bit.t list
+type t
 
 (** The zero bit list *)
 val zero : Qptypes.N_int_number.t -> t

ocaml/Determinant.ml

@@ -25,19 +25,6 @@ let to_bitlist_couple x =
   in (xa,xb)
 
 
-let bitlist_to_string ~mo_num x =
-  let len =
-    MO_number.to_int mo_num
-  in
-  let s =
-    List.map (function
-        | Bit.Zero -> "-"
-        | Bit.One  -> "+"
-    ) x
-    |> String.concat ""
-  in
-  String.sub s 0 len
-
 
 
 let of_int64_array ~n_int ~alpha ~beta x =
@@ -48,37 +35,29 @@ let of_int64_array ~n_int ~alpha ~beta x =
    in
    if ( (Bitlist.popcnt a) <> alpha) then
      begin
-       let mo_num = MO_number.get_max () in
-       let mo_num = MO_number.of_int mo_num ~max:mo_num  in
        failwith (Printf.sprintf "Expected %d electrons in alpha determinant
-%s" alpha (bitlist_to_string ~mo_num:mo_num a) )
+%s" alpha (Bitlist.to_string a) )
      end;
    if ( (Bitlist.popcnt b) <> beta ) then
      begin
-       let mo_num = MO_number.get_max () in
-       let mo_num = MO_number.of_int mo_num ~max:mo_num  in
        failwith (Printf.sprintf "Expected %d electrons in beta determinant
-%s" beta (bitlist_to_string ~mo_num:mo_num b) )
+%s" beta (Bitlist.to_string b) )
      end;
    x
 
 
-let of_bitlist_couple ?n_int ~alpha ~beta (xa,xb) =
+let of_bitlist_couple ~n_int ~alpha ~beta (xa,xb) =
   let ba, bb =
     Bitlist.to_int64_array xa ,
     Bitlist.to_int64_array xb
-  and n_int =
-    match n_int with
-    | Some x -> x
-    | None -> Bitlist.n_int_of_mo_num (List.length xa)
   in
   of_int64_array ~n_int ~alpha ~beta (Array.concat [ba;bb])
 
 
 let to_string ~mo_num x =
   let (xa,xb) = to_bitlist_couple x in
-  [ "  " ; bitlist_to_string ~mo_num xa ; "\n" ;
+  [ "  " ; Bitlist.to_string xa ; "\n" ;
-    "  " ; bitlist_to_string ~mo_num xb ]
+    "  " ; Bitlist.to_string xb ]
   |> String.concat ""
 
 

ocaml/Determinant.mli

@@ -24,7 +24,7 @@ val to_alpha_beta : t -> (int64 array)*(int64 array)
 val to_bitlist_couple : t -> Bitlist.t * Bitlist.t
 
 (** Create from a bit list *)
-val of_bitlist_couple : ?n_int:Qptypes.N_int_number.t ->
+val of_bitlist_couple : n_int:Qptypes.N_int_number.t ->
   alpha:Qptypes.Elec_alpha_number.t ->
   beta:Qptypes.Elec_beta_number.t ->
   Bitlist.t * Bitlist.t -> t

ocaml/Input_determinants_by_hand.ml

@@ -472,6 +472,7 @@ psi_det                = %s
 
     (* Handle determinants *)
     let psi_det =
+      let n_int = N_int_number.of_int @@ (MO_number.get_max () - 1) / 64 + 1 in
       let n_alpha = Ezfio.get_electrons_elec_alpha_num ()
         |> Elec_alpha_number.of_int
       and n_beta = Ezfio.get_electrons_elec_beta_num ()
@@ -484,7 +485,7 @@ psi_det                = %s
             let newdet =
                (Bitlist.of_string ~zero:'-' ~one:'+' alpha ,
                 Bitlist.of_string ~zero:'-' ~one:'+' beta)
-               |> Determinant.of_bitlist_couple  ~alpha:n_alpha ~beta:n_beta
+               |> Determinant.of_bitlist_couple ~n_int ~alpha:n_alpha ~beta:n_beta
                |> Determinant.sexp_of_t
                |> Sexplib.Sexp.to_string
             in
@@ -492,9 +493,11 @@ psi_det                = %s
           end
       | _::tail -> read_dets accu tail
       in
+      (*
       let dets =
         List.map String_ext.rev dets
       in
+      *)
       let a =
         read_dets [] dets
         |> String.concat ""

ocaml/_tags

@@ -1,4 +1,4 @@
-true: package(cryptokit,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
+true: package(cryptokit,zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt)
 true: thread
 false: profile
 <*byte>   : linkdep(c_bindings.o), custom

ocaml/qp_set_mo_class.ml

@@ -112,9 +112,9 @@ let set ~core ~inact ~act ~virt ~del =
   and av = Excitation.create_single act virt
   in
   let single_excitations = [ ia ; aa ; av ]
-                           |> List.map (fun x ->
+                           |> List.map (fun z ->
                                let open Excitation in
-                               match x with
+                               match z with
                                | Single (x,y) ->
                                  ( MO_class.to_bitlist n_int (Hole.to_mo_class x),
                                    MO_class.to_bitlist n_int (Particle.to_mo_class y) )
@@ -187,9 +187,10 @@ let set ~core ~inact ~act ~virt ~del =
     match aa with
     | Double _ -> assert false
     | Single (x,y) ->
-      ( MO_class.to_bitlist n_int (Hole.to_mo_class x) ) @
+      Bitlist.to_int64_list
+        ( MO_class.to_bitlist n_int (    Hole.to_mo_class x) )  @
+      Bitlist.to_int64_list
         ( MO_class.to_bitlist n_int (Particle.to_mo_class y) )
-      |> Bitlist.to_int64_list
   in
   Ezfio.set_bitmasks_n_mask_cas 1;
   Ezfio.ezfio_array_of_list ~rank:3 ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result

src/cipsi/run_selection_slave.irp.f

@@ -85,7 +85,7 @@ subroutine run_selection_slave(thread,iproc,energy)
       if(ctask > 0) then
         call sort_selection_buffer(buf)
 !        call merge_selection_buffers(buf,buf2)
-print *, task_id(1), pt2(1), buf%cur, ctask
+!print *, task_id(1), pt2(1), buf%cur, ctask
         call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask)
 !        buf%mini = buf2%mini
         pt2(:) = 0d0

src/fci/40.fci.bats

@@ -36,24 +36,24 @@ function run_stoch() {
   [[ -n $TRAVIS ]] && skip
   qp set_file f2.ezfio
   qp set_frozen_core
-  run_stoch -199.30496 1.e-4
+  run_stoch -199.30486   1.e-4
 }
 
 @test "NH3" { # 10.6657s
   qp set_file nh3.ezfio
   qp set_mo_class --core="[1-4]" --act="[5-72]"
-  run -56.244753429144986  1.e-5
+  run -56.244753429144986  1.e-4
 }
 
 @test "DHNO" { # 11.4721s
   qp set_file dhno.ezfio
   qp set_mo_class --core="[1-7]" --act="[8-64]" 
-  run -130.459020029816 1.e-5
+  run -130.459020029816 1.e-4
 }
 
 @test "HCO" { # 12.2868s
   qp set_file hco.ezfio
-  run -113.297494345682 2.e-05 
+  run -113.297494345682 1.e-4
 }
 
 @test "H2O2" { # 12.9214s
@@ -65,82 +65,82 @@ function run_stoch() {
 @test "HBO" { # 13.3144s
   [[ -n $TRAVIS ]] && skip
   qp set_file hbo.ezfio
-  run -100.214185815312 1.e-5
+  run -100.212829869715 1.e-4
 }
 
 @test "H2O" { # 11.3727s
   [[ -n $TRAVIS ]] && skip
   qp set_file h2o.ezfio
-  run -76.2359268957699 2.e-5
+  run -76.2359268957699 1.e-4
 }
 
 @test "ClO" { # 13.3755s
   [[ -n $TRAVIS ]] && skip
   qp set_file clo.ezfio
-  run -534.546005867797 5.e-5
+  run -534.545881614967 1.e-4
 }
 
 @test "SO" { # 13.4952s
   [[ -n $TRAVIS ]] && skip
   qp set_file so.ezfio
-  run -26.0124797722154 1.e-5
+  run -26.0158153138924 1.e-4
 }
 
 @test "H2S" { # 13.6745s
   [[ -n $TRAVIS ]] && skip
   qp set_file h2s.ezfio
-  run -398.859480581924 1.e-5
+  run -398.859168655255 1.e-4
 }
 
 @test "OH" { # 13.865s 
   [[ -n $TRAVIS ]] && skip
   qp set_file oh.ezfio
-  run -75.6119887538831 1.e-05 
+  run -75.6120779012574 1.e-4 
 }
 
 @test "SiH2_3B1" { # 13.938ss
   [[ -n $TRAVIS ]] && skip
   qp set_file sih2_3b1.ezfio
-  run -290.017539006762 1.e-5
+  run -290.017539006762 1.e-4
 }
 
 @test "H3COH" { # 14.7299s
   [[ -n $TRAVIS ]] && skip
   qp set_file h3coh.ezfio
-  run -115.205054063687 1.e-5
+  run -115.205941463667 1.e-4
 }
 
 @test "SiH3" { # 15.99s
   [[ -n $TRAVIS ]] && skip
   qp set_file sih3.ezfio
-  run -5.57269434557089 2.e-05 
+  run -5.57241217753818 1.e-4 
 }
 
 @test "CH4" { # 16.1612s
   [[ -n $TRAVIS ]] && skip
   qp set_file ch4.ezfio
   qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
-  run -40.2409059687324 2.e-5
+  run -40.2409678239136 1.e-4
 }
 
 @test "ClF" { # 16.8864s
   [[ -n $TRAVIS ]] && skip
   qp set_file clf.ezfio
-  run -559.170406471496 1.e-5
+  run -559.170272077166 1.e-4
 }
 
 @test "SO2" { # 17.5645s
   [[ -n $TRAVIS ]] && skip
   qp set_file so2.ezfio
   qp set_mo_class --core="[1-8]" --act="[9-87]" 
-  run -41.5746738713298 5.e-5
+  run -41.5746738713298 1.e-4
 }
 
 @test "C2H2" { # 17.6827s
   [[ -n $TRAVIS ]] && skip
   qp set_file c2h2.ezfio
   qp set_mo_class --act="[1-30]" --del="[31-36]"
-  run -12.3670840202635 2.e-5
+  run -12.3656179738175 1.e-4
 }
 
 @test "N2" { # 18.0198s
@@ -154,14 +154,14 @@ function run_stoch() {
   [[ -n $TRAVIS ]] && skip
   qp set_file n2h4.ezfio
   qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]"
-  run -111.367234092521 2.e-5
+  run -111.367332681559 1.e-4
 }
 
 @test "CO2" { # 21.1748s
   [[ -n $TRAVIS ]] && skip
   qp set_file co2.ezfio
   qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]"
-  run -187.969676381867 1.e-5
+  run -187.968599504402 1.e-4
 }
 
 
@@ -169,13 +169,13 @@ function run_stoch() {
   [[ -n $TRAVIS ]] && skip
   qp set_file cu_nh3_4_2plus.ezfio
   qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
-  run -1862.98610987882  1.e-05
+  run -1862.98614665139  1.e-04
 }
 
 @test "HCN" { # 20.3273s
   [[ -n $TRAVIS ]] && skip
   qp set_file hcn.ezfio
   qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]"
-  run -93.0799328685679 2.e-5
+  run -93.0728641601823 1.e-4
 }