From f4cddaf7351ca2062dfd237c1263c63b965be4da Mon Sep 17 00:00:00 2001
From: orialb <orialb1@gmail.com>
Date: Mon, 16 Sep 2019 08:22:39 +0200
Subject: [PATCH 1/2] exp function for ITensor (with dense storage)

---
 src/itensor.jl       | 19 +++++++++++++++++++
 src/storage/dense.jl |  4 ++++
 test/test_itensor.jl | 24 ++++++++++++++++++++++++
 3 files changed, 47 insertions(+)

diff --git a/src/itensor.jl b/src/itensor.jl
index 293acdc3ce..b5984a9030 100644
--- a/src/itensor.jl
+++ b/src/itensor.jl
@@ -4,6 +4,7 @@ export ITensor,
        combinedindex,
        delta,
        δ,
+       exp,
        replaceindex!,
        inds,
        isNull,
@@ -424,6 +425,24 @@ end
 
 dot(A::ITensor,B::ITensor) = scalar(dag(A)*B)
 
+import LinearAlgebra.exp
+
+"""
+    exp(A::ITensor, Lis::IndexSet)
+Compute the exponent of the tensor `A` by treating it as a matrix ``A_{lr}`` with
+the left index `l` running over all indices in `Lis` and `r` running over all
+indices not in `Lis`. Must have `dim(Lis) == dim(inds(A))/dim(Lis)` for the exponentiation to
+be defined.
+"""
+function exp(A::ITensor, Lis::IndexSet)
+  (dim(Lis) == dim(inds(A))/dim(Lis)) || throw(DimensionMismatch("dimension of the left index set `Lis` must be
+                                                                       equal to `dim(inds(A))/dim(Lis)`"))
+  A, Lis, Ris = _permute_for_factorize(A,Lis)
+  expAs = storage_exp(store(A), Lis,Ris)
+  return ITensor(inds(A),expAs)
+end
+
+
 #######################################################################
 #
 # In-place operations
diff --git a/src/storage/dense.jl b/src/storage/dense.jl
index 691fa21f0f..f6dd7c46ad 100644
--- a/src/storage/dense.jl
+++ b/src/storage/dense.jl
@@ -357,3 +357,7 @@ function storage_polar(Astore::Dense{T},
   return (Qis,Qstore,Pis,Pstore)
 end
 
+function storage_exp(As::Dense{T}, Lis,Ris) where {T}
+  expAdata = exp( reshape(data(As),dim(Lis),dim(Ris)) )
+  return Dense{T}(vec(expAdata))
+end
diff --git a/test/test_itensor.jl b/test/test_itensor.jl
index 481d912245..e2cb518a0a 100644
--- a/test/test_itensor.jl
+++ b/test/test_itensor.jl
@@ -172,6 +172,30 @@ end
   @test dot(A, B) == 11.0
 end
 
+@testset "exponentiate" begin
+  s1 = Index(2,"s1")
+  s2 = Index(2,"s2")
+  i1 = Index(2,"i1")
+  i2 = Index(2,"i2")
+  Amat = rand(2,2,2,2)
+  A = ITensor(Amat, i1,i2,s1,s2)
+
+  Aexp = exp(A,IndexSet(i1,i2))
+  Amatexp = reshape( exp(reshape(Amat,4,4)), 2,2,2,2)
+  Aexp_from_mat = ITensor(Amatexp, i1,i2,s1,s2)
+  @test Aexp ≈ Aexp_from_mat
+
+  #test that exponentiation works when indices need to be permuted
+  Aexp = exp(A,IndexSet(s1,s2))
+  Amatexp = Array( exp(  reshape(Amat,4,4))' )
+  Aexp_from_mat = ITensor(reshape(Amatexp,2,2,2,2), s1,s2,i1,i2)
+  @test Aexp ≈ Aexp_from_mat
+
+  @test_throws DimensionMismatch exp(A,IndexSet(s1))
+
+end
+
+
 @testset "add and axpy" begin
   i = Index(2,"i")
   a = [1.0; 2.0]

From f622863e9d9471d5ac64cb84c84406efec589665 Mon Sep 17 00:00:00 2001
From: orialb <orialb1@gmail.com>
Date: Tue, 17 Sep 2019 17:33:39 +0200
Subject: [PATCH 2/2] add hermitian keyword to exp

---
 src/itensor.jl       | 12 ++++++++----
 src/storage/dense.jl |  5 +++--
 test/test_itensor.jl | 11 +++++++++++
 3 files changed, 22 insertions(+), 6 deletions(-)

diff --git a/src/itensor.jl b/src/itensor.jl
index b5984a9030..7c9e374f63 100644
--- a/src/itensor.jl
+++ b/src/itensor.jl
@@ -428,21 +428,25 @@ dot(A::ITensor,B::ITensor) = scalar(dag(A)*B)
 import LinearAlgebra.exp
 
 """
-    exp(A::ITensor, Lis::IndexSet)
-Compute the exponent of the tensor `A` by treating it as a matrix ``A_{lr}`` with
+    exp(A::ITensor, Lis::IndexSet; hermitian = false)
+Compute the exponential of the tensor `A` by treating it as a matrix ``A_{lr}`` with
 the left index `l` running over all indices in `Lis` and `r` running over all
 indices not in `Lis`. Must have `dim(Lis) == dim(inds(A))/dim(Lis)` for the exponentiation to
 be defined.
+When `hermitian=true` the exponential of `Hermitian(reshape(A, dim(Lis), :))` is
+computed internally.
 """
-function exp(A::ITensor, Lis::IndexSet)
+function exp(A::ITensor, Lis::IndexSet; hermitian = false)
   (dim(Lis) == dim(inds(A))/dim(Lis)) || throw(DimensionMismatch("dimension of the left index set `Lis` must be
                                                                        equal to `dim(inds(A))/dim(Lis)`"))
   A, Lis, Ris = _permute_for_factorize(A,Lis)
-  expAs = storage_exp(store(A), Lis,Ris)
+  expAs = storage_exp(store(A), Lis,Ris, hermitian = hermitian)
   return ITensor(inds(A),expAs)
 end
 
 
+expHermitian(A::ITensor,Lis::IndexSet) = exp(A,Lis, hermitian=true)
+
 #######################################################################
 #
 # In-place operations
diff --git a/src/storage/dense.jl b/src/storage/dense.jl
index f6dd7c46ad..2f90c590eb 100644
--- a/src/storage/dense.jl
+++ b/src/storage/dense.jl
@@ -357,7 +357,8 @@ function storage_polar(Astore::Dense{T},
   return (Qis,Qstore,Pis,Pstore)
 end
 
-function storage_exp(As::Dense{T}, Lis,Ris) where {T}
-  expAdata = exp( reshape(data(As),dim(Lis),dim(Ris)) )
+function storage_exp(As::Dense{T}, Lis,Ris; hermitian=false) where {T}
+  expAdata = ( hermitian ? Array(exp(Hermitian(reshape(data(As),dim(Lis),dim(Ris))))) :
+                           exp(reshape(data(As),dim(Lis),dim(Ris))) )
   return Dense{T}(vec(expAdata))
 end
diff --git a/test/test_itensor.jl b/test/test_itensor.jl
index e2cb518a0a..430b2e4e9a 100644
--- a/test/test_itensor.jl
+++ b/test/test_itensor.jl
@@ -191,6 +191,17 @@ end
   Aexp_from_mat = ITensor(reshape(Amatexp,2,2,2,2), s1,s2,i1,i2)
   @test Aexp ≈ Aexp_from_mat
 
+  #test exponentiation when hermitian=true is used
+  Amat = reshape(Amat, 4,4)
+  Amat = reshape( Amat + Amat' + randn(4,4)*1e-10 , 2,2,2,2)
+  A = ITensor(Amat, i1,i2,s1,s2)
+  Aexp = exp(A,IndexSet(i1,i2), hermitian=true)
+  Amatexp = Array(reshape( exp(Hermitian(reshape(Amat,4,4))), 2,2,2,2))
+  Aexp_from_mat = ITensor(Amatexp, i1,i2,s1,s2)
+  @test Aexp ≈ Aexp_from_mat
+
+
+
   @test_throws DimensionMismatch exp(A,IndexSet(s1))
 
 end