From 79f3b5ef10ab074197b59f04d935297533c353ea Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 2 Oct 2020 18:55:19 +0200 Subject: [PATCH] Added phantom types to matrices --- common/lib/util.mli | 2 + gaussian_integrals/lib/electron_nucleus.ml | 8 +- gaussian_integrals/lib/electron_nucleus.mli | 4 +- gaussian_integrals/lib/kinetic.ml | 6 +- gaussian_integrals/lib/matrix_on_basis.mli | 10 +- gaussian_integrals/lib/multipole.ml | 2 +- gaussian_integrals/lib/multipole.mli | 30 +-- gaussian_integrals/lib/orthonormalization.ml | 15 +- gaussian_integrals/lib/orthonormalization.mli | 2 +- gaussian_integrals/lib/overlap.ml | 6 +- .../lib/two_electron_integrals.mli | 2 +- linear_algebra/lib/diis.ml | 15 +- linear_algebra/lib/four_idx_storage.ml | 16 +- linear_algebra/lib/four_idx_storage.mli | 36 ++-- linear_algebra/lib/matrix.ml | 47 ++++- linear_algebra/lib/matrix.mli | 191 ++++++++++++------ linear_algebra/lib/orthonormalization.ml | 5 +- linear_algebra/lib/orthonormalization.mli | 5 +- linear_algebra/lib/spherical_to_cartesian.ml | 2 + linear_algebra/lib/spherical_to_cartesian.mli | 5 +- linear_algebra/lib/vector.ml | 1 + linear_algebra/lib/vector.mli | 4 +- 22 files changed, 265 insertions(+), 149 deletions(-) diff --git a/common/lib/util.mli b/common/lib/util.mli index 89cbb75..316bb64 100644 --- a/common/lib/util.mli +++ b/common/lib/util.mli @@ -1,5 +1,7 @@ (** All utilities which should be included in all source files are defined here *) + + (** {2 Functions from libm} *) external erf_float : float -> float = "erf_float_bytecode" "erf_float" diff --git a/gaussian_integrals/lib/electron_nucleus.ml b/gaussian_integrals/lib/electron_nucleus.ml index fef1db9..71edfc0 100644 --- a/gaussian_integrals/lib/electron_nucleus.ml +++ b/gaussian_integrals/lib/electron_nucleus.ml @@ -7,14 +7,14 @@ include Qcaml_gaussian_basis open Util open Constants -type t = Matrix.t -external matrix : t -> Matrix.t = "%identity" -external of_matrix : Matrix.t -> t = "%identity" - module Am = Angular_momentum module Bs = Basis module Cs = Contracted_shell +type t = (Bs.t, Bs.t) Matrix.t +external matrix : t -> (Bs.t, Bs.t) Matrix.t = "%identity" +external of_matrix : (Bs.t, Bs.t) Matrix.t -> t = "%identity" + (** (0|0)^m : Fundamental electron-nucleus repulsion integral $ \int \phi_p(r1) 1/r_{C} dr_1 $ diff --git a/gaussian_integrals/lib/electron_nucleus.mli b/gaussian_integrals/lib/electron_nucleus.mli index 3d10e26..12ee05f 100644 --- a/gaussian_integrals/lib/electron_nucleus.mli +++ b/gaussian_integrals/lib/electron_nucleus.mli @@ -6,11 +6,11 @@ $$ %} *) +include module type of Matrix_on_basis + open Qcaml_particles open Qcaml_gaussian_basis -include module type of Matrix_on_basis - val of_basis_nuclei : basis:Basis.t -> Nuclei.t -> t (** Build from a {Basis.t} and the nuclei (coordinates and charges). *) diff --git a/gaussian_integrals/lib/kinetic.ml b/gaussian_integrals/lib/kinetic.ml index 5e28e5a..b2bfc7a 100644 --- a/gaussian_integrals/lib/kinetic.ml +++ b/gaussian_integrals/lib/kinetic.ml @@ -13,9 +13,9 @@ module Po = Powers module Psp = Primitive_shell_pair module Ps = Primitive_shell -type t = Matrix.t -external matrix : t -> Matrix.t = "%identity" -external of_matrix : Matrix.t -> t = "%identity" +type t = (Basis.t, Basis.t) Matrix.t +external matrix : t -> (Basis.t, Basis.t) Matrix.t = "%identity" +external of_matrix : (Basis.t, Basis.t) Matrix.t -> t = "%identity" let cutoff = integrals_cutoff diff --git a/gaussian_integrals/lib/matrix_on_basis.mli b/gaussian_integrals/lib/matrix_on_basis.mli index 16f4edd..fe89439 100644 --- a/gaussian_integrals/lib/matrix_on_basis.mli +++ b/gaussian_integrals/lib/matrix_on_basis.mli @@ -3,7 +3,7 @@ open Qcaml_gaussian_basis open Qcaml_linear_algebra -type t = Matrix.t +type t = (Basis.t, Basis.t) Matrix.t val of_basis : Basis.t -> t (** Build from a {!Basis.t}. *) @@ -14,8 +14,8 @@ val of_basis_pair : Basis.t -> Basis.t -> t val to_file : filename:string -> t -> unit (** Write the integrals in a file. *) -val matrix : t -> Matrix.t -(** Returns the matrix suitable for Lacaml. *) +val matrix : t -> (Basis.t, Basis.t) Matrix.t +(** Returns the matrix. *) -val of_matrix : Matrix.t -> t -(** Build from a Lacaml matrix. *) +val of_matrix : (Basis.t, Basis.t) Matrix.t -> t +(** Build from a matrix. *) diff --git a/gaussian_integrals/lib/multipole.ml b/gaussian_integrals/lib/multipole.ml index 5cb99ef..7668fbb 100644 --- a/gaussian_integrals/lib/multipole.ml +++ b/gaussian_integrals/lib/multipole.ml @@ -3,7 +3,7 @@ open Qcaml_linear_algebra open Qcaml_gaussian_basis open Constants -type t = Matrix.t array +type t = (Basis.t, Basis.t) Matrix.t array (* [| "x"; "y"; "z"; "x2"; "y2"; "z2" |] *) diff --git a/gaussian_integrals/lib/multipole.mli b/gaussian_integrals/lib/multipole.mli index 8e576b4..5672048 100644 --- a/gaussian_integrals/lib/multipole.mli +++ b/gaussian_integrals/lib/multipole.mli @@ -14,49 +14,49 @@ open Qcaml_gaussian_basis type t -val matrix_x : t -> Matrix.t +val matrix_x : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *) -val matrix_y : t -> Matrix.t +val matrix_y : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *) -val matrix_z : t -> Matrix.t +val matrix_z : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *) -val matrix_x2 : t -> Matrix.t +val matrix_x2 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *) -val matrix_xy : t -> Matrix.t +val matrix_xy : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | xy | \chi_j \rangle $$ %} *) -val matrix_yz : t -> Matrix.t +val matrix_yz : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | yz | \chi_j \rangle $$ %} *) -val matrix_xz : t -> Matrix.t +val matrix_xz : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | xz | \chi_j \rangle $$ %} *) -val matrix_y2 : t -> Matrix.t +val matrix_y2 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *) -val matrix_z2 : t -> Matrix.t +val matrix_z2 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *) -val matrix_x3 : t -> Matrix.t +val matrix_x3 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | x^3 | \chi_j \rangle $$ %} *) -val matrix_y3 : t -> Matrix.t +val matrix_y3 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | y^3 | \chi_j \rangle $$ %} *) -val matrix_z3 : t -> Matrix.t +val matrix_z3 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | z^3 | \chi_j \rangle $$ %} *) -val matrix_x4 : t -> Matrix.t +val matrix_x4 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | x^4 | \chi_j \rangle $$ %} *) -val matrix_y4 : t -> Matrix.t +val matrix_y4 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | y^4 | \chi_j \rangle $$ %} *) -val matrix_z4 : t -> Matrix.t +val matrix_z4 : t -> (Basis.t, Basis.t) Matrix.t (** {% $$ \langle \chi_i | z^4 | \chi_j \rangle $$ %} *) val of_basis : Basis.t -> t diff --git a/gaussian_integrals/lib/orthonormalization.ml b/gaussian_integrals/lib/orthonormalization.ml index 4fd4c89..85b478a 100644 --- a/gaussian_integrals/lib/orthonormalization.ml +++ b/gaussian_integrals/lib/orthonormalization.ml @@ -2,13 +2,13 @@ open Qcaml_common open Qcaml_linear_algebra open Qcaml_gaussian_basis -type t = Matrix.t - module Am = Angular_momentum module Bs = Basis module Cs = Contracted_shell module Ov = Overlap +type t = (Bs.t, Bs.t) Matrix.t + let make_canonical ~thresh ~basis ~cartesian ~overlap = @@ -17,21 +17,22 @@ let make_canonical ~thresh ~basis ~cartesian ~overlap = let make_canonical_spherical basis = let ao_num = Bs.size basis in - let cart_sphe = Matrix.make ao_num ao_num 0. + let cart_sphe : (Bs.t, Bs.t) Matrix.t = Matrix.make ao_num ao_num 0. and i = ref 0 and n = ref 0 in Array.iter (fun shell -> let submatrix = Spherical_to_cartesian.matrix (Cs.ang_mom shell) + |> Matrix.relabel in Matrix.copy_inplace ~b:cart_sphe ~br:(!i+1) ~bc:(!n+1) submatrix; i := !i + Matrix.dim1 submatrix; n := !n + Matrix.dim2 submatrix; ) (Bs.contracted_shells basis); - let s = Matrix.gemm ~transa:`T ~m:!n cart_sphe overlap_matrix in - let overlap_matrix = Matrix.gemm s ~n:!n cart_sphe in + let s = Matrix.gemm_tn ~m:!n cart_sphe overlap_matrix in + let overlap_matrix = Matrix.gemm_nn s ~n:!n cart_sphe in let s = Orthonormalization.canonical_ortho ~thresh ~overlap:overlap_matrix (Matrix.identity !n) in - Matrix.gemm cart_sphe ~k:!n s + Matrix.gemm_nn cart_sphe ~k:!n s in if cartesian then @@ -61,7 +62,7 @@ let make_lowdin ~thresh ~overlap = Matrix.init_cols (Matrix.dim1 u_vec) (Matrix.dim2 u_vec) (fun i j -> u_vec_x.{i,j} *. (Vector.at u_val j)) in - Matrix.gemm u_vec' ~transb:`T u_vec + Matrix.gemm_nt u_vec' u_vec diff --git a/gaussian_integrals/lib/orthonormalization.mli b/gaussian_integrals/lib/orthonormalization.mli index 98f4caa..bc69a4c 100644 --- a/gaussian_integrals/lib/orthonormalization.mli +++ b/gaussian_integrals/lib/orthonormalization.mli @@ -3,7 +3,7 @@ open Qcaml_gaussian_basis open Qcaml_linear_algebra -type t = Matrix.t +type t = (Basis.t, Basis.t) Matrix.t val make: ?thresh:float -> ?basis:Basis.t -> cartesian:bool -> Overlap.t -> t (** Returns a matrix or orthonormal vectors in the basis. The vectors are diff --git a/gaussian_integrals/lib/overlap.ml b/gaussian_integrals/lib/overlap.ml index 16a0803..cd1f3f1 100644 --- a/gaussian_integrals/lib/overlap.ml +++ b/gaussian_integrals/lib/overlap.ml @@ -12,9 +12,9 @@ module Csp = Contracted_shell_pair module Po = Powers module Psp = Primitive_shell_pair -type t = Matrix.t -external matrix : t -> Matrix.t = "%identity" -external of_matrix : Matrix.t -> t = "%identity" +type t = (Basis.t, Basis.t) Matrix.t +external matrix : t -> (Basis.t, Basis.t) Matrix.t = "%identity" +external of_matrix : (Basis.t, Basis.t) Matrix.t -> t = "%identity" let cutoff = integrals_cutoff diff --git a/gaussian_integrals/lib/two_electron_integrals.mli b/gaussian_integrals/lib/two_electron_integrals.mli index 83ce1f9..af36f9e 100644 --- a/gaussian_integrals/lib/two_electron_integrals.mli +++ b/gaussian_integrals/lib/two_electron_integrals.mli @@ -28,7 +28,7 @@ module Make : functor (T : Two_ei_structure) -> sig include module type of Four_idx_storage - val of_basis : ?operator:Operator.t -> Basis.t -> t + val of_basis : ?operator:Operator.t -> Basis.t -> Basis.t t (** Compute all ERI's for a given {!Basis.t}. *) end diff --git a/linear_algebra/lib/diis.ml b/linear_algebra/lib/diis.ml index 5ea7f4d..9bd02c7 100644 --- a/linear_algebra/lib/diis.ml +++ b/linear_algebra/lib/diis.ml @@ -32,6 +32,9 @@ let append ~p ~e diis = } +type nvec +type one + let next diis = let e = Matrix.of_col_vecs_list diis.e and p = Matrix.of_col_vecs_list diis.p @@ -41,7 +44,7 @@ let next diis = let a = Matrix.make (m+1) (m+1) 1. in let ax = Matrix.to_bigarray_inplace a in ax.{m+1,m+1} <- 0.; - Matrix.gemm_inplace ~c:a ~transa:`T ~m ~n:m e e; + Matrix.gemm_tn_inplace ~c:a ~m ~n:m e e; if m > 1 && Matrix.sycon a > 1.e-14 then (aux [@tailcall]) (m-1) else a @@ -49,12 +52,16 @@ let next diis = aux diis.m in let m = Matrix.dim1 a - 1 in - let c = Matrix.make0 (m+1) 1 in + let c : (nvec,one) Matrix.t = Matrix.make0 (m+1) 1 in let cx = Matrix.to_bigarray_inplace c in cx.{m+1,1} <- 1.; - Matrix.sysv_inplace ~b:c a; - Matrix.gemm p c ~k:m + Matrix.sysv_inplace a ~b:c ; + Matrix.gemm_nn ~k:m p c |> Matrix.as_vec_inplace + |> Vector.relabel +(* + |> Vector.relabel + *) diff --git a/linear_algebra/lib/four_idx_storage.ml b/linear_algebra/lib/four_idx_storage.ml index 3d3f517..8af21e7 100644 --- a/linear_algebra/lib/four_idx_storage.ml +++ b/linear_algebra/lib/four_idx_storage.ml @@ -4,18 +4,18 @@ let max_index = 1 lsl 14 type index_pair = { first : int ; second : int } -type storage_t = - | Dense of Matrix.t +type 'a storage_t = + | Dense of ('a,'a) Matrix.t | Sparse of (int, float) Hashtbl.t -type t = +type 'a t = { size : int ; - two_index : Matrix.t; - two_index_anti : Matrix.t; - three_index : Matrix.t; - three_index_anti : Matrix.t; - four_index : storage_t ; + two_index : ('a,'a) Matrix.t; + two_index_anti : ('a,'a) Matrix.t; + three_index : ('a,'a) Matrix.t; + three_index_anti : ('a,'a) Matrix.t; + four_index : 'a storage_t ; } let key_of_indices ~r1 ~r2 = diff --git a/linear_algebra/lib/four_idx_storage.mli b/linear_algebra/lib/four_idx_storage.mli index 24e9e1c..66cd1d8 100644 --- a/linear_algebra/lib/four_idx_storage.mli +++ b/linear_algebra/lib/four_idx_storage.mli @@ -7,7 +7,7 @@ There are two kinds of ordering of indices: *) -type t +type 'a t (** Element for the stream *) type element = @@ -19,61 +19,61 @@ type element = value: float } -val create : size:int -> [< `Dense | `Sparse ] -> t +val create : size:int -> [< `Dense | `Sparse ] -> 'a t (** If [`Dense] is chosen, internally the data is stored as a 4-dimensional [Bigarray]. Else, it is stored as a hash table. *) (** {2 Accessors} *) -val size : t -> int +val size : 'a t -> int (** Number of stored elements *) -val get_chem : t -> int -> int -> int -> int -> float +val get_chem : 'a t -> int -> int -> int -> int -> float (** Get an integral using the Chemist's convention {% $(ij|kl)$ %}. *) -val get_phys : t -> int -> int -> int -> int -> float +val get_phys : 'a t -> int -> int -> int -> int -> float (** Get an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) -val set_chem : t -> int -> int -> int -> int -> float -> unit +val set_chem : 'a t -> int -> int -> int -> int -> float -> unit (** Set an integral using the Chemist's convention {% $(ij|kl)$ %}. *) -val set_phys : t -> int -> int -> int -> int -> float -> unit +val set_phys : 'a t -> int -> int -> int -> int -> float -> unit (** Set an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) -val increment_chem : t -> int -> int -> int -> int -> float -> unit +val increment_chem : 'a t -> int -> int -> int -> int -> float -> unit (** Increments an integral using the Chemist's convention {% $(ij|kl)$ %}. *) -val increment_phys : t -> int -> int -> int -> int -> float -> unit +val increment_phys : 'a t -> int -> int -> int -> int -> float -> unit (** Increments an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) -val get_chem_all_i : t -> j:int -> k:int -> l:int -> 'a Vector.t +val get_chem_all_i : 'a t -> j:int -> k:int -> l:int -> 'a Vector.t (** Get all integrals in an array [a.{i} =] {% $(\cdot j|kl)$ %} . *) -val get_phys_all_i : t -> j:int -> k:int -> l:int -> 'a Vector.t +val get_phys_all_i : 'a t -> j:int -> k:int -> l:int -> 'a Vector.t (** Get all integrals in an array [a.{i} =] {% $\langle \cdot j|kl \rangle$ %} . *) -val get_chem_all_ij : t -> k:int -> l:int -> Matrix.t +val get_chem_all_ij : 'a t -> k:int -> l:int -> ('a,'a) Matrix.t (** Get all integrals in an array [a.{i,j} =] {% $(\cdot \cdot|kl)$ %} . *) -val get_phys_all_ij : t -> k:int -> l:int -> Matrix.t +val get_phys_all_ij : 'a t -> k:int -> l:int -> ('a,'a) Matrix.t (** Get all integrals in an array [a.{i,j} =] {% $\langle \cdot \cdot|kl \rangle$ %} . *) -val to_stream : t -> element Stream.t +val to_stream : 'a t -> element Stream.t (** Retrun the data structure as a stream. *) -val to_list : t -> element list +val to_list : 'a t -> element list (** Retrun the data structure as a list. *) -val four_index_transform : Matrix.t -> t -> t +val four_index_transform : ('a,'b) Matrix.t -> 'a t -> 'b t (** Perform a four-index transformation *) (** {2 I/O} *) -val to_file : ?cutoff:float -> filename:string -> t -> unit +val to_file : ?cutoff:float -> filename:string -> 'a t -> unit (** Write the data to file, using the physicist's ordering. *) val of_file : size:int -> sparsity:[< `Dense | `Sparse ] - -> Scanf.Scanning.file_name -> t + -> Scanf.Scanning.file_name -> 'a t (** Read from a text file with format ["%d %d %d %d %f"]. *) diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml index df13691..c634799 100644 --- a/linear_algebra/lib/matrix.ml +++ b/linear_algebra/lib/matrix.ml @@ -1,6 +1,6 @@ open Lacaml.D -type t = Mat.t +type ('a, 'b) t = Mat.t let dim1 t = Mat.dim1 t let dim2 t = Mat.dim2 t @@ -10,6 +10,8 @@ let out_of_place f t = f result; 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 @@ -41,9 +43,9 @@ let add_const x a = let at t i j = t.{i,j} -external to_bigarray_inplace : t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t = "%identity" +external to_bigarray_inplace : ('a,'b) t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t = "%identity" -external of_bigarray_inplace : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> t = "%identity" +external of_bigarray_inplace : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> ('a,'b) t = "%identity" let to_bigarray t = lacpy t let of_bigarray t = lacpy t @@ -99,9 +101,23 @@ let scale_inplace x t = let scale x t = out_of_place (fun t -> scale_inplace x t) t + let gemm_inplace ?m ?n ?k ?(beta=0.) ~c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b = ignore @@ gemm ?m ?n ?k ~beta ~c ~transa ~alpha a ~transb b +let gemm_nn_inplace ?m ?n ?k ?(beta=0.) ~c ?(alpha=1.0) a b = + gemm_inplace ?m ?n ?k ~beta ~c ~alpha a b ~transa:`N ~transb:`N + +let gemm_nt_inplace ?m ?n ?k ?(beta=0.) ~c ?(alpha=1.0) a b = + gemm_inplace ?m ?n ?k ~beta ~c ~alpha a b ~transa:`N ~transb:`T + +let gemm_tn_inplace ?m ?n ?k ?(beta=0.) ~c ?(alpha=1.0) a b = + gemm_inplace ?m ?n ?k ~beta ~c ~alpha a b ~transa:`T ~transb:`N + +let gemm_tt_inplace ?m ?n ?k ?(beta=0.) ~c ?(alpha=1.0) a b = + gemm_inplace ?m ?n ?k ~beta ~c ~alpha a b ~transa:`T ~transb:`T + + let gemm ?m ?n ?k ?beta ?c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b = let c = match c with @@ -110,10 +126,27 @@ let gemm ?m ?n ?k ?beta ?c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b = in gemm ?m ?n ?k ?beta ?c ~transa ~alpha a ~transb b +let gemm_nn ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`N ~transb:`N + +let gemm_nt ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`N ~transb:`T + +let gemm_tn ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`T ~transb:`N + +let gemm_tt ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`T ~transb:`T + 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 init_cols = Mat.init_cols let of_col_vecs a = @@ -244,7 +277,13 @@ let scale_cols a v = a' -let svd a = +let svd a = + let d, u, vt = + gesvd (lacpy a) + in + u, (Vector.of_bigarray d), vt + +let svd_t a = let d, u, vt = gesvd (lacpy a) in diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli index 712f606..35ff43e 100644 --- a/linear_algebra/lib/matrix.mli +++ b/linear_algebra/lib/matrix.mli @@ -1,205 +1,262 @@ -type t +(** Type for matrices. The ['a] and ['b] types are labels for the rows and columns. *) -val dim1: t -> int +type ('a,'b) t + +val dim1: ('a,'b) t -> int (** First dimension of the matrix *) -val dim2: t -> int +val dim2: ('a,'b) t -> int (** Second dimension of the matrix *) -val make : int -> int -> float -> t +val make : int -> int -> float -> ('a,'b) t (** Creates a matrix initialized with the given value. *) -val make0 : int -> int -> t +val make0 : int -> int -> ('a,'b) t (** Creates a zero-initialized matrix. *) -val create : int -> int -> t +val create : int -> int -> ('a,'b) t (** Creates an uninitialized matrix. *) -val reshape : t -> int -> int -> t +val reshape : ('a,'b) t -> int -> int -> ('c,'d) t (** Changes the dimensions of the matrix *) -val init_cols : int -> int -> (int -> int -> float) -> t +val init_cols : int -> int -> (int -> int -> float) -> ('a,'b) t (** Creates an uninitialized matrix. *) -val identity: int -> t +val identity: int -> ('a,'b) t (** Creates an identity matrix. *) -val fill_inplace: t -> float -> unit +val fill_inplace: ('a,'b) t -> float -> unit (** Fills the matrix with the give value. *) -val add_const_diag : float -> t -> t +val add_const_diag : float -> ('a,'b) t -> ('a,'b) t (** Adds a constant to the diagonal *) -val add_const_diag_inplace : float -> t -> unit +val add_const_diag_inplace : float -> ('a,'b) t -> unit (** Adds a constant to the diagonal *) -val add_const_inplace : float -> t -> unit +val add_const_inplace : float -> ('a,'b) t -> unit (** Adds a constant to the diagonal *) -val add_const : float -> t -> t +val add_const : float -> ('a,'b) t -> ('a,'b) t (** Adds a constant to the diagonal *) -val add : t -> t -> t +val add : ('a,'b) t -> ('a,'b) t -> ('a,'b) t (** Adds two matrices *) -val sub : t -> t -> t +val sub : ('a,'b) t -> ('a,'b) t -> ('a,'b) t (** Subtracts two matrices *) -val mul : t -> t -> t +val mul : ('a,'b) t -> ('a,'b) t -> ('a,'b) t (** Multiplies two matrices element-wise *) -val div : t -> t -> t +val div : ('a,'b) t -> ('a,'b) t -> ('a,'b) t (** Divides two matrices element-wise *) -val add_inplace : c:t -> t -> t -> unit +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:t -> t -> t -> unit +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:t -> t -> t -> unit +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:t -> t -> t -> unit +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 at : t -> int -> int -> float +val at : ('a,'b) t -> int -> int -> float (** [at i j] returns the element at i,j. *) -val to_bigarray : t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t +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 to_bigarray_inplace : - t -> (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 (** Converts the matrix into a Bigarray in Fortran layout in place*) -val to_col_vecs : t -> 'a Vector.t array +val to_col_vecs : ('a,'b) t -> 'a Vector.t array (** Converts the matrix into an array of vectors *) -val to_col_vecs_list : t -> 'a Vector.t list +val to_col_vecs_list : ('a,'b) t -> 'a Vector.t list (** Converts the matrix into a list of vectors *) -val of_col_vecs : 'a Vector.t array -> t +val of_col_vecs : 'a Vector.t array -> ('a,'b) t (** Converts an array of vectors into a matrix *) -val of_col_vecs_list : 'a Vector.t list -> t +val of_col_vecs_list : 'a Vector.t list -> ('a,'b) t (** Converts a list of vectors into a matrix *) -val of_bigarray : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> 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 *) val of_bigarray_inplace : - (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> 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 copy: ?m:int -> ?n:int -> ?br:int -> ?bc:int -> ?ar:int -> ?ac:int -> t -> t +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 *) -val copy_inplace: ?m:int -> ?n:int -> ?br:int -> ?bc:int -> b:t -> ?ar:int -> ?ac:int -> t -> unit +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 *) -val col: t -> int -> 'a Vector.t +val col: ('a,'b) t -> int -> 'a Vector.t (** Returns a column of the matrix as a vector *) -val detri: t -> 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: 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: t -> 'a Vector.t +val as_vec_inplace: ('a,'b) t -> ('a*'b) Vector.t (** Interpret the matrix as a vector (reshape). *) -val as_vec: t -> 'a Vector.t +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 -> t +val random: ?rnd_state:Random.State.t -> ?from:float -> ?range:float -> int -> int -> ('a,'b) t (** Creates a random matrix, similarly to [Vector.random] *) -val map: (float -> float) -> t -> t +val map: (float -> float) -> ('a,'b) t -> ('a,'b) t (** Apply the function to all elements of the matrix *) -val map_inplace: (float -> float) -> b:t -> t -> unit +val map_inplace: (float -> float) -> b:('a,'b) t -> ('a,'b) t -> unit (** [map_inplace f b a] : Apply the function to all elements of the matrix [a] and store the results in [b] *) -val scale: float -> t -> t +val scale: float -> ('a,'b) t -> ('a,'b) t (** Multiplies the matrix by a constant *) -val scale_inplace: float -> t -> unit +val scale_inplace: float -> ('a,'b) t -> unit (** Multiplies the matrix by a constant *) -val scale_cols: t -> 'a Vector.t -> t +val scale_cols: ('a,'b) t -> 'b Vector.t -> ('a,'b) t (** Multiplies the matrix by a constant *) -val scale_cols_inplace: t -> 'a Vector.t -> unit +val scale_cols_inplace: ('a,'b) t -> 'b Vector.t -> unit (** Multiplies the matrix by a constant *) -val sycon: t -> float +val sycon: ('a,'b) t -> float (** Returns the condition number of a matrix *) -val outer_product : ?alpha:float -> 'a Vector.t -> 'b Vector.t -> t +val outer_product : ?alpha:float -> 'a Vector.t -> 'b Vector.t -> ('a,'b) t (** Computes M = %{ $\alpha u.v^t$ %} *) -val outer_product_inplace : t -> ?alpha:float -> 'a Vector.t -> 'b Vector.t -> unit +val outer_product_inplace : ('a,'b) t -> ?alpha:float -> 'a Vector.t -> 'b Vector.t -> unit (** Computes M = %{ $\alpha u.v^t$ %} *) -val gemm_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float -> c:t -> ?transa:[`N | `T] -> ?alpha:float -> - t -> ?transb:[`N | `T] -> t -> unit + +val gemm_inplace : ?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 -> unit (** Performs the Lapack GEMM operation. Default values: [beta=0.] [transa=`N] [alpha=1.0] [transb=`N]. [gemm ~beta c ~alpha a b]: %{ $C = \beta C + \alpha A B$ *) -val gemm: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> ?c:t -> ?transa:[`N | `T] -> ?alpha:float -> - t -> ?transb:[`N | `T] -> t -> t +val gemm_nn_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float -> + c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('b,'c) t -> unit +(** Performs gemm with [~transa=`N] and [~transb=`N]. *) + +val gemm_nt_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float -> + c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> unit +(** Performs gemm with [~transa=`N] and [~transb=`T]. *) + +val gemm_tt_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float -> + c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> unit +(** Performs gemm with [~transa=`T] and [~transb=`T]. *) + +val gemm_tn_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float -> + c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('b,'c) t -> unit +(** Performs gemm with [~transa=`T] and [~transb=`N]. *) + + +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 (** Performs the Lapack GEMM operation. Default values: [beta=0.] [transa=`N] [alpha=1.0] [transb=`N] [gemm ~beta ~alpha a b]: %{ $C = \beta C + \alpha A B$ *) -val gemm_trace: ?transa:[`N | `T] -> t -> ?transb:[`N | `T] -> t -> float +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 +(** 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 -> ('c,'a) 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 +(** 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 -> ('c,'b) 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 (** Computes the trace of a GEMM. Default values: [transa=`N] [transb=`N] [gemm_trace a b]: %{ $C = Tr(A B)$ *) -val svd: t -> t * 'b Vector.t * t -(** Singular value decomposition. *) +val gemm_nn_trace: ('a,'b) t -> ('b,'c) t -> float +(** Computes the trace of a GEMM with [~transa=`N] and [~transb=`N]. *) -val qr: t -> t * t +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 +(** Computes the trace of a GEMM with [~transa=`T] and [~transb=`N]. *) + +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 +(** Singular value decomposition of A(m,n) when m >= n. *) + +val svd_t: ('a,'b) t -> ('a,'a) t * 'a Vector.t * ('a,'b) t +(** Singular value decomposition of A(m,n) when m < n. *) + +val qr: ('a,'b) t -> ('a,'b) t * ('b,'b) t (** QR factorization *) -val normalize_mat : t -> t +val normalize_mat : ('a,'b) t -> ('a,'b) t (** Returns the matrix with all the column vectors normalized *) -val normalize_mat_inplace : t -> unit +val normalize_mat_inplace : ('a,'b) t -> unit (** Returns the matrix with all the column vectors normalized *) -val diagonalize_symm : t -> t * 'a Vector.t +val diagonalize_symm : ('a,'a) t -> ('a,'a) t * 'a Vector.t (** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *) -val xt_o_x : o:t -> x:t -> t +val xt_o_x : o:('a,'a) t -> x:('a,'b) t -> ('b,'b) t (** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *) -val x_o_xt : o:t -> x:t -> 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 -> t -> unit +val debug_matrix: string -> ('a,'b) t -> unit (** Prints a matrix in stdout for debug *) -val matrix_of_file : string -> t +val matrix_of_file : string -> ('a,'b) t (** Reads a matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) -val sym_matrix_of_file : string -> t +val relabel : ('a,'b) t -> ('c,'d) t + +val sym_matrix_of_file : string -> ('a,'b) t (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) -val sysv_inplace : b:t -> t -> unit +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:t -> t -> t +val sysv : b:('a,'b) t -> ('a,'a) t -> ('a,'b) t (** Solves %{ $AX=B$ %} when A is symmetric *) -val pp : Format.formatter -> t -> unit +val pp : Format.formatter -> ('a,'b) t -> unit diff --git a/linear_algebra/lib/orthonormalization.ml b/linear_algebra/lib/orthonormalization.ml index 6fbe7a1..10937bb 100644 --- a/linear_algebra/lib/orthonormalization.ml +++ b/linear_algebra/lib/orthonormalization.ml @@ -1,4 +1,5 @@ -let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = + +let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = let u, d, _ = Matrix.svd overlap in let d_sqrt = Vector.sqrt d in let n = (* Number of non-negligible singular vectors *) @@ -14,7 +15,7 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = Printf.printf "Removed linear dependencies below %f\n" (1. /. dx.{n}) ; Matrix.scale_cols_inplace u d_inv_sq ; - Matrix.gemm c u + Matrix.gemm_nn u c let qr_ortho m = diff --git a/linear_algebra/lib/orthonormalization.mli b/linear_algebra/lib/orthonormalization.mli index 1a3e253..db20f9a 100644 --- a/linear_algebra/lib/orthonormalization.mli +++ b/linear_algebra/lib/orthonormalization.mli @@ -1,4 +1,5 @@ -val canonical_ortho: ?thresh:float -> overlap:Matrix.t -> Matrix.t -> Matrix.t +val canonical_ortho: ?thresh:float -> overlap:('a,'a) Matrix.t -> ('a,'b) Matrix.t -> + ('a,'b) Matrix.t (** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %}, and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize. @@ -12,7 +13,7 @@ $$ *) -val qr_ortho: Matrix.t -> Matrix.t +val qr_ortho: ('a,'b) Matrix.t -> ('a,'b) Matrix.t (** QR orthogonalization of the input matrix *) diff --git a/linear_algebra/lib/spherical_to_cartesian.ml b/linear_algebra/lib/spherical_to_cartesian.ml index aa7a52f..9c76407 100644 --- a/linear_algebra/lib/spherical_to_cartesian.ml +++ b/linear_algebra/lib/spherical_to_cartesian.ml @@ -2,6 +2,8 @@ open Qcaml_common module Am = Angular_momentum +type num_cartesian_ao +type num_spherical_ao let matrix = function | Am.S -> Matrix.make 1 1 1. diff --git a/linear_algebra/lib/spherical_to_cartesian.mli b/linear_algebra/lib/spherical_to_cartesian.mli index 504dc07..ae62078 100644 --- a/linear_algebra/lib/spherical_to_cartesian.mli +++ b/linear_algebra/lib/spherical_to_cartesian.mli @@ -2,7 +2,10 @@ open Qcaml_common -val matrix : Angular_momentum.t -> Matrix.t +type num_cartesian_ao +type num_spherical_ao + +val matrix : Angular_momentum.t -> (num_cartesian_ao, num_spherical_ao) Matrix.t (** Returns a transformation matrix to rotate between the basis of atom-centered spherical coordinates to x,y,z coordinates. diff --git a/linear_algebra/lib/vector.ml b/linear_algebra/lib/vector.ml index 9bd1387..d4f0c99 100644 --- a/linear_algebra/lib/vector.ml +++ b/linear_algebra/lib/vector.ml @@ -2,6 +2,7 @@ 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 diff --git a/linear_algebra/lib/vector.mli b/linear_algebra/lib/vector.mli index e9117aa..4d158d8 100644 --- a/linear_algebra/lib/vector.mli +++ b/linear_algebra/lib/vector.mli @@ -11,6 +11,8 @@ open Lacaml.D type 'a t (* Parameter ['a] defines the basis on which the vector is expanded. *) +val relabel : 'a t -> 'b t + val dim : 'a t -> int (** Returns the dimension of the vector *) @@ -111,7 +113,7 @@ val make : int -> float -> 'a t val make0 : int -> 'a t (** Creates a zero-initialized vector.*) -val fold : ('a -> float -> 'a) -> 'a -> 'a t -> 'a +val fold : ('a -> float -> 'a) -> 'a -> 'b t -> 'a (** Equivalent to [Array.fold] *) val random : ?rnd_state:Random.State.t -> ?from:float -> ?range:float -> int -> 'a t