TO BE AMENDED

Helpers_mod.f90

@@ -1,4 +1,4 @@
-module Utils
+module Helpers
     implicit none
     contains
         subroutine Read_dataset(filename, &
@@ -53,4 +53,19 @@ module Utils
             !! End of reading the dataset from file
            close(1000)
         end subroutine Read_dataset
-end module Utils
+
+        subroutine Transpose(S_inv, S_inv_transpose, dim)
+            use, intrinsic :: iso_c_binding, only : c_int, c_double
+            implicit none
+            
+            integer, intent(in) :: dim
+            real(c_double), allocatable, intent(in) :: S_inv(:,:)
+            real(c_double), allocatable, intent(inout) :: S_inv_transpose(:,:)
+            integer :: i, j
+            do i=1,dim
+                do j=1,dim
+                    S_inv_transpose(i,j) = S_inv(j,i)
+                end do
+            end do
+        end subroutine Transpose
+    end module Helpers

Makefile

@@ -2,12 +2,12 @@
 
 ## Used compilers
 H5CXX = h5c++
-CXX = clang++
-FC = flang
+CXX = icpc
+FC = ifort
 
 ## Compiler flags & common obs & libs
-CXXFLAGS = -O0 #-debug full -traceback
-FFLAGS = -O0 #-debug full -traceback
+CXXFLAGS = -O0 -debug full -traceback
+FFLAGS = -O0 -debug full -traceback
 FLIBS = -lstdc++
 OBJS = SM_MaponiA3.o
 
@@ -15,7 +15,7 @@ OBJS = SM_MaponiA3.o
 cMaponiA3_test_3x3_3OBJ = cMaponiA3_test_3x3_3.o
 fMaponiA3_test_3x3_3OBJ = SM_MaponiA3_mod.o fMaponiA3_test_3x3_3.o
 fMaponiA3_test_4x4_2OBJ = SM_MaponiA3_mod.o fMaponiA3_test_4x4_2.o
-QMCChem_dataset_testOBJ = Utils_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o
+QMCChem_dataset_testOBJ = Helpers_mod.o SM_MaponiA3_mod.o QMCChem_dataset_test.o
 
 
 ## Default build target: build everything

QMCChem_dataset_test.f90

@@ -1,16 +1,15 @@
 program QMCChem_dataset_test
-    use Sherman_Morrison, only : MaponiA3
-    use Utils, only : Read_dataset
-    use, intrinsic :: iso_c_binding, only : c_int, c_double
+    use Sherman_Morrison
+    use Helpers
     implicit none
 
     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 :: S, S_inv, S_inv_trans
+    real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t
 
-    call Read_dataset("update_cycle_13.dat", &
+    call Read_dataset("datasets/update_cycle_13.dat", &
                        cycle_id, &
                        dim, &
                        n_updates, &
@@ -18,7 +17,9 @@ program QMCChem_dataset_test
                        S_inv, &
                        Updates_index, &
                        Updates)
-    
+    allocate(S_inv_t(dim,dim))
+    call Transpose(S_inv, S_inv_t, dim)
+
     !! Write current S and S_inv to file for check in Octave
     open(unit = 2000, file = "Slater_old.dat")
     open(unit = 3000, file = "Slater_old_inv.dat")
@@ -52,7 +53,8 @@ program QMCChem_dataset_test
     end do
 
     !! Update S_inv
-    call MaponiA3(S_inv, dim, n_updates, Updates, Updates_index)
+    call MaponiA3(S_inv_t, dim, n_updates, Updates, Updates_index)
+    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")
@@ -68,5 +70,5 @@ program QMCChem_dataset_test
     close(4000)
     close(5000)
 
-    deallocate(S, S_inv, Updates, Updates_index)
+    deallocate(S, S_inv, S_inv_t, Updates, Updates_index)
 end program

SM_MaponiA3.cpp

@@ -30,7 +30,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   for (k = 1; k < N_updates + 1; k++) {
     for (i = 1; i < Dim + 1; i++) {
       for (j = 1; j < Dim + 1; j++) {
-        ylk[0][k][i] += Slater_inv[(i-1) + (j-1)*Dim]
+        ylk[0][k][i] += Slater_inv[(i-1)*Dim + (j-1)]
                       * Updates[(k-1)*Dim + (j-1)];
       }
     }
@@ -70,7 +70,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   // a new pointer 'copy' that points to whereever 'Slater_inv' points to now.
   double *copy = Slater_inv;
 
-  Slater_inv = transpose(Slater_inv, Dim);
+  // Slater_inv = transpose(Slater_inv, Dim);
 
   // Construct A-inverse from A0-inverse and the ylk
   for (l = 0; l < N_updates; l++) {       // l    = 0, 1
@@ -85,7 +85,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     Slater_inv = matMul(Al, Slater_inv, Dim);
   }
 
-  Slater_inv = transpose(Slater_inv, Dim);
+  // Slater_inv = transpose(Slater_inv, Dim);
 
   // Assign the new values of 'Slater_inv' to the old values in 'copy[][]'
   for (i = 0; i < Dim; i++) {

SM_MaponiA3_mod.f90

@@ -1,4 +1,5 @@
 module Sherman_Morrison
+    use, intrinsic :: iso_c_binding, only : c_int, c_double
     interface
         subroutine MaponiA3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_f")
             use, intrinsic :: iso_c_binding, only : c_int, c_double

tests/convert-to-h5.py

@@ -7,7 +7,7 @@ def rl(rf):
 
 
 with h5py.File('datasets.hdf5', 'w') as f:
-    with open('dataset.dat', 'r') as rf:
+    with open('datasets.dat', 'r') as rf:
         while(1):
             line = rl(rf)
             if not line or not line.startswith('#START_PACKET'):