mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-08-31 00:33:40 +02:00
Davidson with sparse matrix
This commit is contained in:
parent
c5a16d0765
commit
a976bacfdf
40
CI/CI.ml
40
CI/CI.ml
@ -5,8 +5,8 @@ module Ds = Determinant_space
|
|||||||
type t =
|
type t =
|
||||||
{
|
{
|
||||||
det_space : Ds.t ;
|
det_space : Ds.t ;
|
||||||
m_H : Mat.t lazy_t ;
|
m_H : Matrix.t lazy_t ;
|
||||||
m_S2 : Mat.t lazy_t ;
|
m_S2 : Matrix.t lazy_t ;
|
||||||
eigensystem : (Mat.t * Vec.t) lazy_t;
|
eigensystem : (Mat.t * Vec.t) lazy_t;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -53,18 +53,19 @@ let make det_space =
|
|||||||
let ndet = Ds.size det_space in
|
let ndet = Ds.size det_space in
|
||||||
let det = Ds.determinants det_space in
|
let det = Ds.determinants det_space in
|
||||||
let mo_basis = Ds.mo_basis det_space in
|
let mo_basis = Ds.mo_basis det_space in
|
||||||
(*
|
|
||||||
let m_H = lazy (
|
let m_H = lazy (
|
||||||
Util.list_range 0 (ndet-1)
|
Array.init ndet (fun i -> let ki = det.(i) in
|
||||||
|> List.map (fun i -> let ki = det.(i) in
|
Vec.init ndet (fun j -> let kj = det.(j-1) in
|
||||||
Array.init ndet (fun j -> let kj = det.(j) in
|
|
||||||
h_ij mo_basis ki kj
|
h_ij mo_basis ki kj
|
||||||
)
|
)
|
||||||
|> Vec.of_array)
|
|> Vector.sparse_of_vec
|
||||||
|> Mat.of_col_vecs_list
|
)
|
||||||
|
|> Matrix.sparse_of_vector_array
|
||||||
)
|
)
|
||||||
in
|
in
|
||||||
*)
|
|
||||||
|
|
||||||
(*
|
(*
|
||||||
let m_H = lazy (
|
let m_H = lazy (
|
||||||
let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in
|
let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in
|
||||||
@ -93,6 +94,8 @@ let make det_space =
|
|||||||
|> Mat.of_col_vecs_list
|
|> Mat.of_col_vecs_list
|
||||||
) in
|
) in
|
||||||
*)
|
*)
|
||||||
|
|
||||||
|
(* Parallel
|
||||||
let m_H = lazy (
|
let m_H = lazy (
|
||||||
let h =
|
let h =
|
||||||
if Parallel.master then
|
if Parallel.master then
|
||||||
@ -120,12 +123,15 @@ let make det_space =
|
|||||||
List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l);
|
List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l);
|
||||||
Parallel.broadcast (lazy h)
|
Parallel.broadcast (lazy h)
|
||||||
) in
|
) in
|
||||||
|
*)
|
||||||
let m_S2 = lazy (
|
let m_S2 = lazy (
|
||||||
Array.init ndet (fun i -> let ki = det.(i) in
|
Array.init ndet (fun i -> let ki = det.(i) in
|
||||||
Array.init ndet (fun j -> let kj = det.(j) in
|
Vec.init ndet (fun j -> let kj = det.(j-1) in
|
||||||
CIMatrixElement.make_s2 ki kj
|
CIMatrixElement.make_s2 ki kj
|
||||||
))
|
)
|
||||||
|> Mat.of_array
|
|> Vector.sparse_of_vec
|
||||||
|
)
|
||||||
|
|> Matrix.sparse_of_vector_array
|
||||||
)
|
)
|
||||||
in
|
in
|
||||||
let eigensystem = lazy (
|
let eigensystem = lazy (
|
||||||
@ -134,9 +140,13 @@ let make det_space =
|
|||||||
Parallel.broadcast @@
|
Parallel.broadcast @@
|
||||||
lazy (Util.diagonalize_symm m_H)
|
lazy (Util.diagonalize_symm m_H)
|
||||||
*)
|
*)
|
||||||
let diagonal = Vec.init (Mat.dim1 m_H) (fun i -> m_H.{i,i}) in
|
let diagonal =
|
||||||
let matrix_vector psi = symv m_H psi in
|
Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
|
||||||
Davidson.make diagonal matrix_vector
|
in
|
||||||
|
let matrix_prod psi =
|
||||||
|
Matrix.mm ~transa:`T m_H psi
|
||||||
|
in
|
||||||
|
Davidson.make diagonal matrix_prod
|
||||||
)
|
)
|
||||||
in
|
in
|
||||||
{ det_space ; m_H ; m_S2 ; eigensystem }
|
{ det_space ; m_H ; m_S2 ; eigensystem }
|
||||||
|
@ -8,7 +8,7 @@ let make
|
|||||||
?(n_iter=10)
|
?(n_iter=10)
|
||||||
?(threshold=1.e-8)
|
?(threshold=1.e-8)
|
||||||
diagonal
|
diagonal
|
||||||
matrix_vector
|
matrix_prod
|
||||||
=
|
=
|
||||||
|
|
||||||
let n = Vec.dim diagonal in (* Size of the matrix to diagonalize *)
|
let n = Vec.dim diagonal in (* Size of the matrix to diagonalize *)
|
||||||
@ -70,7 +70,12 @@ let make
|
|||||||
|
|
||||||
(* Apply the operator the m last vectors *)
|
(* Apply the operator the m last vectors *)
|
||||||
let w_new =
|
let w_new =
|
||||||
List.map matrix_vector u_new_ortho
|
matrix_prod (
|
||||||
|
u_new_ortho
|
||||||
|
|> Mat.of_col_vecs_list
|
||||||
|
|> Matrix.dense_of_mat )
|
||||||
|
|> Matrix.to_mat
|
||||||
|
|> Mat.to_col_vecs_list
|
||||||
in
|
in
|
||||||
|
|
||||||
(* Data for the next iteration *)
|
(* Data for the next iteration *)
|
||||||
|
@ -60,6 +60,10 @@ let dense_of_sparse = function
|
|||||||
|
|
||||||
let dense_of_mat m = Dense m
|
let dense_of_mat m = Dense m
|
||||||
|
|
||||||
|
let rec to_vector_array ?(threshold=epsilon) = function
|
||||||
|
| Sparse {m ; n ; v} -> v
|
||||||
|
| Dense m -> to_vector_array (sparse_of_dense ~threshold (Dense m))
|
||||||
|
|
||||||
|
|
||||||
let sparse_of_mat ?(threshold=epsilon) m =
|
let sparse_of_mat ?(threshold=epsilon) m =
|
||||||
dense_of_mat m
|
dense_of_mat m
|
||||||
|
@ -27,7 +27,7 @@ val get : t -> int -> int -> float
|
|||||||
val to_mat : t -> Mat.t
|
val to_mat : t -> Mat.t
|
||||||
(** Convert into a Lacaml Mat. *)
|
(** Convert into a Lacaml Mat. *)
|
||||||
|
|
||||||
val to_vector_array : t -> Vector.t array
|
val to_vector_array : ?threshold:float -> t -> Vector.t array
|
||||||
(** Convert the matrix into an array of column vectors. *)
|
(** Convert the matrix into an array of column vectors. *)
|
||||||
|
|
||||||
val sparse_of_dense : ?threshold:float -> t -> t
|
val sparse_of_dense : ?threshold:float -> t -> t
|
||||||
|
Loading…
Reference in New Issue
Block a user