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

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

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

smvarsrc

@@ -1,3 +1,6 @@
 #!/usr/bin/env bash
 
+export HDF5_CXX=icpc
+export HDF5_CXXLINKER=icpc
+
 export PATH=$PWD/bin:$PATH

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
 

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