Merge pull request #15 from fmgjcoppens/replacement_updates

Replacement updates

Makefile

@@ -1,18 +1,18 @@
 ## Compilers
 ARCH = 
-H5CXX = h5c++ -std=gnu++11
+CXX = icpc
-CXX = g++
+FC = ifort
-FC = gfortran
+H5CXX = h5c++
 
 ## Compiler flags
-H5CXXFLAGS = -O0 -g
+CXXFLAGS = -O0 -debug full -traceback
-CXXFLAGS = -O0 -g
+FFLAGS = $(CXXFLAGS)
-FFLAGS = -O0 -g
+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

datasets/update_cycle_8169_repl.dat

@@ -0,0 +1,32 @@
+#START_PACKET
+#CYCLE_ID:   8169
+#SLATER_MATRIX_DIM:  4
+#NUPDATES:  2
+#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
+(01,01)    0.100000000000000E+01    0.100000000000000E+01
+(01,02)    0.000000000000000E+00    0.100000000000000E+01
+(01,03)    0.100000000000000E+01   -0.100000000000000E+01
+(01,04)   -0.100000000000000E+01   -0.100000000000000E+01
+(02,01)    0.000000000000000E+00   -0.100000000000000E+01
+(02,02)    0.100000000000000E+01    0.000000000000000E+00
+(02,03)    0.100000000000000E+01    0.100000000000000E+01
+(02,04)    0.000000000000000E+00    0.000000000000000E+00
+(03,01)   -0.100000000000000E+01    0.100000000000000E+01
+(03,02)    0.000000000000000E+00    0.200000000000000E+01
+(03,03)   -0.100000000000000E+01   -0.200000000000000E+01
+(03,04)    0.000000000000000E+00   -0.100000000000000E+01
+(04,01)    0.100000000000000E+01    0.100000000000000E+01
+(04,02)    0.100000000000000E+01    0.100000000000000E+01
+(04,03)    0.100000000000000E+01   -0.100000000000000E+01
+(04,04)    0.100000000000000E+01    0.000000000000000E+00
+#COL_UPDATE_INDEX: 2
+#COL_UPDATE_COMP_(01):   0.000000000000000E+00
+#COL_UPDATE_COMP_(02):  -1.000000000000000E+00
+#COL_UPDATE_COMP_(03):   0.000000000000000E+00
+#COL_UPDATE_COMP_(04):   1.000000000000000E+00
+#COL_UPDATE_INDEX: 4
+#COL_UPDATE_COMP_(01):  -1.000000000000000E+00
+#COL_UPDATE_COMP_(02):  -1.000000000000000E+00
+#COL_UPDATE_COMP_(03):   0.000000000000000E+00
+#COL_UPDATE_COMP_(04):  -1.000000000000000E+00
+#END_PACKET

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/SM_MaponiA3.cpp

@@ -17,7 +17,7 @@ void selectBestUpdate(unsigned int l, unsigned int N_updates,
     breakdown = abs(1 + ylk[l][index][component]);
     #ifdef DEBUG
     cout << "Inside selectBestUpdate()" << endl;
-    cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "])" << endl;
+    cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << endl;
     cout << endl;
     #endif
     if (breakdown > max) {
@@ -69,7 +69,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
   for (k = 1; k < N_updates + 1; k++) {
     #ifdef DEBUG
     cout << "Compute y0k: " << endl;
-    cout << "ylk[0][" << k << "][:]" << endl;
+    cout << "ylk[0][" << k << "][:]";
     cout << endl;
     #endif
     for (i = 1; i < Dim + 1; i++) {
@@ -78,6 +78,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
                       * Updates[(k-1)*Dim + (j-1)];
       }
     }
+    #ifdef DEBUG
+    showVector(ylk[0][k], Dim, "");
+    #endif
   }
 
   // Calculate the {y_{l,k}} from the {y_{0,k}}
@@ -96,7 +99,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     #ifdef DEBUG
     cout << "p[l+1] = " << p[l+1] << endl;
     cout << "component = " << component << endl;
-    cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl;
+    cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl;
     cout << endl;
     #endif
     if (fabs(beta) < 1e-6) {

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,10 +6,10 @@ 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", &
+    call Read_dataset("datasets/update_cycle_8169_repl.dat", &
                        cycle_id, &
                        dim, &
                        n_updates, &
@@ -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)
+    call MaponiA3(S_inv_t, dim, n_updates, U, Updates_index)
-    !! S_inv_t needs to be transposed back before we
+    !! S_inv_t needs to be transposed back before it
-    !! can multiply it with S to test unity
+    !! 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")
@@ -74,5 +77,5 @@ program QMCChem_dataset_test
     close(4000)
     close(5000)
 
-    deallocate(S, S_inv, S_inv_t, Updates, Updates_index)
+    deallocate(S, S_inv, S_inv_t, Updates, U, Updates_index)
 end program

tests/test_external_h5.cpp

@@ -47,6 +47,9 @@ int test_cycle(H5File file, int cycle) {
   double * updates = new double[nupdates*dim];
   read_double(file, group + "/updates", updates);
 
+  double * u = new double[nupdates*dim];
+
+
   /* Test */
 #ifdef DEBUG
   showMatrix(slater_matrix, dim, "OLD Slater");
@@ -59,12 +62,13 @@ int test_cycle(H5File file, int cycle) {
   for (j = 0; j < nupdates; j++) {
     for (i = 0; i < dim; i++) {
       col = col_update_index[j];
-      slater_matrix[i*dim + (col - 1)] += updates[i + j*dim];
+      u[i + j*dim] = updates[i + j*dim] - slater_matrix[i*dim + (col - 1)];
+      slater_matrix[i*dim + (col - 1)] = updates[i + j*dim];
     }
   }
 
   //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
-  SM(slater_inverse, dim, nupdates, updates, col_update_index);
+  SM(slater_inverse, dim, nupdates, u, col_update_index);
 
 #ifdef DEBUG
   showMatrix(slater_matrix, dim, "NEW Slater");
@@ -82,7 +86,7 @@ int test_cycle(H5File file, int cycle) {
   showMatrix(res, dim, "Result");
 #endif
 
-  delete [] res, updates, col_update_index,
+  delete [] res, updates, u, col_update_index,
             slater_matrix, slater_inverse;
 
   return ok;

tests/test_internal_h5.cpp

@@ -8,7 +8,7 @@
 #include "Helpers.hpp"
 
 using namespace H5;
-//#define DEBUG
+// #define DEBUG
 
 const H5std_string FILE_NAME( "datasets.hdf5" );
 
@@ -47,6 +47,8 @@ int test_cycle(H5File file, int cycle) {
   double * updates = new double[nupdates*dim];
   read_double(file, group + "/updates", updates);
 
+  double * u = new double[nupdates*dim];
+
   /* Test */
 #ifdef DEBUG
   showMatrix(slater_matrix, dim, "OLD Slater");
@@ -59,11 +61,12 @@ int test_cycle(H5File file, int cycle) {
   for (j = 0; j < nupdates; j++) {
     for (i = 0; i < dim; i++) {
       col = col_update_index[j];
-      slater_matrix[i*dim + (col - 1)] += updates[i + j*dim];
+      u[i + j*dim] = updates[i + j*dim] - slater_matrix[i*dim + (col - 1)];
+      slater_matrix[i*dim + (col - 1)] = updates[i + j*dim];
     }
   }
 
-  MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
+  MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
   //SM(slater_inverse, dim, nupdates, updates, col_update_index);
 
 #ifdef DEBUG
@@ -82,7 +85,7 @@ int test_cycle(H5File file, int cycle) {
   showMatrix(res, dim, "Result");
 #endif
 
-  delete [] res, updates, col_update_index,
+  delete [] res, updates, u, col_update_index,
             slater_matrix, slater_inverse;
 
   return ok;