LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 12:23:31 +01:00
Code Releases Activity

FCI with direct davidson

Browse Source
This commit is contained in:
Anthony Scemama 2019-04-03 18:09:13 +02:00
parent f00a490b5e
commit 4de337619f
4 changed files with 291 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

128
CI/CI.ml
View File

@ -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 ()

14
Utils/Davidson.ml
View File

@ -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 =

205
Utils/Matrix.ml
View File

@ -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 () =

10
Utils/Matrix.mli
View File

@ -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. *)