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

Added CI

Browse Source
This commit is contained in:
Anthony Scemama 2024-07-30 12:55:42 +02:00
parent a9c5fd16a9
commit 5f47e96962
10 changed files with 1859 additions and 115 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1183
ci/lib/ci.ml Normal file
View File

File diff suppressed because it is too large Load Diff

283
ci/lib/ci_matrix_element.ml Normal file
View File

@ -0,0 +1,283 @@
open Common
module De = Determinant
module Ex = Excitation
module Sp = Spindeterminant
type t = float list
(** Computes non-zero integral values.
@param integrals A list of tuples containing one-electron, two-electron integrals, and optionally three-electron integrals.
@param degree_a The degree of excitation in alpha spin orbitals.
@param degree_b The degree of excitation in beta spin orbitals.
@param ki The initial Slater determinant.
@param kj The final Slater determinant.
@return TODO A list of computed integral values based on the specified degrees and molecular orbitals.
This function performs the following operations:
- Converts spin determinants to lists of molecular orbitals for alpha and beta spins.
- Defines helper functions for singly and doubly excited determinants, and diagonal elements
- Uses lazy evaluation to defer computations until required.
- Supports both two-electron and three-electron integrals, handling phase factors appropriately.
Example usage:
{[
let integrals = [(one_e_integral, two_e_integral, None); (one_e_integral, two_e_integral, Some three_e_integral)] in
let result = non_zero integrals 1 1 ki kj
]}
*)
let non_zero integrals degree_a degree_b ki kj =
let kia = De.alfa ki and kib = De.beta ki
and kja = De.alfa kj and kjb = De.beta kj
in
let single h p spin same opposite =
let same_spin_mo_list =
Sp.to_list same
and opposite_spin_mo_list =
Sp.to_list opposite
in
fun one_e two_e ->
let same_spin =
List.fold_left (fun accu i -> accu +. two_e h i p i spin spin) 0. same_spin_mo_list
and opposite_spin =
List.fold_left (fun accu i -> accu +. two_e h i p i spin (Spin.other spin) ) 0. opposite_spin_mo_list
in (one_e h p spin) +. same_spin +. opposite_spin
in
let diag_element =
let mo_a = Sp.to_list kia
and mo_b = Sp.to_list kib
in
fun one_e two_e ->
let one =
(List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a)
+.
(List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b)
in
let two =
let rec aux_same spin accu = function
| [] -> accu
| i :: rest ->
let new_accu =
List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest
in
(aux_same [@tailcall]) spin new_accu rest
in
let rec aux_opposite accu other = function
| [] -> accu
| i :: rest ->
let new_accu =
List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other
in
(aux_opposite [@tailcall]) new_accu other rest
in
(aux_same Spin.Alfa 0. mo_a) +. (aux_same Spin.Beta 0. mo_b) +.
(aux_opposite 0. mo_a mo_b)
in
one +. two
in
let result_2e = lazy (
match degree_a, degree_b with
| 1, 1 -> (* alpha-beta double *)
begin
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
match phase_a, phase_b with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
end
| 2, 0 -> (* alpha double *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
match phase with
| Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
| Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
end
| 0, 2 -> (* beta double *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in
match phase with
| Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
| Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
end
| 1, 0 -> (* alpha single *)
begin
let h, p, phase = Ex.single_of_spindet kia kja in
match phase with
| Phase.Pos -> fun one_e two_e -> single h p Spin.Alfa kia kib one_e two_e
| Phase.Neg -> fun one_e two_e -> -. single h p Spin.Alfa kia kib one_e two_e
end
| 0, 1 -> (* beta single *)
begin
let h, p, phase = Ex.single_of_spindet kib kjb in
match phase with
| Phase.Pos -> fun one_e two_e -> single h p Spin.Beta kib kia one_e two_e
| Phase.Neg -> fun one_e two_e -> -. single h p Spin.Beta kib kia one_e two_e
end
| 0, 0 -> (* diagonal element *)
diag_element
| _ -> assert false
) in
let result_3e = lazy (
match degree_a, degree_b with
| 1, 1 -> (* alpha-beta double *)
begin
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
match phase_a, phase_b with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg -> fun _ two_e _ -> two_e ha hb pa pb Spin.Alfa Spin.Beta
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg -> fun _ two_e _ -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
end
| 2, 0 -> (* alpha double *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
match phase with
| Phase.Pos -> fun _ two_e _ -> two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
| Phase.Neg -> fun _ two_e _ -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
end
| 0, 2 -> (* beta double *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in
match phase with
| Phase.Pos -> fun _ two_e _ -> two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
| Phase.Neg -> fun _ two_e _ -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
end
| 1, 0 -> (* alpha single *)
begin
let h, p, phase = Ex.single_of_spindet kia kja in
match phase with
| Phase.Pos -> fun one_e two_e _ -> single h p Spin.Alfa kia kib one_e two_e
| Phase.Neg -> fun one_e two_e _ -> -. single h p Spin.Alfa kia kib one_e two_e
end
| 0, 1 -> (* beta single *)
begin
let h, p, phase = Ex.single_of_spindet kib kjb in
match phase with
| Phase.Pos -> fun one_e two_e _ -> single h p Spin.Beta kib kia one_e two_e
| Phase.Neg -> fun one_e two_e _ -> -. single h p Spin.Beta kib kia one_e two_e
end
| 0, 0 -> (* diagonal element *)
fun one_e two_e _ -> diag_element one_e two_e
| 3, 0 -> (* alpha triple *)
begin
let h1, p1, h2, p2, h3, p3, phase = Ex.triple_of_spindet kia kja in
match phase with
| Phase.Pos -> fun _ _ three_e -> three_e h1 h2 h3 p1 p2 p3 Spin.Alfa Spin.Alfa Spin.Alfa
| Phase.Neg -> fun _ _ three_e -> -. three_e h1 h2 h3 p1 p2 p3 Spin.Alfa Spin.Alfa Spin.Alfa
end
| 0, 3 -> (* beta triple *)
begin
let h1, p1, h2, p2, h3, p3, phase = Ex.triple_of_spindet kib kja in
match phase with
| Phase.Pos -> fun _ _ three_e -> three_e h1 h2 h3 p1 p2 p3 Spin.Beta Spin.Beta Spin.Beta
| Phase.Neg -> fun _ _ three_e -> -. three_e h1 h2 h3 p1 p2 p3 Spin.Beta Spin.Beta Spin.Beta
end
| 2, 1 -> (* alpha2 beta triple *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
let h3, p3, phase' = Ex.single_of_spindet kib kjb in
match phase, phase' with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg ->
fun _ _ three_e -> three_e h1 h2 h3 p1 p2 p3 Spin.Alfa Spin.Alfa Spin.Beta
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg ->
fun _ _ three_e -> -. three_e h1 h2 h3 p1 p2 p3 Spin.Alfa Spin.Alfa Spin.Beta
end
| 1, 2 -> (* alpha beta2 triple *)
begin
let h1, p1, phase = Ex.single_of_spindet kia kja in
let h2, p2, h3, p3, phase' = Ex.double_of_spindet kib kjb in
match phase, phase' with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg ->
fun _ _ three_e -> three_e h1 h2 h3 p1 p2 p3 Spin.Alfa Spin.Beta Spin.Beta
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg ->
fun _ _ three_e -> -. three_e h1 h2 h3 p1 p2 p3 Spin.Alfa Spin.Beta Spin.Beta
end
| _ -> fun _ _ _ -> 0.
) in
List.map (fun (one_e, two_e, x) ->
match x with
| None -> (Lazy.force result_2e) one_e two_e
| Some three_e -> (Lazy.force result_3e) one_e two_e three_e
) integrals
let make integrals ki kj =
let degree_a, degree_b =
De.excitation_levels ki kj
in
if degree_a+degree_b > 2 then
List.map (fun _ -> 0.) integrals
else
non_zero integrals degree_a degree_b ki kj
let make_s2 ki kj =
let degree_a = De.excitation_level_alfa ki kj in
let kia = De.alfa ki in
let kja = De.alfa kj in
if degree_a > 1 then 0.
else
let degree_b = De.excitation_level_beta ki kj in
let kib = De.beta ki in
let kjb = De.beta kj in
match degree_a, degree_b with
| 1, 1 -> (* alpha-beta double *)
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
if ha = pb && hb = pa then
begin
match phase_a, phase_b with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg -> -1.
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg -> 1.
end
else 0.
| 0, 0 ->
let ba = Sp.bitstring kia and bb = Sp.bitstring kib in
let tmp = Bitstring.logxor ba bb in
let n_a = Bitstring.logand ba tmp |> Bitstring.popcount in
let n_b = Bitstring.logand bb tmp |> Bitstring.popcount in
let s_z = 0.5 *. float_of_int (n_a - n_b) in
float_of_int n_a +. s_z *. (s_z -. 1.)
| _ -> 0.

67
ci/lib/ci_matrix_element.mli Normal file
View File

@ -0,0 +1,67 @@
open Common
module De = Determinant
type t
val make :
((int -> int -> Spin.t -> float) *
(int -> int -> int -> int -> Spin.t -> Spin.t -> float) *
(int -> int -> int -> int -> int -> int -> Spin.t -> Spin.t -> Spin.t -> float) option) list ->
De.t -> De.t -> float list
(** [make integrals ki kj] Computes matrix elements for multiple operators between two Slater
determinants, or returns zeros if the total degree of excitation exceeds 2.
@param integrals A list of tuples containing one-electron, two-electron
integrals, and optionally three-electron integrals, for each operator
@param ki The initial Slater determinant.
@param kj The final Slater determinant.
@return A list of computed matrix elements or zeroes if the total excitation
lebel (degree_a + degree_b) is greater than 2.
Example usage:
{[
let integrals = [(one_e_integral, two_e_integral, None); (one_e_integral', two_e_integral', Some three_e_integral')] in
let result = make integrals ki kj
]}
*)
val make_s2 : De.t -> De.t -> float
(** [make_s2 ki kj] computes the value of the $S^2$ operator for two determinants.
@param ki The initial spin determinant.
@param kj The final spin determinant.
@return The computed value of the $S^2$ operator.
Example usage:
{[
let s2_value = make_s2 ki kj
]}
*)
(** Computes matrix elements when the user knows they are non-zero.
@param integrals A list of tuples containing one-electron, two-electron integrals, and optionally three-electron integrals.
@param degree_a The degree of excitation in alpha spin orbitals.
@param degree_b The degree of excitation in beta spin orbitals.
@param ki The initial Slater determinant.
@param kj The final Slater determinant.
@return A list of matrix elements for multiple operators
Example usage:
{[
let integrals = [(one_e_integral, two_e_integral, None); (one_e_integral, two_e_integral, Some three_e_integral)] in
let result = non_zero integrals 1 1 ki kj
]}
*)
val non_zero :
((int -> int -> Spin.t -> float) *
(int -> int -> int -> int -> Spin.t -> Spin.t -> float) *
(int -> int -> int -> int -> int -> int -> Spin.t -> Spin.t -> Spin.t -> float) option) list ->
int -> int -> De.t -> De.t -> float list

162
linear_algebra/lib/davidson.ml Normal file
View File

@ -0,0 +1,162 @@
open Common
let (%.) = Vector.(%.)
let (%:) = Matrix.(%:)
type t
let make
?guess
?(n_states=1)
?(n_iter=8)
?(threshold=1.e-7)
diagonal
matrix_prod
=
(* Size of the matrix to diagonalize *)
let n = Vector.dim diagonal in
let m = n_states in
(* Create guess vectors u, with unknown vectors initialized to unity. *)
let init_vectors m =
let result = Matrix.make n m 0. in
for i=0 to m-1 do
Matrix.set result i i 1.;
done;
result
in
let guess =
match guess with
| Some vectors -> vectors
| None -> init_vectors m
in
let guess =
if Matrix.dim2 guess = n_states then guess
else
(Matrix.to_col_vecs_list guess) @
(Matrix.to_col_vecs_list (init_vectors (m-(Matrix.dim2 guess))) )
|> Matrix.of_col_vecs_list
in
let pick_new u =
Matrix.to_col_vecs_list u
|> Util.list_pack m
|> List.rev
|> List.hd
in
let u_new = Matrix.to_col_vecs_list guess in
let rec iteration u u_new w iter macro =
(* u is a list of orthonormal vectors, on which the operator has
been applied : w = op.u
u_new is a list of vectors which will increase the size of the
space.
*)
(* Orthonormalize input vectors u_new *)
let u_new_ortho =
List.concat [u ; u_new]
|> Matrix.of_col_vecs_list
|> Orthonormalization.qr_ortho
|> pick_new
in
(* Apply the operator to the m last vectors *)
let w_new =
matrix_prod (
u_new_ortho
|> Matrix.of_col_vecs_list)
|> Matrix.to_col_vecs_list
in
(* Data for the next iteration *)
let u_next =
List.concat [ u ; u_new_ortho ]
and w_next =
List.concat [ w ; w_new ]
in
(* Build the small matrix h = <U_k | W_l> *)
let m_U =
Matrix.of_col_vecs_list u_next
and m_W =
Matrix.of_col_vecs_list w_next
in
let m_h =
Matrix.gemm_tn m_U m_W
in
(* Diagonalize h *)
let y, lambda =
Matrix.diagonalize_symm m_h
in
(* Express m lowest eigenvectors of h in the large basis *)
let m_new_U =
Matrix.gemm ~n:m m_U y
and m_new_W =
Matrix.gemm ~n:m m_W y
in
(* Compute the residual as proposed new vectors *)
let u_proposed =
Matrix.init_cols n m (fun i k ->
let delta = lambda%.(k) -. diagonal%.(i) in
let delta =
if abs_float delta > 1.e-2 then delta
else if delta > 0. then 1.e-2
else (-1.e-2)
in
(lambda%.(k) *. (m_new_U%:(i,k)) -. (m_new_W%:(i,k)) ) /. delta
)
|> Matrix.to_col_vecs_list
in
let residual_norms = List.rev @@ List.rev_map Vector.norm u_proposed in
let residual_norm =
List.fold_left (fun accu i -> accu +. i *. i) 0. residual_norms
|> sqrt
in
Printf.printf "%3d " iter;
Vector.iteri (fun i x -> if (i<=m) then Printf.printf "%16.10f " x) lambda;
Printf.printf "%16.8e%!\n" residual_norm;
(* Make new vectors sparse *)
let u_proposed =
Matrix.of_col_vecs_list u_proposed
in
let maxu = Matrix.norm ~l:`Linf u_proposed in
let thr = maxu *. 0.00001 in
let u_proposed =
Matrix.map (fun x -> if abs_float x < thr then 0. else x) u_proposed
|> Matrix.to_col_vecs_list
in
(*
Format.printf "%a@." Matrix.pp_matrix @@ m_new_U;
*)
if residual_norm > threshold && macro > 0 then
let u_next, w_next, iter, macro =
if iter = n_iter then
m_new_U |> pick_new,
m_new_W |> pick_new,
0, (macro-1)
else
u_next, w_next, (iter+1), macro
in
(iteration [@tailcall]) u_next u_proposed w_next iter macro
else
(m_new_U |> pick_new |> Matrix.of_col_vecs_list), lambda
in
iteration [] u_new [] 1 30

21
linear_algebra/lib/davidson.mli Normal file
View File

@ -0,0 +1,21 @@
type t
val make :
?guess:('a, 'b) Matrix.t ->
?n_states:int ->
?n_iter:int ->
?threshold:float ->
'a Vector.t ->
(('a, 'b) Matrix.t -> ('a, 'b) Matrix.t) ->
('a, 'b) Matrix.t * 'b Vector.t
(** Performs a Davidson diagonalization. Example:
let eigenvectors, eigenvalues =
Davidson.make diagonal mat_prod
- [diagonal] contains the diagonal of the matrix to diagonalize
- [mat_prod] is a function performing a matrix multiplication of the matrix to
diagonalize with the matrix containing the current set of vectors
*)

1
linear_algebra/lib/linear_algebra.ml
View File

@ -8,3 +8,4 @@ module Matrix = Matrix
module Orthonormalization = Orthonormalization
module Spherical_to_cartesian = Spherical_to_cartesian
module Vector = Vector
module Davidson = Davidson

97
linear_algebra/lib/matrix.ml
View File

@ -12,7 +12,7 @@ let out_of_place f t =
result
let relabel t = t
let make m n x = Mat.make m n x
let make0 m n = Mat.make0 m n
let create m n = Mat.create m n
@ -53,16 +53,16 @@ let of_bigarray t = lacpy t
let reshape_inplace m n a =
assert ( (dim1 a) * (dim2 a) = m*n);
let b =
to_bigarray a
let b =
to_bigarray a
|> Bigarray.genarray_of_array2
in
Bigarray.reshape_2 b m n
let reshape_vec_inplace m n v =
assert ( Vector.dim v = m*n);
let b =
Vector.to_bigarray_inplace v
let b =
Vector.to_bigarray_inplace v
|> Bigarray.genarray_of_array1
in
Bigarray.reshape_2 b m n
@ -96,7 +96,7 @@ let to_col_vecs_list t =
|> List.rev_map Vector.of_bigarray_inplace
|> List.rev
let detri_inplace t =
let detri_inplace t =
Mat.detri t
let detri t =
@ -106,7 +106,7 @@ let as_vec_inplace t =
Mat.as_vec t
|> Vector.of_bigarray_inplace
let as_vec t =
let as_vec t =
lacpy t
|> Mat.as_vec
|> Vector.of_bigarray_inplace
@ -114,26 +114,26 @@ let as_vec t =
let random ?rnd_state ?(from= -. 1.0) ?(range=2.0) m n =
Mat.random ?rnd_state ~from ~range m n
let map_inplace f ~b a =
let map_inplace f ~b a =
ignore @@ Mat.map f ~b a
let map f a =
let map f a =
Mat.map f a
let scale_inplace x t =
Mat.scal x t
let scale_inplace x t =
Mat.scal x t
let scale x t =
let scale x t =
out_of_place (fun t -> scale_inplace x t) t
let of_diag v =
Vector.to_bigarray_inplace v
|> Mat.of_diag
|> Mat.of_diag
let diag t =
Mat.copy_diag t
|> Vector.of_bigarray_inplace
let gemv_n_inplace ?m ?n ?(beta=0.) y ?(alpha=1.) ?(ar=1) ?(ac=1) t v =
let y = Vector.to_bigarray_inplace y in
let v = Vector.to_bigarray_inplace v in
@ -203,26 +203,26 @@ let gemm_tt ?m ?n ?k ?beta ?c ?(alpha=1.0) a b =
let gemm_trace ?(transa=`N) a ?(transb=`N) b =
Mat.gemm_trace ~transa a ~transb b
let gemm_nn_trace a b = gemm_trace a b ~transa:`N ~transb:`N
let gemm_nt_trace a b = gemm_trace a b ~transa:`N ~transb:`T
let gemm_tn_trace a b = gemm_trace a b ~transa:`T ~transb:`N
let gemm_tt_trace a b = gemm_trace a b ~transa:`T ~transb:`T
let gemm_nn_trace a b = gemm_trace a b ~transa:`N ~transb:`N
let gemm_nt_trace a b = gemm_trace a b ~transa:`N ~transb:`T
let gemm_tn_trace a b = gemm_trace a b ~transa:`T ~transb:`N
let gemm_tt_trace a b = gemm_trace a b ~transa:`T ~transb:`T
let init_cols = Mat.init_cols
let of_col_vecs a =
Array.map Vector.to_bigarray_inplace a
|> Mat.of_col_vecs
|> Mat.of_col_vecs
let of_col_vecs_list a =
List.rev_map Vector.to_bigarray_inplace a
|> List.rev
|> Mat.of_col_vecs_list
|> Mat.of_col_vecs_list
let of_array a =
Bigarray.Array2.of_array Bigarray.Float64 Bigarray.fortran_layout a
let to_array a =
let result = Array.make_matrix (Mat.dim1 a) (Mat.dim2 a) 0. in
for i=1 to Mat.dim1 a do
@ -232,12 +232,19 @@ let to_array a =
done
done;
result
let normalize_mat_inplace t =
let norm ?(l=`L2) t =
match l with
| `L2 -> lange ~norm:`F t
| `L1 -> lange ~norm:`O t
| `Linf -> lange ~norm:`I t
let normalize_mat_inplace t =
let norm = Mat.as_vec t |> nrm2 in
Mat.scal norm t
let normalize_mat t =
let normalize_mat t =
out_of_place normalize_mat_inplace t
let diagonalize_symm m_H =
@ -261,8 +268,8 @@ let x_o_xt ~o ~x =
|> gemm x
let amax t =
Mat.as_vec t |> amax
Mat.as_vec t |> amax
let pp ppf m =
let rows = Mat.dim1 m
and cols = Mat.dim2 m
@ -285,7 +292,7 @@ let pp ppf m =
let sysv_inplace ~b a =
sysv a b
let sysv ~b a =
let sysv ~b a =
out_of_place (fun b -> sysv_inplace ~b a) b
let debug_matrix name a =
@ -339,10 +346,20 @@ let sym_of_file filename =
let copy ?m ?n ?br ?bc ?ar ?ac a =
lacpy a ?m ?n ?br ?bc ?ar ?ac
lacpy a ?m ?n ?br ?bc ?ar ?ac
let copy_inplace ?m ?n ?br ?bc ~b ?ar ?ac a =
ignore @@ lacpy ?m ?n ?br ?bc ~b ?ar ?ac a
ignore @@ lacpy ?m ?n ?br ?bc ~b ?ar ?ac a
let copy_row ?vec a i =
let result =
match vec with
| None -> Mat.copy_row a i
| Some v ->
let vec = Vector.to_bigarray_inplace v in
Mat.copy_row ~vec a i
in
Vector.of_bigarray_inplace result
let scale_cols_inplace a v =
Vector.to_bigarray_inplace v
@ -366,13 +383,13 @@ let scale_rows v a =
a'
let svd a =
let d, u, vt =
let d, u, vt =
gesvd (lacpy a)
in
u, (Vector.of_bigarray_inplace d), vt
let svd' a =
let d, u, vt =
let d, u, vt =
gesvd (lacpy a)
in
u, (Vector.of_bigarray_inplace d), vt
@ -394,7 +411,7 @@ let qr a =
let q = result in
q, r
let exponential a =
assert (dim1 a = dim2 a);
let rec loop result accu n =
@ -416,7 +433,7 @@ let exponential a =
let exponential_antisymmetric a =
let n = dim1 a in
assert (n = dim2 a);
let a2 = gemm a a in
@ -424,13 +441,13 @@ let exponential_antisymmetric a =
let tau = Vector.map (fun x -> -. sqrt x) w in
let cos_tau = Vector.cos tau in
let sin_tau_tau = Vector.mul (Vector.sin tau) (Vector.reci tau) in
let result =
let result =
add (gemm (scale_cols u cos_tau) vt) (gemm (scale_cols u sin_tau_tau) @@ gemm vt a)
in
(* Post-condition: Check if exp(-A) * exp(A) = I *)
let id = identity n in
let test =
let test =
gemm_tn ~beta:(-.1.0) ~c:id result result
|> amax
|> abs_float
@ -449,7 +466,7 @@ let to_file ~filename ?(sym=false) ?(cutoff=0.) t =
let n = Mat.dim1 t in
let m = Mat.dim2 t in
if sym then
if sym then
for j=1 to n do
for i=1 to j do
if (abs_float (t.{i,j}) > cutoff) then
@ -467,8 +484,8 @@ let to_file ~filename ?(sym=false) ?(cutoff=0.) t =
let (%:) t (i,j) = t.{i,j}
let set t i j v = t.{i,j} <- v

106
linear_algebra/lib/matrix.mli
View File

@ -3,10 +3,10 @@
type ('a,'b) t
val dim1: ('a,'b) t -> int
(** First dimension of the matrix *)
(** First dimension of the matrix *)
val dim2: ('a,'b) t -> int
(** Second dimension of the matrix *)
(** Second dimension of the matrix *)
val make : int -> int -> float -> ('a,'b) t
(** Creates a matrix initialized with the given value.
@ -36,58 +36,61 @@ val identity: int -> ('a,'b) t
val of_diag: 'a Vector.t -> ('a,'a) t
(** Creates a diagonal matrix. *)
val diag: ('a,'a) t -> 'a Vector.t
val diag: ('a,'a) t -> 'a Vector.t
(** Returns the diagonal of a matrix. *)
val fill_inplace: ('a,'b) t -> float -> unit
(** Fills the matrix with the give value. *)
(*
val add_const_diag : float -> ('a,'b) t -> ('a,'b) t
val add_const_diag : float -> ('a,'b) t -> ('a,'b) t
(** Adds a constant to the diagonal *)
val add_const_diag_inplace : float -> ('a,'b) t -> unit
(** Adds a constant to the diagonal *)
val add_const_inplace : float -> ('a,'b) t -> unit
val add_const_inplace : float -> ('a,'b) t -> unit
(** Adds a constant to the diagonal *)
*)
val add_const : float -> ('a,'b) t -> ('a,'b) t
(** Adds a constant to the diagonal *)
val add : ('a,'b) t -> ('a,'b) t -> ('a,'b) t
(** Adds two matrices *)
val sub : ('a,'b) t -> ('a,'b) t -> ('a,'b) t
(** Subtracts two matrices *)
val mul : ('a,'b) t -> ('a,'b) t -> ('a,'b) t
(** Multiplies two matrices element-wise *)
val div : ('a,'b) t -> ('a,'b) t -> ('a,'b) t
(** Divides two matrices element-wise *)
val amax : ('a,'b) t -> float
(** Maximum of the absolute values of the elements of the matrix. *)
val norm : ?l:[< `L1 | `L2 | `Linf > `L2 ] -> ('a,'b) t -> float
(** $L^1$, $L^2$ or $L^\infty$ norm of the matrix *)
val add_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit
(** [add_inplace c a b] : performs [c = a+b] in-place. *)
val sub_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit
(** [sub_inplace c a b] : performs [c = a+b] in-place. *)
val mul_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit
(** [mul_inplace c a b] : performs [c = a*b] element-wise in-place. *)
val div_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit
(** [div_inplace c a b] : performs [c = a/b] element-wise in-place. *)
(*
val to_bigarray : ('a,'b) t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t
(** Converts the matrix into a Bigarray in Fortran layout *)
val of_bigarray : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> ('a,'b) t
val of_bigarray : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> ('a,'b) t
(** Converts a [Bigarray.Array2] in Fortran layout into a matrix *)
*)
@ -96,7 +99,7 @@ val to_bigarray_inplace :
(** Converts the matrix into a Bigarray in Fortran layout in place*)
val of_bigarray_inplace :
(float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> ('a,'b) t
(float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> ('a,'b) t
(** Converts a [Bigarray.Array2] in Fortran layout into a matrix in place*)
val to_col_vecs : ('a,'b) t -> 'a Vector.t array
@ -118,10 +121,13 @@ val of_array : float array array -> ('a,'b) t
(** Converts an array of arrays into a matrix *)
val copy: ?m:int -> ?n:int -> ?br:int -> ?bc:int -> ?ar:int -> ?ac:int -> ('a,'b) t -> ('a,'b) t
(** Copies all or part of a two-dimensional matrix A to a new matrix B *)
(** Copies all or part of a matrix A to a new matrix B *)
val copy_inplace: ?m:int -> ?n:int -> ?br:int -> ?bc:int -> b:('a,'b) t -> ?ar:int -> ?ac:int -> ('a,'b) t -> unit
(** Copies all or part of a two-dimensional matrix A to an existing matrix B *)
(** Copies all or part of a matrix A to an existing matrix B *)
val copy_row : ?vec:('b Vector.t) -> ('a,'b) t -> int -> 'b Vector.t
(** Copies a given row of a matrix into a vector *)
(*
val col: ('a,'b) t -> int -> 'a Vector.t
@ -137,20 +143,20 @@ val transpose: ('a,'b) t -> ('b,'a) t
val transpose_inplace: ('a,'a) t -> unit
(** Transposes the matrix in place. *)
val detri: ('a,'b) t -> ('a,'b) t
val detri: ('a,'b) t -> ('a,'b) t
(** Takes an upper-triangular matrix, and makes it a symmetric matrix
by mirroring the defined triangle along the diagonal. *)
val detri_inplace: ('a,'b) t -> unit
val detri_inplace: ('a,'b) t -> unit
(** Takes an upper-triangular matrix, and makes it a symmetric matrix
by mirroring the defined triangle along the diagonal. *)
val as_vec_inplace: ('a,'b) t -> ('a*'b) Vector.t
(** Interpret the matrix as a vector (reshape). *)
val as_vec: ('a,'b) t -> ('a*'b) Vector.t
(** Return a copy of the reshaped matrix into a vector *)
val random: ?rnd_state:Random.State.t -> ?from:float -> ?range:float -> int -> int -> ('a,'b) t
(** Creates a random matrix, similarly to [Vector.random] *)
@ -163,25 +169,25 @@ val map_inplace: (float -> float) -> b:('a,'b) t -> ('a,'b) t -> unit
val scale: float -> ('a,'b) t -> ('a,'b) t
(** Multiplies the matrix by a constant *)
val scale_inplace: float -> ('a,'b) t -> unit
(** Multiplies the matrix by a constant *)
val scale_cols: ('a,'b) t -> 'b Vector.t -> ('a,'b) t
(** Multiplies the matrix by a constant *)
val scale_cols_inplace: ('a,'b) t -> 'b Vector.t -> unit
(** Multiplies the matrix by a constant *)
val scale_rows: 'a Vector.t -> ('a,'b) t -> ('a,'b) t
(** Multiplies the matrix by a constant *)
val scale_rows_inplace: 'a Vector.t -> ('a,'b) t -> unit
(** Multiplies the matrix by a constant *)
val sycon: ('a,'b) t -> float
(** Returns the condition number of a matrix *)
val outer_product : ?alpha:float -> 'a Vector.t -> 'b Vector.t -> ('a,'b) t
(** Computes M = %{ $\alpha u.v^t$ %} *)
@ -190,7 +196,7 @@ val outer_product_inplace : ('a,'b) t -> ?alpha:float -> 'a Vector.t -> 'b Vecto
val gemv_n_inplace : ?m:int -> ?n:int -> ?beta:float -> 'a Vector.t ->
?alpha:float -> ?ar:int -> ?ac:int -> ('a,'b) t -> 'b Vector.t ->
?alpha:float -> ?ar:int -> ?ac:int -> ('a,'b) t -> 'b Vector.t ->
unit
(** Performs the Lapack GEMV operation. Default values:
[beta=0.] [alpha=1.0].
@ -248,43 +254,43 @@ val gemm_tn_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
val gemm: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
?c:('a,'b) t -> ?transa:[`N | `T] -> ?alpha:float ->
('c,'d) t -> ?transb:[`N | `T] -> ('e,'f) t -> ('a,'b) t
('c,'d) t -> ?transb:[`N | `T] -> ('e,'f) t -> ('a,'b) t
(** Performs the Lapack GEMM operation. Default values:
[beta=0.] [transa=`N] [alpha=1.0] [transb=`N]
[beta=0.] [transa=`N] [alpha=1.0] [transb=`N]
[gemm ~beta ~alpha a b]: %{ $C = \beta C + \alpha A B$ *)
val gemm_nn: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('b,'c) t -> ('a,'c) t
?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('b,'c) t -> ('a,'c) t
(** Performs gemm with [~transa=`N] and [~transb=`N]. *)
val gemm_nt: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> ('a,'c) t
?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> ('a,'c) t
(** Performs gemm with [~transa=`N] and [~transb=`T]. *)
val gemm_tn: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('b,'c) t -> ('a,'c) t
val gemm_tn: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('b,'c) t -> ('a,'c) t
(** Performs gemm with [~transa=`T] and [~transb=`N]. *)
val gemm_tt: ?m:int -> ?n:int -> ?k:int -> ?beta:float ->
?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> ('a,'c) t
?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> ('a,'c) t
(** Performs gemm with [~transa=`T] and [~transb=`T]. *)
val gemm_trace: ?transa:[`N | `T] -> ('a,'b) t -> ?transb:[`N | `T] -> ('c,'d) t -> float
val gemm_trace: ?transa:[`N | `T] -> ('a,'b) t -> ?transb:[`N | `T] -> ('c,'d) t -> float
(** Computes the trace of a GEMM. Default values:
[transa=`N] [transb=`N]
[transa=`N] [transb=`N]
[gemm_trace a b]: %{ $C = Tr(A B)$ *)
val gemm_nn_trace: ('a,'b) t -> ('b,'c) t -> float
val gemm_nn_trace: ('a,'b) t -> ('b,'c) t -> float
(** Computes the trace of a GEMM with [~transa=`N] and [~transb=`N]. *)
val gemm_nt_trace: ('a,'b) t -> ('c,'b) t -> float
val gemm_nt_trace: ('a,'b) t -> ('c,'b) t -> float
(** Computes the trace of a GEMM with [~transa=`N] and [~transb=`T]. *)
val gemm_tn_trace: ('b,'a) t -> ('b,'c) t -> float
val gemm_tn_trace: ('b,'a) t -> ('b,'c) t -> float
(** Computes the trace of a GEMM with [~transa=`T] and [~transb=`N]. *)
val gemm_tt_trace: ('b,'a) t -> ('c,'b) t -> float
val gemm_tt_trace: ('b,'a) t -> ('c,'b) t -> float
(** Computes the trace of a GEMM with [~transa=`T] and [~transb=`T]. *)
val svd: ('a,'b) t -> ('a,'b) t * 'b Vector.t * ('b,'b) t
@ -295,7 +301,7 @@ val svd': ('a,'b) t -> ('a,'a) t * 'a Vector.t * ('a,'b) t
val qr: ('a,'b) t -> ('a,'b) t * ('b,'b) t
(** QR factorization *)
val normalize_mat : ('a,'b) t -> ('a,'b) t
(** Returns the matrix with all the column vectors normalized *)
@ -317,7 +323,7 @@ val xt_o_x : o:('a,'a) t -> x:('a,'b) t -> ('b,'b) t
val x_o_xt : o:('b,'b) t -> x:('a,'b) t -> ('a,'a) t
(** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *)
val debug_matrix: string -> ('a,'b) t -> unit
val debug_matrix: string -> ('a,'b) t -> unit
(** Prints a matrix in stdout for debug *)
val of_file : string -> ('a,'b) t
@ -332,10 +338,10 @@ val sym_of_file : string -> ('a,'b) t
val sysv_inplace : b:('a,'b) t -> ('a,'a) t -> unit
(** Solves %{ $AX=B$ %} when A is symmetric, and stores the result in B. *)
val sysv : b:('a,'b) t -> ('a,'a) t -> ('a,'b) t
(** Solves %{ $AX=B$ %} when A is symmetric *)
val (%:) : ('a,'b) t -> (int*int) -> float
(** [t%:(i,j)] returns the element at i,j. *)
@ -344,6 +350,6 @@ val set : ('a,'b) t -> int -> int -> float -> unit
val to_file : filename:string -> ?sym:bool -> ?cutoff:float -> ('a,'b) t -> unit
(** Write the matrix to a file. *)
val pp : Format.formatter -> ('a,'b) t -> unit

16
linear_algebra/lib/vector.ml
View File

@ -1,10 +1,14 @@
open Lacaml.D
type 'a t = Vec.t
let relabel t = t
let copy ?n ?ofsy ?incy ?y ?ofsx ?incx t = copy ?n ?ofsy ?incy ?y ?ofsx ?incx t
let norm t = nrm2 t
let copy ?n ?ofsy ?incy ?y ?ofsx ?incx t = copy ?n ?ofsy ?incy ?y ?ofsx ?incx t
let norm ?(l=`L2) t =
match l with
| `L2 -> nrm2 t
| `L1 -> asum t
| `Linf -> amax t
let dim t = Vec.dim t
let neg t = Vec.neg t
@ -50,8 +54,8 @@ let to_bigarray t = copy t
let of_bigarray_inplace t = t
let to_bigarray_inplace t = t
let random ?rnd_state ?(from= -. 1.0) ?(range=2.0) n =
let random ?rnd_state ?(from= -. 1.0) ?(range=2.0) n =
let state =
match rnd_state with
| None -> Random.get_state ()
@ -67,5 +71,5 @@ let normalize v =
let (%.) t i = t.{i}
let set t i v = t.{i} <- v

38
linear_algebra/lib/vector.mli
View File

@ -31,58 +31,58 @@ val div : 'a t -> 'a t -> 'a t
val dot : 'a t -> 'a t -> float
(** [dot v1 v2] : Dot product between v1 and v2 *)
val norm : 'a t -> float
(** Norm of the vector : %{ ||v|| = $\sqrt{v.v}$ %} *)
val norm : ?l:[< `L1 | `L2 | `Linf > `L2 ] -> 'a t -> float
(** $L^1$, $L^2$ or $L^\infty$ norm of the vector *)
val sqr : 'a t -> 'a t
(** [sqr t = map (fun x -> x *. x) t] *)
val sqrt : 'a t -> 'a t
(** [sqrt t = map sqrt t] *)
val neg : 'a t -> 'a t
(** [neg t = map (f x -> -. x) t *)
val reci : 'a t -> 'a t
(** [reci t = map (f x -> 1./x) t *)
val abs : 'a t -> 'a t
(** [abs t = map (f x -> abs_float x) t *)
val cos : 'a t -> 'a t
(** [cos t = map (f x -> cos x) t *)
val sin : 'a t -> 'a t
(** [sin t = map (f x -> sin x) t *)
val tan : 'a t -> 'a t
(** [tan t = map (f x -> tan x) t *)
val acos : 'a t -> 'a t
(** [acos t = map (f x -> acos x) t *)
val asin : 'a t -> 'a t
(** [asin t = map (f x -> asin x) t *)
val atan : 'a t -> 'a t
(** [atan t = map (f x -> atan x) t *)
val amax : 'a t -> float
(** Maximum of the absolute values of the elements of the vector. *)
val normalize : 'a t -> 'a t
(** Returns the vector normalized *)
val init : int -> (int -> float) -> 'a t
val init : int -> (int -> float) -> 'a t
(** Initialize the vector with a function taking the index as parameter.*)
val sum : 'a t -> float
(** Returns the sum of the elements of the vector *)
val copy : ?n:int -> ?ofsy:int -> ?incy:int -> ?y:vec -> ?ofsx:int -> ?incx:int -> 'a t -> 'a t
(** Returns a copy of the vector X into Y. [ofs] controls the offset and [inc]
the increment. *)
val map : (float -> float) -> 'a t -> 'a t
(** Equivalent to [Array.map] *)
@ -106,13 +106,13 @@ val to_list : 'a t -> float list
val create : int -> 'a t
(** Creates an uninitialized vector *)
val make : int -> float -> 'a t
(** Creates a vector initialized with the given [float].*)
val make0 : int -> 'a t
(** Creates a zero-initialized vector.*)
val fold : ('a -> float -> 'a) -> 'a -> 'b t -> 'a
(** Equivalent to [Array.fold] *)
@ -141,4 +141,4 @@ val (%.) : 'a t -> int -> float
val set : 'a t -> int -> float -> unit
(** Modifies the value in-place at the i-th position *)