From 983f87d504ab7da3770511903864e1dc17442aba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Mon, 1 Mar 2021 16:08:14 +0100 Subject: [PATCH 01/11] Rewrote matMul function so it doesn't declare memory inside. --- Makefile | 8 ++++---- include/Helpers.hpp | 13 +++++++++++++ src/SM_MaponiA3.cpp | 32 +++++++++++++------------------- tests/QMCChem_dataset_test.f90 | 2 +- 4 files changed, 31 insertions(+), 24 deletions(-) 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..9f3426b 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -3,6 +3,7 @@ #include #include #include +#include using namespace std; template @@ -76,6 +77,18 @@ T *matMul(T *A, T *B, unsigned int M) { return C; } +template +void matMul2(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++) { + C[i*M+j] += A[i*M+k] * B[k*M+j]; + } + } + } +} + template T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) { T1 *C = new T1[M*M]; diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index c09f40f..dfdd093 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -68,16 +68,13 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } - // 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; - + double *next = new double[Dim*Dim] {0}; + double *last; + last = 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 + for (l = 0; l < N_updates; l++) { + k = l + 1; + component = Updates_index[p[k] - 1]; beta = 1 + ylk[l][p[k]][component]; for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { @@ -85,24 +82,21 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, * ylk[l][p[k]][i + 1] / beta; } } - Slater_inv = matMul(Al, Slater_inv, Dim); - } - // 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]; - } + matMul2(Al, last, next, Dim); + double *tmp = next; + next = last; + last = tmp; } - + memcpy(Slater_inv, last, Dim*Dim*sizeof(double)); + 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[] Al, next; delete[] p, breakdown; } 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, & From 3f607797001786e215a3026c7e7eddd1faa35ea9 Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Tue, 2 Mar 2021 10:24:14 +0100 Subject: [PATCH 02/11] Cleanup --- include/Helpers.hpp | 2 +- src/SM_MaponiA3.cpp | 8 +++----- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/include/Helpers.hpp b/include/Helpers.hpp index 9f3426b..5f52bee 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -79,7 +79,7 @@ T *matMul(T *A, T *B, unsigned int M) { template void matMul2(T *A, T *B, T *C, unsigned int M) { - memset(C,0,M*M*sizeof(T)); + 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++) { diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index dfdd093..9fe3bc9 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -67,11 +67,10 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } } - - double *next = new double[Dim*Dim] {0}; - double *last; - last = Slater_inv; + // Construct A-inverse from A0-inverse and the ylk + double *last = Slater_inv; + double *next = new double[Dim*Dim] {0}; for (l = 0; l < N_updates; l++) { k = l + 1; component = Updates_index[p[k] - 1]; @@ -82,7 +81,6 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, * ylk[l][p[k]][i + 1] / beta; } } - matMul2(Al, last, next, Dim); double *tmp = next; next = last; From b2c6847ed8924c9df5b810118495ee3003de2e4e Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Tue, 2 Mar 2021 11:50:08 +0100 Subject: [PATCH 03/11] Add debug comments --- src/SM_MaponiA3.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index 9fe3bc9..e18a5a3 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -61,11 +61,14 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } for (k = l + 1; k < N_updates + 1; k++) { alpha = ylk[l - 1][p[k]][component] / beta; + cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; for (i = 1; i < Dim + 1; i++) { + cout << "ylk[" << l << "][p[" << k << "]][" << i << "] = ylk[" << l - 1 << "][p[" << k << "]][" << i << "] - alpha * ylk[" << l - 1 << "][p[" << l << "]][" << i << "]" << endl; ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - alpha * ylk[l - 1][p[l]][i]; } } + } // Construct A-inverse from A0-inverse and the ylk @@ -75,6 +78,8 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, k = l + 1; component = Updates_index[p[k] - 1]; beta = 1 + ylk[l][p[k]][component]; + cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; + cout << "ylk[" << l << "][" << p[k] << "][i + 1]" << endl; for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { Al[i*Dim + j] = (i == j) - (j == component-1) @@ -88,6 +93,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } memcpy(Slater_inv, last, Dim*Dim*sizeof(double)); + // Free memory for (l = 0; l < N_updates; l++) { for (k = 0; k < N_updates + 1; k++) { delete[] ylk[l][k]; From 9c82092cff91148eab8b6ec716d3a51d6db7cc7b Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Thu, 4 Mar 2021 11:10:59 +0100 Subject: [PATCH 04/11] Removed dependence on a breakdown array. --- include/Helpers.hpp | 14 ----------- src/SM_MaponiA3.cpp | 59 +++++++++++++++++++++++++++++++-------------- 2 files changed, 41 insertions(+), 32 deletions(-) diff --git a/include/Helpers.hpp b/include/Helpers.hpp index 5f52bee..958b9ed 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -6,20 +6,6 @@ #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; diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index e18a5a3..05fb614 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -8,12 +8,15 @@ 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, lbar, i, j, tmp, component, max, breakdown; 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]; - + // Populate update-order vector for (i = 0; i < N_updates; i++) { p[i + 1] = i + 1; @@ -27,6 +30,10 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, ylk[l][k] = new double[Dim + 1] {0}; } } + + /* + START THE ALGORITHM + */ // Calculate the y0k for (k = 1; k < N_updates + 1; k++) { @@ -38,21 +45,36 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } + // // Compute A_1_inv and A_1_inv * Slater_inv + // double *lastIntermRes = Slater_inv; + // double *Res = new double[Dim*Dim]{0}; + // component = Updates_index[0]; + // beta = 1 + ylk[0][1][component]; + // for (i = 0; i < Dim; i++) { + // for (j = 0; j < Dim; j++) { + // Al[i*Dim + j] = (i == j) - (j == component-1) * ylk[0][1][i + 1] / beta; + // } + // } + // matMul2(Al, Slater_inv, Res, Dim); + // Slater_inv = Res; + // Res = lastIntermRes; + // lastIntermRes = Slater_inv; + // showMatrix(Slater_inv, Dim, "Slater_inv"); + // Calculate all the ylk from the y0k for (l = 1; l < N_updates; l++) { + + // Select update with largest break-down val + lbar = l; max = 0; 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; + breakdown = abs(1 + ylk[l - 1][p[j]][component]); + if (breakdown > max) { + max = breakdown; + lbar = j; + } + } tmp = p[l]; p[l] = p[lbar]; p[lbar] = tmp; + component = Updates_index[p[l] - 1]; beta = 1 + ylk[l - 1][p[l]][component]; if (fabs(beta) < 1e-6) { @@ -68,7 +90,6 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, - alpha * ylk[l - 1][p[l]][i]; } } - } // Construct A-inverse from A0-inverse and the ylk @@ -93,15 +114,17 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } memcpy(Slater_inv, last, Dim*Dim*sizeof(double)); - // Free memory + /* + 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, next; - delete[] p, breakdown; + delete[] Al, next, p; } extern "C" { From 066e50914cd03aa5ded197b08fa53d1a8adc3bd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Thu, 4 Mar 2021 18:33:28 +0100 Subject: [PATCH 05/11] Regularised main loop. --- src/SM_MaponiA3.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index 05fb614..d73cf6b 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -62,32 +62,32 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, // showMatrix(Slater_inv, Dim, "Slater_inv"); // Calculate all the ylk from the y0k - for (l = 1; l < N_updates; l++) { + for (l = 0; l < N_updates; l++) { // Select update with largest break-down val - lbar = l; max = 0; - for (j = l; j < N_updates + 1; j++) { + lbar = l+1; max = 0; + for (j = l+1; j < N_updates + 1; j++) { component = Updates_index[p[j] - 1]; - breakdown = abs(1 + ylk[l - 1][p[j]][component]); + breakdown = abs(1 + ylk[l+1 - 1][p[j]][component]); if (breakdown > max) { max = breakdown; lbar = j; } - } tmp = p[l]; p[l] = p[lbar]; p[lbar] = tmp; + } tmp = p[l+1]; p[l+1] = p[lbar]; p[lbar] = tmp; - component = Updates_index[p[l] - 1]; - beta = 1 + ylk[l - 1][p[l]][component]; - if (fabs(beta) < 1e-6) { + component = Updates_index[p[l+1] - 1]; + beta = 1 + ylk[l+1 - 1][p[l+1]][component]; + if (beta == 0) { 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 (k = l+1 + 1; k < N_updates + 1; k++) { + alpha = ylk[l+1 - 1][p[k]][component] / beta; cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; for (i = 1; i < Dim + 1; i++) { cout << "ylk[" << l << "][p[" << k << "]][" << i << "] = ylk[" << l - 1 << "][p[" << k << "]][" << i << "] - alpha * ylk[" << l - 1 << "][p[" << l << "]][" << i << "]" << endl; - ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - - alpha * ylk[l - 1][p[l]][i]; + ylk[l+1][p[k]][i] = ylk[l+1 - 1][p[k]][i] + - alpha * ylk[l+1 - 1][p[l+1]][i]; } } } From a9604b3838d708e0721fcda87a164edd6689fa3c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Thu, 4 Mar 2021 18:36:51 +0100 Subject: [PATCH 06/11] Cleaned-up main loop. --- src/SM_MaponiA3.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index d73cf6b..2ab3c16 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -76,18 +76,18 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } tmp = p[l+1]; p[l+1] = p[lbar]; p[lbar] = tmp; component = Updates_index[p[l+1] - 1]; - beta = 1 + ylk[l+1 - 1][p[l+1]][component]; + beta = 1 + ylk[l][p[l+1]][component]; if (beta == 0) { cout << "Break-down occured. Exiting..." << endl; exit(1); } - for (k = l+1 + 1; k < N_updates + 1; k++) { - alpha = ylk[l+1 - 1][p[k]][component] / beta; + for (k = l+2; k < N_updates + 1; k++) { + alpha = ylk[l][p[k]][component] / beta; cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; for (i = 1; i < Dim + 1; i++) { cout << "ylk[" << l << "][p[" << k << "]][" << i << "] = ylk[" << l - 1 << "][p[" << k << "]][" << i << "] - alpha * ylk[" << l - 1 << "][p[" << l << "]][" << i << "]" << endl; - ylk[l+1][p[k]][i] = ylk[l+1 - 1][p[k]][i] - - alpha * ylk[l+1 - 1][p[l+1]][i]; + ylk[l+1][p[k]][i] = ylk[l][p[k]][i] + - alpha * ylk[l][p[l+1]][i]; } } } From 3b3567040b8627bcd44e6fef67a827efea91bf84 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 10 Mar 2021 11:41:34 +0100 Subject: [PATCH 07/11] Added more debug output. --- Makefile | 2 +- src/SM_MaponiA3.cpp | 79 +++++++++++++++++++++++++++------------------ 2 files changed, 48 insertions(+), 33 deletions(-) diff --git a/Makefile b/Makefile index c4fc58b..3afb9d3 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ FC = gfortran ## Compiler flags H5CXXFLAGS = -O0 -g -CXXFLAGS = -O0 -g +CXXFLAGS = -O0 -g -DDEBUG FFLAGS = -O0 -g INCLUDE = -I $(INC_DIR)/ diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index 2ab3c16..ce55ff1 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -4,6 +4,34 @@ #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; + unsigned int max = 0; + unsigned int component = 0; + unsigned int index = 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) { @@ -12,7 +40,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, DECLARE AND INITIALISE ARRAYS */ - unsigned int k, l, lbar, i, j, tmp, component, max, breakdown; + unsigned int k, l, i, j, component; unsigned int *p = new unsigned int[N_updates + 1] {0}; double alpha, beta; double *Al = new double[Dim * Dim]; @@ -45,47 +73,29 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } - // // Compute A_1_inv and A_1_inv * Slater_inv - // double *lastIntermRes = Slater_inv; - // double *Res = new double[Dim*Dim]{0}; - // component = Updates_index[0]; - // beta = 1 + ylk[0][1][component]; - // for (i = 0; i < Dim; i++) { - // for (j = 0; j < Dim; j++) { - // Al[i*Dim + j] = (i == j) - (j == component-1) * ylk[0][1][i + 1] / beta; - // } - // } - // matMul2(Al, Slater_inv, Res, Dim); - // Slater_inv = Res; - // Res = lastIntermRes; - // lastIntermRes = Slater_inv; - // showMatrix(Slater_inv, Dim, "Slater_inv"); - // Calculate all the ylk from the y0k for (l = 0; l < N_updates; l++) { - // Select update with largest break-down val - lbar = l+1; max = 0; - for (j = l+1; j < N_updates + 1; j++) { - component = Updates_index[p[j] - 1]; - breakdown = abs(1 + ylk[l+1 - 1][p[j]][component]); - if (breakdown > max) { - max = breakdown; - lbar = j; - } - } tmp = p[l+1]; p[l+1] = p[lbar]; p[lbar] = tmp; + // Select lk-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]; - if (beta == 0) { +#ifdef DEBUG + cout << "In outer compute-ylk-loop:" << 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+2; k < N_updates + 1; k++) { + + // Compute ylk + for (k = l+2; k < N_updates+1; k++) { alpha = ylk[l][p[k]][component] / beta; - cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; for (i = 1; i < Dim + 1; i++) { - cout << "ylk[" << l << "][p[" << k << "]][" << i << "] = ylk[" << l - 1 << "][p[" << k << "]][" << i << "] - alpha * ylk[" << l - 1 << "][p[" << l << "]][" << i << "]" << endl; ylk[l+1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l+1]][i]; } @@ -99,8 +109,13 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, k = l + 1; component = Updates_index[p[k] - 1]; beta = 1 + ylk[l][p[k]][component]; - cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; +#ifdef DEBUG + cout << "Compute inverse. Inside l-loop: l = " << l << endl; + cout << "component = Updates_index[p[" << k << "] - 1] = Updates_index[" << p[k] - 1 << "] = " << Updates_index[p[k] - 1] << endl; + cout << "beta = 1 + ylk[" << l << "][" << p[k] << "][" << component << "]" << endl; cout << "ylk[" << l << "][" << p[k] << "][i + 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) From de8e89df75d92494fb6891a9432189011eb2abb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 10 Mar 2021 11:54:43 +0100 Subject: [PATCH 08/11] Change in debug output. --- src/SM_MaponiA3.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index ce55ff1..4226713 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -95,6 +95,11 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, // Compute ylk for (k = l+2; k < N_updates+1; k++) { alpha = ylk[l][p[k]][component] / beta; +#ifdef DEBUG + cout << "Inside k-loop of ylk-loop:" << 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]; @@ -113,7 +118,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, cout << "Compute inverse. Inside l-loop: l = " << l << endl; cout << "component = Updates_index[p[" << k << "] - 1] = Updates_index[" << p[k] - 1 << "] = " << Updates_index[p[k] - 1] << endl; cout << "beta = 1 + ylk[" << l << "][" << p[k] << "][" << component << "]" << endl; - cout << "ylk[" << l << "][" << p[k] << "][i + 1]" << endl; + cout << "ylk[" << l << "][" << p[k] << "][:]" << endl; cout << endl; #endif for (i = 0; i < Dim; i++) { From ffbde8f88c5b44c54b688b832c302bd4ec07beb9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 10 Mar 2021 12:18:08 +0100 Subject: [PATCH 09/11] More changes to debug output. --- src/SM_MaponiA3.cpp | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index 4226713..a989423 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -17,11 +17,11 @@ void selectBestUpdate(unsigned int l, unsigned int N_updates, index = p[j]; component = Updates_index[index - 1]; breakdown = abs(1 + ylk[l][index][component]); -#ifdef DEBUG + #ifdef DEBUG cout << "Inside selectBestUpdate()" << endl; cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "])" << endl; cout << endl; -#endif + #endif if (breakdown > max) { max = breakdown; lbar = j; @@ -65,6 +65,11 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, // Calculate the y0k 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)] @@ -82,11 +87,11 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, // Select component and comp. bd-condition. component = Updates_index[p[l+1] - 1]; beta = 1 + ylk[l][p[l+1]][component]; -#ifdef DEBUG + #ifdef DEBUG cout << "In outer compute-ylk-loop:" << endl; cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl; cout << endl; -#endif + #endif if (fabs(beta) < 1e-6) { cout << "Break-down occured. Exiting..." << endl; exit(1); @@ -95,11 +100,11 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, // Compute ylk for (k = l+2; k < N_updates+1; k++) { alpha = ylk[l][p[k]][component] / beta; -#ifdef DEBUG + #ifdef DEBUG cout << "Inside k-loop of ylk-loop:" << endl; cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << endl; cout << endl; -#endif + #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]; @@ -114,13 +119,13 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, k = l + 1; component = Updates_index[p[k] - 1]; beta = 1 + ylk[l][p[k]][component]; -#ifdef DEBUG + #ifdef DEBUG cout << "Compute inverse. Inside l-loop: l = " << l << endl; cout << "component = Updates_index[p[" << k << "] - 1] = Updates_index[" << p[k] - 1 << "] = " << Updates_index[p[k] - 1] << endl; cout << "beta = 1 + ylk[" << l << "][" << p[k] << "][" << component << "]" << endl; - cout << "ylk[" << l << "][" << p[k] << "][:]" << endl; + cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[k] << "][:]" << endl; cout << endl; -#endif + #endif for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { Al[i*Dim + j] = (i == j) - (j == component-1) From e18e80ff5c315d8623b44b0bcb7cff27a4814858 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 10 Mar 2021 15:27:53 +0100 Subject: [PATCH 10/11] Removed matMul() from Helpers.hpp and renamed matMul2() to matMul(). --- include/Helpers.hpp | 16 ++-------------- src/SM_MaponiA3.cpp | 2 +- 2 files changed, 3 insertions(+), 15 deletions(-) diff --git a/include/Helpers.hpp b/include/Helpers.hpp index 958b9ed..06ea441 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -37,6 +37,7 @@ void showMatrix(T *matrix, unsigned int M, string name) { cout << " ]," << endl; } cout << "]" << endl; + cout << endl; } template @@ -51,20 +52,7 @@ T *transpose(T *A, unsigned int M) { } template -T *matMul(T *A, T *B, unsigned int M) { - T *C = new T[M*M] {0}; - for (unsigned int i = 0; i < M; i++) { - for (unsigned int j = 0; j < M; j++) { - for (unsigned int k = 0; k < M; k++) { - C[i*M+j] += A[i*M+k] * B[k*M+j]; - } - } - } - return C; -} - -template -void matMul2(T *A, T *B, T *C, unsigned int M) { +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++) { diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index a989423..8aca638 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -132,7 +132,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, * ylk[l][p[k]][i + 1] / beta; } } - matMul2(Al, last, next, Dim); + matMul(Al, last, next, Dim); double *tmp = next; next = last; last = tmp; From 081fbfc1d8cf11f7c9943d7a68e9ac2293fd9701 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 10 Mar 2021 16:56:22 +0100 Subject: [PATCH 11/11] Loop merge successful. Next step: reduce the number of indices of ylk by not keeping all the previous solutions. --- Makefile | 2 +- src/SM_MaponiA3.cpp | 87 ++++++++++++++++++++++++--------------------- 2 files changed, 47 insertions(+), 42 deletions(-) diff --git a/Makefile b/Makefile index 3afb9d3..c4fc58b 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ FC = gfortran ## Compiler flags H5CXXFLAGS = -O0 -g -CXXFLAGS = -O0 -g -DDEBUG +CXXFLAGS = -O0 -g FFLAGS = -O0 -g INCLUDE = -I $(INC_DIR)/ diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index 8aca638..ed64470 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -7,10 +7,8 @@ void selectBestUpdate(unsigned int l, unsigned int N_updates, unsigned int *Updates_index, unsigned int *p, double ***ylk) { - unsigned int lbar = l+1; - unsigned int max = 0; - unsigned int component = 0; - unsigned int index = 0; + 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++) { @@ -44,13 +42,15 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int *p = new unsigned int[N_updates + 1] {0}; double alpha, beta; 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]; @@ -58,12 +58,14 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, ylk[l][k] = new double[Dim + 1] {0}; } } - + + + /* - START THE ALGORITHM + START ALGORITHM */ - // Calculate the y0k + // Calculate the {y_{0,k}} for (k = 1; k < N_updates + 1; k++) { #ifdef DEBUG cout << "Compute y0k: " << endl; @@ -78,17 +80,22 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } - // Calculate all the ylk from the y0k + // 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 - // Select lk-update with largest break-down val + // 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 << "In outer compute-ylk-loop:" << endl; + 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 @@ -97,11 +104,33 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, exit(1); } - // Compute ylk + // 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[l+1]][i + 1] / beta; + } + } + matMul(Al, last, next, Dim); + tmp = next; + next = last; + last = tmp; + #ifdef DEBUG + showMatrix(last, Dim, "last"); + #endif + + // 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 of ylk-loop:" << endl; + cout << "Inside k-loop: k = " << k << endl; cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << endl; cout << endl; #endif @@ -111,34 +140,10 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } } } - - // Construct A-inverse from A0-inverse and the ylk - double *last = Slater_inv; - double *next = new double[Dim*Dim] {0}; - for (l = 0; l < N_updates; l++) { - k = l + 1; - component = Updates_index[p[k] - 1]; - beta = 1 + ylk[l][p[k]][component]; - #ifdef DEBUG - cout << "Compute inverse. Inside l-loop: l = " << l << endl; - cout << "component = Updates_index[p[" << k << "] - 1] = Updates_index[" << p[k] - 1 << "] = " << Updates_index[p[k] - 1] << endl; - cout << "beta = 1 + ylk[" << l << "][" << p[k] << "][" << component << "]" << endl; - cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[k] << "][:]" << 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; - } - } - matMul(Al, last, next, Dim); - double *tmp = next; - next = last; - last = tmp; - } memcpy(Slater_inv, last, Dim*Dim*sizeof(double)); - + + + /* CLEANUP MEMORY */