10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-13 17:43:41 +01:00

Added phantom types to matrices

This commit is contained in:
Anthony Scemama 2020-10-02 18:55:19 +02:00
parent 6122bf79d8
commit 79f3b5ef10
22 changed files with 265 additions and 149 deletions

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -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" |]
*)

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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