Add test infrastructure for datasets

Helpers.hpp

@@ -119,3 +119,15 @@ T matDet(T **A, unsigned int M) {
     }
     delete [] temp;
 }
+
+
+template<typename T>
+bool is_identity(T *A, unsigned int M, double tolerance) {
+   for (unsigned int i = 0; i < M; i++) {
+        for (unsigned int j = 0; j < M; j++) {
+          if (i==j && fabs(A[i*M+j]-1) > tolerance) return false;
+          if (i!=j && fabs(A[i*M+j]) > tolerance) return false;
+        }
+   }
+   return true;
+}

Makefile

@@ -1,11 +1,11 @@
 ## Used compilers
-CXX = icpc
-FC = ifort
+H5CXX = h5c++
+CXX = clang++
+FC = flang
 
 ## Compiler flags
-CXXFLAGS = -O0 #-debug full -traceback
-FFLAGS = -O0 #-debug full -traceback
-# ARCH = -xCORE-AVX2
+CXXFLAGS = -O0
+FFLAGS = -O0
 
 ## Deps & objs for C++ cMaponiA3_test
 cMaponiA3_testDEP = cMaponiA3_test.cpp SM_MaponiA3.cpp SM_MaponiA3.hpp Helpers.hpp
@@ -32,7 +32,7 @@ QMCChem_dataset_testLIB = -lstdc++
 ## Build tagets
 .PHONY: all clean distclean
 
-all: cMaponiA3_test fMaponiA3_test QMCChem_dataset_test
+all: cMaponiA3_test fMaponiA3_test QMCChem_dataset_test tests/test
 
 clean:
 	@rm -vf *.o *.mod
@@ -51,3 +51,6 @@ fMaponiA3_test: $(fMaponiA3_testOBJ)
 ## Linking Fortran example program calling the C++ function 'Sherman_Morrison()'
 QMCChem_dataset_test: $(QMCChem_dataset_testOBJ)
 	$(FC) $(ARCH) $(FFLAGS) $(QMCChem_dataset_testLIB) -o $@ $^
+
+tests/test: tests/test.cpp SM_MaponiA3.o
+	$(H5CXX) $(ARCH) $(CXXFLAGS) -o $@ $^

tests/convert-to-h5.py

@@ -0,0 +1,46 @@
+import h5py
+import numpy as np
+from parse import parse
+
+def rl(rf):
+    return " ".join(rf.readline().split())
+
+
+with h5py.File('datasets.hdf5', 'w') as f:
+    with open('dataset.dat', 'r') as rf:
+        while(1):
+            line = rl(rf)
+            if not line or not line.startswith('#START_PACKET'):
+                break
+            cycle_id = parse('#CYCLE_ID: {:d}', rl(rf))[0]
+            slater_matrix_dim = parse('#SLATER_MATRIX_DIM: {:d}', rl(rf))[0]
+            nupdates = parse('#NUPDATES: {:d}', rl(rf))[0]
+            assert(rf.readline().startswith('#SLATER_MATRIX'))
+
+            # Read matrices
+            slater_matrix = np.zeros((slater_matrix_dim,slater_matrix_dim))
+            slater_inverse = np.zeros((slater_matrix_dim,slater_matrix_dim))
+            for i in range(slater_matrix_dim*slater_matrix_dim):
+                res = parse('({i:d},{j:d}) {sla:e} {inv:e}', rl(rf))
+                slater_matrix[res['i']-1, res['j']-1] = res['sla']
+                slater_inverse[res['i']-1, res['j']-1] = res['inv']
+
+            # Read updates
+            col_update_index = np.zeros(nupdates, dtype='i')
+            updates = np.zeros((nupdates, slater_matrix_dim))
+            for n in range(nupdates):
+                col_update_index[n] = parse('#COL_UPDATE_INDEX: {:d}', rl(rf))[0]
+                for i in range(slater_matrix_dim):
+                    res = parse('#COL_UPDATE_COMP_({i:d}): {x:e}', rl(rf))
+                    updates[n][res['i']-1] = res['x']
+
+            assert(rf.readline().startswith('#END_PACKET'))
+            rf.readline()
+
+            cycle = f.create_group('cycle_{}'.format(cycle_id))
+            cycle.create_dataset("slater_matrix_dim", data=slater_matrix_dim)
+            cycle.create_dataset("nupdates", data=nupdates)
+            cycle.create_dataset("slater_matrix", data=slater_matrix, compression='gzip')
+            cycle.create_dataset("slater_inverse", data=slater_inverse, compression='gzip')
+            cycle.create_dataset("col_update_index", data=col_update_index)
+            cycle.create_dataset("updates", data=updates, compression='gzip')

tests/test.cpp

@@ -0,0 +1,89 @@
+#include <iostream>
+#include <string>
+#include "hdf5/serial/hdf5.h"
+#include "H5Cpp.h"
+
+#include "../SM_MaponiA3.hpp"
+#include "../Helpers.hpp"
+
+using namespace H5;
+#define DEBUG 1
+
+const H5std_string FILE_NAME( "datasets.hdf5" );
+
+void read_int(H5File file, std::string key, unsigned int * data) {
+  DataSet ds = file.openDataSet(key);
+  ds.read(data, PredType::STD_U32LE);
+  ds.close();
+}
+
+void read_double(H5File file, std::string key, double * data) {
+  DataSet ds = file.openDataSet(key);
+  ds.read(data, PredType::IEEE_F64LE);
+  ds.close();
+}
+
+
+int test_cycle(H5File file, int cycle) {
+
+  /* Read the data */
+
+  std::string group = "cycle_" + std::to_string(cycle);
+
+  unsigned int dim, nupdates;
+  read_int(file, group + "/slater_matrix_dim", &dim);
+  read_int(file, group + "/nupdates", &nupdates);
+
+  double * slater_matrix = new double[dim*dim];
+  read_double(file, group + "/slater_matrix", slater_matrix);
+
+  double * slater_inverse = new double[dim*dim];
+  read_double(file, group + "/slater_inverse", slater_inverse);
+
+  unsigned int * col_update_index = new unsigned int[nupdates];
+  read_int(file, group + "/col_update_index", col_update_index);
+
+  double * updates = new double[nupdates*dim];
+  read_double(file, group + "/updates", updates);
+
+  /* Test */
+#ifdef DEBUG
+  showMatrix(slater_matrix, dim, "Slater");
+#endif
+
+  MaponiA3(slater_matrix, slater_inverse, dim, nupdates, updates, col_update_index);
+
+#ifdef DEBUG
+  showMatrix(slater_inverse, dim, "Inverse");
+#endif
+
+  double * res = matMul(slater_matrix, slater_inverse, dim);
+  bool ok = is_identity(res, dim, 1.0e-8);
+
+#ifdef DEBUG
+  showMatrix(res, dim, "Result");
+#endif
+
+  delete [] res, updates, col_update_index, slater_matrix, slater_inverse;
+
+  return ok;
+}
+
+int main(int argc, char **argv) {
+  if (argc != 2) {
+    std::cerr << "usage: ./test_dataset <cycle>" << std::endl;
+    return 1;
+  }
+  int cycle = std::stoi(argv[1]);
+  H5File file(FILE_NAME, H5F_ACC_RDONLY);
+
+  bool ok = test_cycle(file, 21);
+
+  if (ok) {
+    std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;
+  } else {
+    std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl;
+  }
+  return ok;
+}
+