From 3c27ec0c1057ca09c6e74df1ba1e24173c06f965 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 17 Dec 2020 12:55:11 +0100
Subject: [PATCH] Added exponential of matrix

---
 linear_algebra/lib/matrix.ml  | 20 +++++++++++++++++++-
 linear_algebra/lib/matrix.mli |  4 +++-
 mo/lib/basis.ml               |  1 +
 3 files changed, 23 insertions(+), 2 deletions(-)

diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml
index 74095a3..ee77ed0 100644
--- a/linear_algebra/lib/matrix.ml
+++ b/linear_algebra/lib/matrix.ml
@@ -382,7 +382,24 @@ let qr a =
   q, r
 
   
-let (%:) t (i,j) = t.{i,j}
+let exponential a =
+  assert (dim1 a = dim2 a);
+  let a = to_bigarray_inplace a in
+  let n = Mat.dim1 a in
+  let a2 = Lacaml.D.gemm a a in
+  let (lv, wr, _wi, rv) = Lacaml.D.geev a2 in
+  let tau = Vec.map (fun x -> -. sqrt x) wr in
+  let cos_tau =
+    Mat.init_cols n n (fun i j ->
+                    if i<>j then 0. else cos tau.{i})
+  in
+  let sin_tau =
+    Mat.init_cols n n (fun i j ->
+                    if i<>j then 0. else (sin tau.{i}) /. tau.{i})
+  in
+  let g = Lacaml.D.gemm in
+  Mat.add (g lv @@ g cos_tau rv) (g lv @@ g sin_tau @@ g rv a)
+  |> of_bigarray_inplace
 
 
 let to_file ~filename ?(sym=false) ?(cutoff=0.) t =
@@ -411,6 +428,7 @@ 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
 
diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli
index 46f3090..55466b8 100644
--- a/linear_algebra/lib/matrix.mli
+++ b/linear_algebra/lib/matrix.mli
@@ -279,7 +279,6 @@ val gemm_tn_trace: ('b,'a) t -> ('b,'c) 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
 (** Singular value decomposition of A(m,n) when m >= n. *)
 
@@ -298,6 +297,9 @@ val normalize_mat_inplace : ('a,'b) t -> unit
 val diagonalize_symm : ('a,'a) t -> ('a,'a) t * 'a Vector.t
 (** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *)
 
+val exponential : ('a,'a) t -> ('a,'a) t
+(** Computes the exponential of a square matrix. *)
+
 val xt_o_x : o:('a,'a) t -> x:('a,'b) t -> ('b,'b) t
 (** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *)
 
diff --git a/mo/lib/basis.ml b/mo/lib/basis.ml
index 7e8c619..a675274 100644
--- a/mo/lib/basis.ml
+++ b/mo/lib/basis.ml
@@ -109,6 +109,7 @@ let values t point =
   let a = Ao.Basis.values (Simulation.ao_basis t.simulation) point in
   Matrix.gemv_t c a
 
+
 let of_hartree_fock hf =
   let mo_coef       = HF.eigenvectors hf in
   let simulation    = HF.simulation   hf in