From ef97b40196e86007ae9562994371174dbbdef355 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <francois.coppens@uvsq.fr>
Date: Fri, 12 Mar 2021 12:38:50 +0100
Subject: [PATCH] Changed code to treat the dataset as replacement updates,
 instead of additive updates. Before the updates are used to compute the
 inverse they are transformed accordingly. However, this breaks the 4x4
 example in update cycle 8169.

---
 Makefile                            | 16 ++++++++--------
 octave-tests/QMCChem_dataset_test.m | 20 +++++++++++++++++++-
 smvarsrc                            |  3 +++
 src/det.irp.f                       |  5 +----
 tests/QMCChem_dataset_test.f90      | 17 ++++++++++-------
 5 files changed, 41 insertions(+), 20 deletions(-)

diff --git a/Makefile b/Makefile
index c4fc58b..1249960 100644
--- a/Makefile
+++ b/Makefile
@@ -1,18 +1,18 @@
 ## Compilers
-ARCH =
-H5CXX = h5c++ -std=gnu++11
-CXX = g++
-FC = gfortran
+ARCH = 
+CXX = icpc
+FC = ifort
+H5CXX = h5c++
 
 ## Compiler flags
-H5CXXFLAGS = -O0 -g
-CXXFLAGS = -O0 -g
-FFLAGS = -O0 -g
+CXXFLAGS = -O0 -debug full -traceback
+FFLAGS = $(CXXFLAGS)
+H5CXXFLAGS = $(CXXFLAGS) -fPIC
+FLIBS = -lstdc++
 
 INCLUDE = -I $(INC_DIR)/
 DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o
 DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
-FLIBS = -lstdc++
 
 SRC_DIR := src
 TST_DIR := tests
diff --git a/octave-tests/QMCChem_dataset_test.m b/octave-tests/QMCChem_dataset_test.m
index 2861551..b43ab74 100644
--- a/octave-tests/QMCChem_dataset_test.m
+++ b/octave-tests/QMCChem_dataset_test.m
@@ -17,4 +17,22 @@ printf("--------------------------------------\n")
 printf("Determinant of S x S_inv : %f\n", det(S*S_inv))
 printf("Trace of S x S_inv       : %f\n", trace(S*S_inv))
 printf("Norm of S x S_inv        : %f\n", norm(S*S_inv))
-printf("======================================\n")
+
+cutoff = 1e-6;
+printf("\n")
+printf("Cutoff set to %e: S x S_inv = \n", cutoff)
+printf(" -----------------------------------------\n")
+dim = columns(S);
+res=S*S_inv;
+for i = 1:dim
+    for j = 1:dim
+        if (res(i,j) < cutoff)
+            res(i,j) = 0;
+        elseif (S(i,j)-1 < cutoff)
+            res(i,j) = 1;
+        endif
+    endfor
+endfor
+format free;
+disp(res);
+printf(" =========================================\n")
diff --git a/smvarsrc b/smvarsrc
index 0265f11..a4434e3 100644
--- a/smvarsrc
+++ b/smvarsrc
@@ -1,3 +1,6 @@
 #!/usr/bin/env bash
 
+export HDF5_CXX=icpc
+export HDF5_CXXLINKER=icpc
+
 export PATH=$PWD/bin:$PATH
diff --git a/src/det.irp.f b/src/det.irp.f
index 10a0968..e83d64a 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -970,7 +970,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
   !DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0)
   !DIR$ ASSUME (n>150)
 
-  integer                        :: i,j,n4
+  integer          :: i,j,n4
   double precision :: zl
 
   !DIR$ NOPREFETCH
@@ -1309,9 +1309,6 @@ END_PROVIDER
   write(501,10007)
   write(502,10007)
   cycle_id = cycle_id + 1
-  !close(501) !! Close file
-  !close(502) !! Close file
-  !stop
   
   else
 
diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90
index 508c0f3..2f537a8 100644
--- a/tests/QMCChem_dataset_test.f90
+++ b/tests/QMCChem_dataset_test.f90
@@ -6,7 +6,7 @@ program QMCChem_dataset_test
     integer :: i, j, col !! Iterators
     integer :: cycle_id, dim, n_updates
     integer(c_int), dimension(:), allocatable :: Updates_index
-    real(c_double), dimension(:,:), allocatable :: Updates
+    real(c_double), dimension(:,:), allocatable :: Updates, U
     real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t
 
     call Read_dataset("datasets/update_cycle_8169.dat", &
@@ -17,7 +17,7 @@ program QMCChem_dataset_test
                        S_inv, &
                        Updates_index, &
                        Updates)
-    allocate(S_inv_t(dim,dim))
+    allocate(S_inv_t(dim,dim), U(dim,n_updates))
 
     !! Write current S and S_inv to file for check in Octave
     open(unit = 2000, file = "Slater_old.dat")
@@ -43,11 +43,13 @@ program QMCChem_dataset_test
     end do
     close(2000)
 
-    !! Update S
+    !! Update S & transform replacement updates 'Updates'
+    !! into additive updates 'U' to compute the inverse
     do j=1,n_updates
         do i=1,dim
             col = Updates_index(j)
-            S(i,col) = S(i,col) + Updates(i,j)
+            U(i,j) = Updates(i,j) - S(i,col)
+            S(i,col) = Updates(i,j)
         end do
     end do
 
@@ -55,11 +57,12 @@ program QMCChem_dataset_test
     !! S_inv needs to be transposed first before it
     !! goes to MaponiA3
     call Transpose(S_inv, S_inv_t, dim)
-    call MaponiA3(S_inv_t, dim, n_updates, Updates, Updates_index)
-    !! S_inv_t needs to be transposed back before we
-    !! can multiply it with S to test unity
+    call MaponiA3(S_inv_t, dim, n_updates, U, Updates_index)
+    !! S_inv_t needs to be transposed back before it
+    !! can be multiplied with S to test unity
     call Transpose(S_inv_t, S_inv, dim)
 
+
     !! Write new S and S_inv to file for check in Octave
     open(unit = 4000, file = "Slater.dat")
     open(unit = 5000, file = "Slater_inv.dat")