diff --git a/Makefile b/Makefile index 12deb07..c4fc58b 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,8 @@ ## Compilers ARCH = -H5CXX = h5c++ -CXX = clang++ -FC = flang +H5CXX = h5c++ -std=gnu++11 +CXX = g++ +FC = gfortran ## Compiler flags H5CXXFLAGS = -O0 -g @@ -55,7 +55,7 @@ $(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR) ## Fortran modules $(OBJ_DIR)/%_mod.o: $(SRC_DIR)/%_mod.f90 | $(OBJ_DIR) - $(FC) $(FFLAGS) $(ARCH) -module $(OBJ_DIR)/ -c -o $@ $< + $(FC) $(FFLAGS) $(ARCH) -J$(OBJ_DIR)/ -c -o $@ $< ## Fortran objects $(OBJ_DIR)/%.o: $(TST_DIR)/%.f90 | $(OBJ_DIR) diff --git a/include/Helpers.hpp b/include/Helpers.hpp index 466478d..06ea441 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -3,22 +3,9 @@ #include #include #include +#include using namespace std; -template -unsigned int getMaxIndex(T *vector, unsigned int size) { - unsigned int i = 0; - unsigned int maxi = i; - unsigned int max = vector[maxi]; - for (i = 1; i < size; i++) { - if (vector[i] > max) { - max = vector[i]; - maxi = i; - } - } - return maxi; -} - template void showScalar(T scalar, string name) { cout << name << " = " << scalar << endl << endl; @@ -50,6 +37,7 @@ void showMatrix(T *matrix, unsigned int M, string name) { cout << " ]," << endl; } cout << "]" << endl; + cout << endl; } template @@ -64,8 +52,8 @@ T *transpose(T *A, unsigned int M) { } template -T *matMul(T *A, T *B, unsigned int M) { - T *C = new T[M*M] {0}; +void matMul(T *A, T *B, T *C, unsigned int M) { + memset(C, 0, M*M*sizeof(T)); for (unsigned int i = 0; i < M; i++) { for (unsigned int j = 0; j < M; j++) { for (unsigned int k = 0; k < M; k++) { @@ -73,7 +61,6 @@ T *matMul(T *A, T *B, unsigned int M) { } } } - return C; } template diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index c09f40f..ed64470 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -4,22 +4,53 @@ #include "SM_MaponiA3.hpp" #include "Helpers.hpp" +void selectBestUpdate(unsigned int l, unsigned int N_updates, + unsigned int *Updates_index, unsigned int *p, + double ***ylk) { + unsigned int lbar = l+1, max =0; + unsigned int index = 0, component = 0; + unsigned int tmp = 0; + double breakdown = 0; + for (unsigned int j = lbar; j < N_updates + 1; j++) { + index = p[j]; + component = Updates_index[index - 1]; + breakdown = abs(1 + ylk[l][index][component]); + #ifdef DEBUG + cout << "Inside selectBestUpdate()" << endl; + cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "])" << endl; + cout << endl; + #endif + if (breakdown > max) { + max = breakdown; + lbar = j; + } + } + tmp = p[l+1]; + p[l+1] = p[lbar]; + p[lbar] = tmp; +} + void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { - unsigned int k, l, lbar, i, j, tmp, component; + /* + DECLARE AND INITIALISE ARRAYS + */ + + unsigned int k, l, i, j, component; unsigned int *p = new unsigned int[N_updates + 1] {0}; double alpha, beta; - double *breakdown = new double[N_updates + 1] {0}; double *Al = new double[Dim * Dim]; + double *next = new double[Dim*Dim] {0}; + double *last = Slater_inv, *tmp; // Populate update-order vector for (i = 0; i < N_updates; i++) { p[i + 1] = i + 1; } - // Declare auxiliary solution matrix ylk + // Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim] double ***ylk = new double **[N_updates]; for (l = 0; l < N_updates; l++) { ylk[l] = new double *[N_updates + 1]; @@ -28,8 +59,19 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } - // Calculate the y0k + + + /* + START ALGORITHM + */ + + // Calculate the {y_{0,k}} for (k = 1; k < N_updates + 1; k++) { + #ifdef DEBUG + cout << "Compute y0k: " << endl; + cout << "ylk[0][" << k << "][:]" << endl; + cout << endl; + #endif for (i = 1; i < Dim + 1; i++) { for (j = 1; j < Dim + 1; j++) { ylk[0][k][i] += Slater_inv[(i-1)*Dim + (j-1)] @@ -38,72 +80,81 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } - // Calculate all the ylk from the y0k - for (l = 1; l < N_updates; l++) { - for (j = l; j < N_updates + 1; j++) { - component = Updates_index[p[j] - 1]; - breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]); - } - lbar = getMaxIndex(breakdown, N_updates + 1); - // Reset breakdown back to 0 for next round to avoid case where - // its first element is always the largest - for (i = 0; i < N_updates + 1; i++) { - breakdown[i] = 0; - } - tmp = p[l]; - p[l] = p[lbar]; - p[lbar] = tmp; - component = Updates_index[p[l] - 1]; - beta = 1 + ylk[l - 1][p[l]][component]; + // Calculate the {y_{l,k}} from the {y_{0,k}} + for (l = 0; l < N_updates; l++) { + #ifdef DEBUG + cout << "In outer compute-ylk-loop: l = " << l << endl; + cout << endl; + #endif + + // For given l select intermediate update with largest break-down val + selectBestUpdate(l, N_updates, Updates_index, p, ylk); + + // Select component and comp. bd-condition. + component = Updates_index[p[l+1] - 1]; + beta = 1 + ylk[l][p[l+1]][component]; + #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 << endl; + #endif if (fabs(beta) < 1e-6) { cout << "Break-down occured. Exiting..." << endl; exit(1); } - for (k = l + 1; k < N_updates + 1; k++) { - alpha = ylk[l - 1][p[k]][component] / beta; - for (i = 1; i < Dim + 1; i++) { - ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - - alpha * ylk[l - 1][p[l]][i]; - } - } - } - // Keep the memory location of the passed array 'Slater_inv' before - // 'Slater_inv' gets reassigned by 'matMul(...)' in the next line, - // by creating a new pointer 'copy' that points to whereever - // 'Slater_inv' points to now. - double *copy = Slater_inv; - - // Construct A-inverse from A0-inverse and the ylk - for (l = 0; l < N_updates; l++) { // l = 0, 1 - k = l + 1; // k = 1, 2 - component = Updates_index[p[k] - 1]; // comp = 2, 4 - beta = 1 + ylk[l][p[k]][component]; + // Compute intermediate update to Slater_inv + #ifdef DEBUG + cout << "Compute intermediate update to Slater_inv" << endl; + cout << "component = " << component << endl; + cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl; + cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << endl; + cout << endl; + #endif for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { Al[i*Dim + j] = (i == j) - (j == component-1) - * ylk[l][p[k]][i + 1] / beta; + * ylk[l][p[l+1]][i + 1] / beta; } } - Slater_inv = matMul(Al, Slater_inv, Dim); - } + matMul(Al, last, next, Dim); + tmp = next; + next = last; + last = tmp; + #ifdef DEBUG + showMatrix(last, Dim, "last"); + #endif - // Assign the new values of 'Slater_inv' to the old values - // in 'copy[][]' - for (i = 0; i < Dim; i++) { - for (j = 0; j < Dim; j++) { - copy[i * Dim + j] = Slater_inv[i * Dim + j]; + // For given l != 0 compute the next {y_{l,k}} + for (k = l+2; k < N_updates+1; k++) { + alpha = ylk[l][p[k]][component] / beta; + #ifdef DEBUG + cout << "Inside k-loop: k = " << k << endl; + cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << endl; + cout << endl; + #endif + for (i = 1; i < Dim + 1; i++) { + ylk[l+1][p[k]][i] = ylk[l][p[k]][i] + - alpha * ylk[l][p[l+1]][i]; + } } } + memcpy(Slater_inv, last, Dim*Dim*sizeof(double)); + + + /* + CLEANUP MEMORY + */ + for (l = 0; l < N_updates; l++) { for (k = 0; k < N_updates + 1; k++) { delete[] ylk[l][k]; } delete[] ylk[l]; } - delete[] Al; - delete[] p, breakdown; + delete[] Al, next, p; } extern "C" { diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90 index e4c4986..508c0f3 100644 --- a/tests/QMCChem_dataset_test.f90 +++ b/tests/QMCChem_dataset_test.f90 @@ -9,7 +9,7 @@ program QMCChem_dataset_test real(c_double), dimension(:,:), allocatable :: Updates real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t - call Read_dataset("datasets/update_cycle_13.dat", & + call Read_dataset("datasets/update_cycle_8169.dat", & cycle_id, & dim, & n_updates, &