From 837e160f170b4f89b10fae27817b4652c175a78a Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Mon, 21 Jun 2021 14:29:24 +0200 Subject: [PATCH] - Added various Sherman-Morrison-Woodbury kernels - Separated debug information into 2 different debug levels. --- Makefile | 2 +- include/Helpers.hpp | 36 +++++++ include/SMWB.hpp | 17 +++- src/Helpers.cpp | 6 +- src/SMWB.cpp | 215 +++++++++++++++++++++++++++++++++++++++++- src/SM_Maponi.cpp | 28 +++--- src/SM_Standard.cpp | 8 ++ src/Woodbury.cpp | 38 +++++--- tests/run_test.sh | 2 +- tests/test_h5.cpp | 19 ++-- tests/vfc_test_h5.cpp | 10 +- 11 files changed, 332 insertions(+), 49 deletions(-) diff --git a/Makefile b/Makefile index cb7ad9b..5579147 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ ifeq ($(ENV),INTEL) CXX = icpc FC = ifort ARCH = -xCORE-AVX2 - OPT = -O0 + OPT = -O2 DEBUG = -debug full else ifeq ($(ENV),LLVM) CXX = clang++ diff --git a/include/Helpers.hpp b/include/Helpers.hpp index f2a4c5d..412263e 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -7,6 +7,7 @@ #ifdef MKL #include #endif + // #define DEBUG #ifndef THRESHOLD @@ -275,3 +276,38 @@ template T residual2(T * A, unsigned int Dim) { } return res; } + +// Computes the condition number of A using a previously computed B=A^{-1} +template T condition1(T *A, T *B, unsigned int Dim) { + T resA = 0; + T resB = 0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + T deltaA = A[i * Dim + j]; + T deltaB = B[i * Dim + j]; + resA += deltaA * deltaA; + resB += deltaB * deltaB; + } + } + return sqrt(resA * resB); +} + +#ifdef MKL +// Computes the condition number of A by first inverting it with LAPACKE +template T condition2(T *A, unsigned int Dim) { + T B[Dim * Dim]; + std::memcpy(B, A, Dim * Dim * sizeof(T)); + inverse(B, Dim); + T resA = 0; + T resB = 0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + T deltaA = A[i * Dim + j]; + T deltaB = B[i * Dim + j]; + resA += deltaA * deltaA; + resB += deltaB * deltaB; + } + } + return sqrt(resA) * sqrt(resB); +} +#endif diff --git a/include/SMWB.hpp b/include/SMWB.hpp index 0de1a49..28adda4 100644 --- a/include/SMWB.hpp +++ b/include/SMWB.hpp @@ -1,4 +1,19 @@ // Sherman-Morrison-Woodbury kernel 1 -// WB2, WB3, SM2 mixing scheme 1 +// WB2, WB3, SM2 mixing scheme void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index); + +// Sherman-Morrison-Woodbury kernel 2 +// WB2, WB3, SM3 mixing scheme +void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); + +// Sherman-Morrison-Woodbury kernel 3 +// WB2, WB3, SM4 mixing scheme +void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); + +// Sherman-Morrison-Woodbury kernel 4 +// WB2, SM2 mixing scheme +void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); diff --git a/src/Helpers.cpp b/src/Helpers.cpp index 12ce105..8f8dac6 100644 --- a/src/Helpers.cpp +++ b/src/Helpers.cpp @@ -3,7 +3,7 @@ // Set common break-down threshold double threshold() { const double threshold = THRESHOLD; -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Break-down threshold set to: " << threshold << std::endl; #endif return threshold; @@ -26,7 +26,7 @@ void selectLargestDenominator(unsigned int l, unsigned int N_updates, index = p[j]; component = Updates_index[index - 1]; breakdown = std::fabs(1 + ylk[l][index][component]); -#ifdef DEBUG +#ifdef DEBUG2 std::cout << "Inside selectLargestDenominator()" << std::endl; std::cout << "breakdown = fabs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << std::endl; @@ -58,4 +58,4 @@ lapack_int inverse(double *A, unsigned n) { ret = LAPACKE_dgetri(LAPACK_COL_MAJOR, n, A, n, ipiv); return ret; } -#endif \ No newline at end of file +#endif diff --git a/src/SMWB.cpp b/src/SMWB.cpp index a112dd7..c353259 100644 --- a/src/SMWB.cpp +++ b/src/SMWB.cpp @@ -4,9 +4,13 @@ #include "Helpers.hpp" // Sherman-Morrison-Woodbury kernel 1 -// WB2, WB3, SM2 mixing scheme 1 +// WB2, WB3, SM2 mixing scheme void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { - std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; +#ifdef DEBUG2 + std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; + showMatrix2(Updates_index, 1, N_updates, "Updates_index"); + showMatrix2(Updates, N_updates, Dim, "Updates"); +#endif unsigned int n_of_3blocks = N_updates / 3; unsigned int remainder = N_updates % 3; @@ -14,8 +18,6 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double unsigned int length_2block = 2 * Dim; unsigned int length_1block = 1 * Dim; - // std::cerr << "Number of blocks: " << n_of_3blocks << ". Remainder: " << remainder << "." << std::endl; - // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel if (n_of_3blocks > 0) { for (unsigned int i = 0; i < n_of_3blocks; i++) { @@ -30,11 +32,16 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double bool ok; ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block); if (!ok) { // Send the entire block to SM2 +#ifdef DEBUG2 std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM2" << std::endl; + showMatrix2(Updates_3block, 3, Dim, "Updates_3block"); + showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block"); +#endif SM2(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block); } } } +/// Convert asrray copies to pointers if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel double Updates_2block[length_2block]; @@ -47,7 +54,9 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double bool ok; ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); if (!ok) { // Send the entire block to SM2 +#ifdef DEBUG2 std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl; +#endif SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); } } else if (remainder == 1) { // Apply last remaining update with SM2 @@ -63,9 +72,207 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double } } +// Sherman-Morrison-Woodbury kernel 2 +// WB2, WB3, SM3 mixing scheme +void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { +#ifdef DEBUG2 + std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; +#endif + + unsigned int n_of_3blocks = N_updates / 3; + unsigned int remainder = N_updates % 3; + unsigned int length_3block = 3 * Dim; + unsigned int length_2block = 2 * Dim; + unsigned int length_1block = 1 * Dim; + +#ifdef DEBUG2 + showMatrix2(Updates_index, 1, N_updates, "Updates_index"); + showMatrix2(Updates, N_updates, Dim, "Updates"); +#endif + + // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel + if (n_of_3blocks > 0) { + for (unsigned int i = 0; i < n_of_3blocks; i++) { + double Updates_3block[length_3block]; + unsigned int Updates_index_3block[3]; + Updates_index_3block[0] = Updates_index[3 * i + 0]; + Updates_index_3block[1] = Updates_index[3 * i + 1]; + Updates_index_3block[2] = Updates_index[3 * i + 2]; + for (unsigned int j = 0; j < length_3block; j++) { + Updates_3block[j] = Updates[i * length_3block + j]; + } + bool ok; + ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block); + if (!ok) { // Send the entire block to SM3 +#ifdef DEBUG2 + std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM3" << std::endl; + showMatrix2(Updates_3block, 3, Dim, "Updates_3block"); + showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block"); +#endif + SM3(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block); + } + } + } + + if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel + double Updates_2block[length_2block]; + unsigned int Updates_index_2block[2]; + Updates_index_2block[0] = Updates_index[3 * n_of_3blocks + 0]; + Updates_index_2block[1] = Updates_index[3 * n_of_3blocks + 1]; + for (unsigned int i = 0; i < length_2block; i++) { + Updates_2block[i] = Updates[n_of_3blocks * length_3block + i]; + } + bool ok; + ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); + if (!ok) { // Send the entire block to SM3 + std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM3" << std::endl; + SM3(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); + } + } else if (remainder == 1) { // Apply last remaining update with SM3 + double Updates_1block[length_1block]; + unsigned int Updates_index_1block[1]; + Updates_index_1block[0] = Updates_index[3 * n_of_3blocks + 0]; + for (unsigned int i = 0; i < length_1block; i++) { + Updates_1block[i] = Updates[n_of_3blocks * length_3block + i]; + } + SM3(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); + } else { // remainder == 0 + // Nothing left to do. + } +} + +// Sherman-Morrison-Woodbury kernel 3 +// WB2, WB3, SM4 mixing scheme +void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { + std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; + + unsigned int n_of_3blocks = N_updates / 3; + unsigned int remainder = N_updates % 3; + unsigned int length_3block = 3 * Dim; + unsigned int length_2block = 2 * Dim; + unsigned int length_1block = 1 * Dim; + +#ifdef DEBUG2 + showMatrix2(Updates_index, 1, N_updates, "Updates_index"); + showMatrix2(Updates, N_updates, Dim, "Updates"); +#endif + + // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel + if (n_of_3blocks > 0) { + for (unsigned int i = 0; i < n_of_3blocks; i++) { + double Updates_3block[length_3block]; + unsigned int Updates_index_3block[3]; + Updates_index_3block[0] = Updates_index[3 * i + 0]; + Updates_index_3block[1] = Updates_index[3 * i + 1]; + Updates_index_3block[2] = Updates_index[3 * i + 2]; + for (unsigned int j = 0; j < length_3block; j++) { + Updates_3block[j] = Updates[i * length_3block + j]; + } + bool ok; + ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block); + if (!ok) { // Send the entire block to SM4 + std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM4" << std::endl; +#ifdef DEBUG2 + showMatrix2(Updates_3block, 3, Dim, "Updates_3block"); + showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block"); +#endif + SM4(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block); + } + } + } + + if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel + double Updates_2block[length_2block]; + unsigned int Updates_index_2block[2]; + Updates_index_2block[0] = Updates_index[3 * n_of_3blocks + 0]; + Updates_index_2block[1] = Updates_index[3 * n_of_3blocks + 1]; + for (unsigned int i = 0; i < length_2block; i++) { + Updates_2block[i] = Updates[n_of_3blocks * length_3block + i]; + } + bool ok; + ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); + if (!ok) { // Send the entire block to SM4 + std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM4" << std::endl; + SM4(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); + } + } else if (remainder == 1) { // Apply last remaining update with SM4 + double Updates_1block[length_1block]; + unsigned int Updates_index_1block[1]; + Updates_index_1block[0] = Updates_index[3 * n_of_3blocks + 0]; + for (unsigned int i = 0; i < length_1block; i++) { + Updates_1block[i] = Updates[n_of_3blocks * length_3block + i]; + } + SM4(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); + } else { // remainder == 0 + // Nothing left to do. + } +} + +// Sherman-Morrison-Woodbury kernel 4 +// WB2, SM2 mixing scheme +void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { +#ifdef DEBUG2 + std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; + showMatrix2(Updates_index, 1, N_updates, "Updates_index"); + showMatrix2(Updates, N_updates, Dim, "Updates"); +#endif + + unsigned int n_of_2blocks = N_updates / 2; + unsigned int remainder = N_updates % 2; + unsigned int length_2block = 2 * Dim; + unsigned int length_1block = 1 * Dim; + + // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel + if (n_of_2blocks > 0) { + for (unsigned int i = 0; i < n_of_2blocks; i++) { + double Updates_2block[length_2block]; + unsigned int Updates_index_2block[2]; + Updates_index_2block[0] = Updates_index[2 * i + 0]; + Updates_index_2block[1] = Updates_index[2 * i + 1]; + for (unsigned int j = 0; j < length_2block; j++) { + Updates_2block[j] = Updates[i * length_2block + j]; + } + bool ok; + ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); + if (!ok) { // Send the entire block to SM2 + std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl; +#ifdef DEBUG2 + showMatrix2(Updates_2block, 2, Dim, "Updates_2block"); + showMatrix2(Updates_index_2block, 1, 2, "Updates_index_2block"); +#endif + SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); + } + } + } + + if (remainder == 1) { // Apply last remaining update with SM4 + double Updates_1block[length_1block]; + unsigned int Updates_index_1block[1]; + Updates_index_1block[0] = Updates_index[2 * n_of_2blocks + 0]; + for (unsigned int i = 0; i < length_1block; i++) { + Updates_1block[i] = Updates[n_of_2blocks * length_2block + i]; + } + SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); + } else { // remainder == 0 + // Nothing left to do. + } +} + extern "C" { void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, double **linUpdates, unsigned int **Updates_index) { SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); } +void SMWB2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, + double **linUpdates, unsigned int **Updates_index) { + SMWB2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} +void SMWB3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, + double **linUpdates, unsigned int **Updates_index) { + SMWB3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} +void SMWB4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, + double **linUpdates, unsigned int **Updates_index) { + SMWB4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} } diff --git a/src/SM_Maponi.cpp b/src/SM_Maponi.cpp index 65afefa..19351c1 100644 --- a/src/SM_Maponi.cpp +++ b/src/SM_Maponi.cpp @@ -40,7 +40,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // Calculate the {y_{0,k}} for (k = 1; k < N_updates + 1; k++) { -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Compute y0k: " << std::endl; std::cerr << "ylk[0][" << k << "][:]"; std::cerr << std::endl; @@ -51,14 +51,14 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, Updates[(k - 1) * Dim + (j - 1)]; } } -#ifdef DEBUG +#ifdef DEBUG2 showVector(ylk[0][k], Dim, ""); #endif } // Calculate the {y_{l,k}} from the {y_{0,k}} for (l = 0; l < N_updates; l++) { -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "In outer compute-ylk-loop: l = " << l << std::endl; std::cerr << std::endl; #endif @@ -69,7 +69,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // Select component and comp. bd-condition. component = Updates_index[p[l + 1] - 1]; beta = 1 + ylk[l][p[l + 1]][component]; -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "p[l+1] = " << p[l + 1] << std::endl; std::cerr << "component = " << component << std::endl; std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component @@ -83,7 +83,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double ibeta = 1.0 / beta; // Compute intermediate update to Slater_inv -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Compute intermediate update to Slater_inv" << std::endl; std::cerr << "component = " << component << std::endl; std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component @@ -102,14 +102,14 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, tmp = next; next = last; last = tmp; -#ifdef DEBUG +#ifdef DEBUG2 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] * ibeta; -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Inside k-loop: k = " << k << std::endl; std::cerr << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl; std::cerr << std::endl; @@ -172,7 +172,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // Calculate the y_{0,k} = S_0^{-1} u_k for (k = 1; k < N_updates + 1; k++) { -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Compute y0k: " << std::endl; std::cerr << "ylk[0][" << k << "][:]"; std::cerr << std::endl; @@ -183,14 +183,14 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, Updates[(k - 1) * Dim + (j - 1)]; } } -#ifdef DEBUG +#ifdef DEBUG2 showVector(ylk[0][k], Dim + 1, ""); #endif } // Calculate the {y_{l,k}} from the {y_{0,k}} for (l = 0; l < N_updates; l++) { -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "In outer compute-ylk-loop: l = " << l << std::endl; std::cerr << std::endl; #endif @@ -201,7 +201,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // Select component and comp. bd-condition. component = Updates_index[p[l + 1] - 1]; beta = 1 + ylk[l][p[l + 1]][component]; -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "p[l+1] = " << p[l + 1] << std::endl; std::cerr << "component = " << component << std::endl; std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component @@ -222,7 +222,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double ibeta = 1.0 / beta; // Compute intermediate update to Slater_inv -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Compute intermediate update to Slater_inv" << std::endl; std::cerr << "component = " << component << std::endl; std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component @@ -241,14 +241,14 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, tmp = next; next = last; last = tmp; -#ifdef DEBUG +#ifdef DEBUG2 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] * ibeta; -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Inside k-loop: k = " << k << std::endl; std::cerr << "ylk[" << l + 1 << "][" << p[k] << "][:]"; std::cerr << std::endl; diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp index c68dd81..42fb75c 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -6,7 +6,10 @@ // Naïve Sherman Morrison void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { +#ifdef DEBUG1 std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl; +#endif + double C[Dim]; double D[Dim]; @@ -50,7 +53,9 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // http://hdl.handle.net/10919/52966 void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { +#ifdef DEBUG1 std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl; +#endif double C[Dim]; double D[Dim]; @@ -72,8 +77,11 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // Denominator double den = 1 + C[Updates_index[l] - 1]; if (std::fabs(den) < threshold()) { +#ifdef DEBUG1 std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; + std::cerr << "Denominator = " << den << std::endl; +#endif // U_l = U_l / 2 (do the split) for (unsigned int i = 0; i < Dim; i++) { diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp index dacf676..474cf3d 100644 --- a/src/Woodbury.cpp +++ b/src/Woodbury.cpp @@ -18,7 +18,9 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, B := 1 + V * C, 2 x 2 D := V * S^{-1}, 2 x dim */ + #ifdef DEBUG1 std::cerr << "Called Woodbury 2x2 kernel" << std::endl; +#endif // Construct V from Updates_index unsigned int V[2 * Dim]; // 2 x Dim matrix stored in row-major order @@ -59,7 +61,10 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, // Check if determinant of inverted matrix is not zero double det = Binv[0] * Binv[3] - Binv[1] * Binv[2]; if (std::fabs(det) < threshold()) { - std::cerr << "Determinant approaching 0!" << std::endl; +#ifdef DEBUG1 + std::cerr << "Determinant too close to zero!" << std::endl; + std::cerr << "Determinant = " << det << std::endl; +#endif return false; } @@ -87,9 +92,10 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, B := 1 + V * C, 3 x 3 D := V * S^{-1}, 3 x dim */ +#ifdef DEBUG1 std::cerr << "Called Woodbury 3x3 kernel" << std::endl; - -#ifdef DEBUG +#endif +#ifdef DEBUG2 showMatrix2(Slater_inv, Dim, Dim, "Slater_inv BEFORE update"); showMatrix2(Updates, 3, Dim, "Updates"); showMatrix2(Updates_index, 1, 3, "Updates_index"); @@ -102,7 +108,7 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, V[Dim + Updates_index[1] - 1] = 1; V[2 * Dim + Updates_index[2] - 1] = 1; -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(V, 3, Dim, "V"); #endif @@ -118,14 +124,14 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, } } -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(C, Dim, 3, "C = S_inv * U"); #endif // Compute D = V * S^{-1} double D[3 * Dim]; matMul2(V, Slater_inv, D, 3, Dim, Dim); -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(D, 3, Dim, "D = V * S_inv"); #endif @@ -136,7 +142,7 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, B[4] += 1; B[8] += 1; -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(B, 3, 3, "B = 1 + V * C"); #endif @@ -147,7 +153,7 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, B[2] * (B[3] * B[7] - B[4] * B[6]); idet = 1.0 / det; -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Determinant of B = " << det << std::endl; #endif @@ -161,7 +167,8 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, Binv[7] = - ( B[0] * B[7] - B[6] * B[1] ) * idet; Binv[8] = ( B[0] * B[4] - B[3] * B[1] ) * idet; -#ifdef DEBUG +#ifdef DEBUG2 + std::cerr << "Conditioning number of B = " << condition1(B, Binv, 3) << std::endl; showMatrix2(Binv, 3, 3, "Binv"); #endif @@ -171,12 +178,15 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, Binv[1] * (Binv[3] * Binv[8] - Binv[5] * Binv[6]) + Binv[2] * (Binv[3] * Binv[7] - Binv[4] * Binv[6]); -#ifdef DEBUG +#ifdef DEBUG2 std::cerr << "Determinant of Binv = " << det << std::endl; #endif if (std::fabs(det) < threshold()) { - std::cerr << "Determinant approached 0!" << std::endl; +#ifdef DEBUG1 + std::cerr << "Determinant too close to zero!" << std::endl; + std::cerr << "Determinant = " << det << std::endl; +#endif return false; } @@ -184,7 +194,7 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, double tmp[3 * Dim]; matMul2(Binv, D, tmp, 3, 3, Dim); -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(tmp, 3, Dim, "tmp = Binv * D"); #endif @@ -192,7 +202,7 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, double tmp2[Dim * Dim]; matMul2(C, tmp, tmp2, Dim, 3, Dim); -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(tmp2, Dim, Dim, "tmp2 = C * tmp"); #endif @@ -201,7 +211,7 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, Slater_inv[i] -= tmp2[i]; } -#ifdef DEBUG +#ifdef DEBUG2 showMatrix2(Slater_inv, Dim, Dim, "Slater_inv AFTER update"); #endif return true; diff --git a/tests/run_test.sh b/tests/run_test.sh index a4b81a9..9b5c185 100755 --- a/tests/run_test.sh +++ b/tests/run_test.sh @@ -7,7 +7,7 @@ ## Call: $ ./run_test.sh <{sm1|sm2|sm3|maponia3}> ## Example: ./run_test.sh sm2 1 500 -TEST=test_h5 +TEST=test_h5.O2 ALGO=$1 START=$2 STOP=$3 diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp index f750f10..f3509c6 100644 --- a/tests/test_h5.cpp +++ b/tests/test_h5.cpp @@ -7,6 +7,7 @@ #include "SMWB.hpp" using namespace H5; + // #define DEBUG const H5std_string FILE_NAME("dataset.hdf5"); @@ -48,7 +49,7 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { double *u = new double[nupdates * dim]; /* Test */ -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_inverse, dim, "OLD Inverse"); #endif @@ -62,11 +63,11 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { } } -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_matrix, dim, "OLD Slater"); #endif -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(u, dim, "Updates"); #endif @@ -84,6 +85,12 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { SM4(slater_inverse, dim, nupdates, u, col_update_index); } else if (version == "smwb1") { SMWB1(slater_inverse, dim, nupdates, u, col_update_index); + } else if (version == "smwb2") { + SMWB2(slater_inverse, dim, nupdates, u, col_update_index); + } else if (version == "smwb3") { + SMWB3(slater_inverse, dim, nupdates, u, col_update_index); + } else if (version == "smwb4") { + SMWB4(slater_inverse, dim, nupdates, u, col_update_index); #ifdef MKL } else if (version == "lapack") { memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double)); @@ -94,11 +101,11 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { exit(1); } -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_matrix, dim, "NEW Slater"); #endif -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_inverse, dim, "NEW Inverse"); #endif @@ -124,7 +131,7 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { std::cout << "Residual = " << version << " " << cycle << " " << res_max << " " << res2 << std::endl; -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(res, dim, "Result"); #endif diff --git a/tests/vfc_test_h5.cpp b/tests/vfc_test_h5.cpp index e567b9f..f30f32d 100644 --- a/tests/vfc_test_h5.cpp +++ b/tests/vfc_test_h5.cpp @@ -89,11 +89,11 @@ int test_cycle(H5File file, int cycle, std::string version, double *u = new double[nupdates * dim]; /* Test */ -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_matrix, dim, "OLD Slater"); #endif -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_inverse, dim, "OLD Inverse"); #endif @@ -121,11 +121,11 @@ int test_cycle(H5File file, int cycle, std::string version, exit(1); } -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_matrix, dim, "NEW Slater"); #endif -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(slater_inverse, dim, "NEW Inverse"); #endif @@ -136,7 +136,7 @@ int test_cycle(H5File file, int cycle, std::string version, double res2 = residual2(res, dim); double frob2 = norm_frobenius2(res, dim); -#ifdef DEBUG +#ifdef DEBUG2 showMatrix(res, dim, "Result"); #endif