Problem solved. Problems was caused by reading the inverse matrix S_inv the same way as the matrix S, but apparently what has been stored is transpose(S_inv) and not S_inv. Therefore it needs to be read transposed compared to S in order to for S*S_inv to give the Identity matrix.

TODO: Calls from C/C++ are still broken because inside the MaponiA3 algo there is compensation for the fact that Fortran sends it's arrays in column-major order. This must be resolved by  transposeing S_inv before it is passed to MaponiA3 and remove the compensations made there.

.gitignore

@@ -1,8 +1,8 @@
 *.o
 *.mod
 .vscode
-cMaponiA3_test
+cMaponiA3_test*
-fMaponiA3_test
+fMaponiA3_test*
 QMCChem_dataset_test
 Slater*
 Updates*

Makefile

@@ -46,7 +46,10 @@ clean:
 	@rm -vf *.o *.mod
 
 distclean: clean
-	@rm -vf cMaponiA3_test_3x3_3 fMaponiA3_test_3x3_3 QMCChem_dataset_test Slater_* Updates.dat
+	@rm -vf cMaponiA3_test_3x3_3 \
+	fMaponiA3_test_3x3_3 fMaponiA3_test_4x4_2 \
+	QMCChem_dataset_test \
+	Slater* Updates.dat
 
 
 ## Linking the C++ example program

QMCChem_dataset_test.f90

@@ -7,7 +7,8 @@ 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 :: S, S_inv, S_inv_trans, Updates
+    real(c_double), dimension(:,:), allocatable :: Updates
+    real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_trans
 
     call Read_dataset("test.dataset.dat", &
                        cycle_id, &
@@ -18,20 +19,13 @@ program QMCChem_dataset_test
                        Updates_index, &
                        Updates)
     
-    allocate(S_inv_trans(dim,dim))
-    do i=1,dim
-        do j=1,dim
-            S_inv_trans(i,j) = S_inv(j,i)
-        end do
-    end do
-
     !! 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")
     do i=1,dim
         do j=1,dim
             write(2000,"(E23.15, 1X)", advance="no") S(i,j)
-            write(3000,"(E23.15, 1X)", advance="no") S_inv_trans(i,j)
+            write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j)
         end do
         write(2000,*)
         write(3000,*)
@@ -57,8 +51,8 @@ program QMCChem_dataset_test
         end do
     end do
 
-    !! Send S_inv and Updates to MaponiA3 algo
+    !! Update S_inv
-    call MaponiA3(S_inv_trans, dim, n_updates, Updates, Updates_index)
+    call MaponiA3(S_inv, dim, n_updates, Updates, Updates_index)
 
     !! Write new S and S_inv to file for check in Octave
     open(unit = 4000, file = "Slater.dat")

QMCChem_dataset_test.m

@@ -3,10 +3,18 @@ Sold_inv = dlmread('Slater_old_inv.dat');
 S = dlmread('Slater.dat');
 S_inv = dlmread('Slater_inv.dat');
 
-det(Sold*Sold_inv)
+printf("\n")
-trace(Sold*Sold_inv)
+printf("======================================\n")
-norm(Sold*Sold_inv)
+printf("OLD Slater-matrix and inverse\n")
+printf("--------------------------------------\n")
+printf("Determinant of S x S_inv : %f\n", det(Sold*Sold_inv))
+printf("Trace of S x S_inv       : %f\n", trace(Sold*Sold_inv))
+printf("Norm of S x S_inv        : %f\n", norm(Sold*Sold_inv))
 
-det(S*S_inv)
+printf("\n")
-trace(S*S_inv)
+printf("NEW Slater-matrix and inverse\n")
-norm(S*S_inv)
+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")

Utils_mod.f90

@@ -1,8 +1,14 @@
 module Utils
     implicit none
     contains
-        subroutine Read_dataset(filename, cycle_id, dim, n_updates, &
+        subroutine Read_dataset(filename, &
-                    S, S_inv, Updates_index, Updates)
+                                cycle_id, &
+                                dim, &
+                                n_updates, &
+                                S, &
+                                S_inv, &
+                                Updates_index, &
+                                Updates)
             use, intrinsic :: iso_c_binding, only : c_int, c_double
             implicit none
             
@@ -14,34 +20,37 @@ module Utils
             integer :: i, j
             character (len = 32) :: ignore
 
-            !! Start of reading the dataset from file
+            !! Start reading dataset from file
             open(unit = 1000, file = filename)
             read(1000,*)
             read(1000,*) ignore, cycle_id
             read(1000,*) ignore, dim
             read(1000,*) ignore,n_updates
 
-            allocate(Updates_index(n_updates), S(dim,dim), &
+            !! Initialise the arrays
-                        S_inv(dim,dim), Updates(dim,n_updates))
+            allocate(Updates_index(n_updates), &
+                     S(dim,dim), &
+                     S_inv(dim,dim), &
+                     Updates(dim,n_updates))
             
             !! Read S and S_inv
             read(1000,*)
             do i=1,dim
                 do j=1,dim
-                    read(1000,*) ignore, ignore, S(i,j), S_inv(i,j)
+                    !! For some reason S_inv needs to be read transposed &
+                    !  for S*S_inv to give the identity matrix
+                    read(1000,*) ignore, ignore, S(i,j), S_inv(j,i)
                 end do
             end do
 
-            !! Read the updates Updates and Updates_index
+            !! Read Updates and Updates_index
             do j=1,n_updates
                 read(1000,*) ignore, Updates_index(j)
                 do i=1,dim
                     read(1000,*) ignore, Updates(i,j)
                 end do
             end do
-
-            read(1000,*) ignore
-            close(1000)
             !! End of reading the dataset from file
+           close(1000)
         end subroutine Read_dataset
 end module Utils