From 6b2ad3db802ea0a2faa755bd919d282ec48e469e Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 11 Aug 2026 16:31:42 -0400
Subject: [PATCH] Inverse improvements. Need to make this properly recursive

---
 Grid/lattice/Lattice_trace.h | 62 ++++++++++++++++++++++++++++++++++++
 1 file changed, 62 insertions(+)

diff --git a/Grid/lattice/Lattice_trace.h b/Grid/lattice/Lattice_trace.h
index 9a0c6f964..db9fe14f4 100644
--- a/Grid/lattice/Lattice_trace.h
+++ b/Grid/lattice/Lattice_trace.h
@@ -124,6 +124,68 @@ Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScala
   return ret;
 }
 
+template<int N>
+Lattice<iMatrix<iScalar<iScalar<iScalar<vComplexD> > > , N> > Inverse(const Lattice<iMatrix<iScalar<iScalar<iScalar<vComplexD> > >, N> > &Umu)
+{
+  GridBase *grid=Umu.Grid();
+  auto lvol = grid->lSites();
+  Lattice<iMatrix<iScalar<iScalar<iScalar<vComplexD> > >, N > > ret(grid);
+  
+  autoView(Umu_v,Umu,CpuRead);
+  autoView(ret_v,ret,CpuWrite);
+  thread_for(site,lvol,{
+    Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
+    Coordinate lcoor;
+    grid->LocalIndexToLocalCoor(site, lcoor);
+    iMatrix<iScalar<iScalar<iScalar<ComplexD> > >, N > Us;
+    iMatrix<iScalar<iScalar<iScalar<ComplexD> > >, N > Ui;
+    peekLocalSite(Us, Umu_v, lcoor);
+    for(int i=0;i<N;i++){
+      for(int j=0;j<N;j++){
+	EigenU(i,j) = Us(i,j)()()();
+      }}
+    Eigen::MatrixXcd EigenUinv = EigenU.inverse();
+    for(int i=0;i<N;i++){
+      for(int j=0;j<N;j++){
+	Ui(i,j)()()() = EigenUinv(i,j);
+      }}
+    pokeLocalSite(Ui,ret_v,lcoor);
+  });
+  return ret;
+}
+
+template<int N>
+Lattice<iMatrix<iScalar<iScalar<iScalar<iScalar<vComplexD> > > > , N> > Inverse(const Lattice<iMatrix<iScalar<iScalar<iScalar<iScalar<vComplexD> > > >, N> > &Umu)
+{
+  GridBase *grid=Umu.Grid();
+  auto lvol = grid->lSites();
+  Lattice<iMatrix<iScalar<iScalar<iScalar<iScalar<vComplexD> > > >, N > > ret(grid);
+  
+  autoView(Umu_v,Umu,CpuRead);
+  autoView(ret_v,ret,CpuWrite);
+  thread_for(site,lvol,{
+    Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
+    Coordinate lcoor;
+    grid->LocalIndexToLocalCoor(site, lcoor);
+    iMatrix<iScalar<iScalar<iScalar<iScalar<ComplexD> > > >, N > Us;
+    iMatrix<iScalar<iScalar<iScalar<iScalar<ComplexD> > > >, N > Ui;
+    peekLocalSite(Us, Umu_v, lcoor);
+    for(int i=0;i<N;i++){
+      for(int j=0;j<N;j++){
+	EigenU(i,j) = Us(i,j)()()()();
+      }}
+    Eigen::MatrixXcd EigenUinv = EigenU.inverse();
+    for(int i=0;i<N;i++){
+      for(int j=0;j<N;j++){
+	Ui(i,j)()()()() = EigenUinv(i,j);
+      }}
+    pokeLocalSite(Ui,ret_v,lcoor);
+  });
+  return ret;
+}
+
+
+
 
 NAMESPACE_END(Grid);
 #endif