From 67f5379beab5db4698d4c45b834a9c4e7fdf7786 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 7 Jul 2021 12:06:31 +0200 Subject: [PATCH] - Moved check on determinants in Woodbury kernels before inversion of B that was there by mistake. - Split SM2 into SM2 and SM2star so that keeps all updates for later when used in combination with the Woodbury kernels to improve the numerical accuracy to the same level as that of SM2. --- include/SMWB.hpp | 16 +-- include/SM_Standard.hpp | 3 + src/SMWB.cpp | 265 ++++++++++++++++++++++------------------ src/SM_Standard.cpp | 68 ++++++++++- src/Woodbury.cpp | 60 ++++----- tests/test_h5.cpp | 18 +-- 6 files changed, 255 insertions(+), 175 deletions(-) diff --git a/include/SMWB.hpp b/include/SMWB.hpp index 28adda4..7ba4fdd 100644 --- a/include/SMWB.hpp +++ b/include/SMWB.hpp @@ -3,15 +3,15 @@ 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 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 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 diff --git a/include/SM_Standard.hpp b/include/SM_Standard.hpp index 4bfabd0..f68b122 100644 --- a/include/SM_Standard.hpp +++ b/include/SM_Standard.hpp @@ -7,6 +7,9 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index); +void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index, double *later_updates, unsigned int* later_index, unsigned int *later); + // Sherman Morrison, leaving zero denominators for later void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index); diff --git a/src/SMWB.cpp b/src/SMWB.cpp index b354a2d..b45b093 100644 --- a/src/SMWB.cpp +++ b/src/SMWB.cpp @@ -3,6 +3,9 @@ #include "Woodbury.hpp" #include "Helpers.hpp" +// #define DEBUG1 +// #define DEBUG2 + // Sherman-Morrison-Woodbury kernel 1 // WB2, WB3, SM2 mixing scheme void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { @@ -19,6 +22,9 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double 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 + double later_updates[Dim * N_updates]; + unsigned int later_index[N_updates]; + unsigned int later = 0; if (n_of_3blocks > 0) { for (unsigned int i = 0; i < n_of_3blocks; i++) { double *Updates_3block = &Updates[i * length_3block]; @@ -26,12 +32,14 @@ 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 + #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); + #endif + unsigned int l = 0; + SM2star(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block, later_updates + (Dim * later), later_index + later, &l); + later = later + l; } } } @@ -42,125 +50,132 @@ 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 + #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); + #endif + unsigned int l = 0; + SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block, later_updates+(Dim * later), later_index + later, &l); + later = later + l; } } else if (remainder == 1) { // Apply last remaining update with SM2 double *Updates_1block = &Updates[n_of_3blocks * length_3block]; unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; - SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); - } else { // remainder == 0 - // Nothing left to do. + unsigned int l = 0; + SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates+(Dim * later), later_index + later, &l); + later = later + l; } + + if (later > 0) { + SM2(Slater_inv, Dim, later, later_updates, later_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) { -#ifdef DEBUG2 - std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; -#endif +// // 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; +// 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 +// #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 = &Updates[i * length_3block]; - unsigned int *Updates_index_3block = &Updates_index[i * 3]; - 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); - } - } - } +// // 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 = &Updates[i * length_3block]; +// unsigned int *Updates_index_3block = &Updates_index[i * 3]; +// 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 = &Updates[n_of_3blocks * length_3block]; - unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; - 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 = &Updates[n_of_3blocks * length_3block]; - unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; - SM3(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); - } else { // remainder == 0 - // Nothing left to do. - } -} +// if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel +// double *Updates_2block = &Updates[n_of_3blocks * length_3block]; +// unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; +// 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 = &Updates[n_of_3blocks * length_3block]; +// unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; +// 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; +// // 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; +// 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 +// #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 = &Updates[i * length_3block]; - unsigned int *Updates_index_3block = &Updates_index[i * 3]; - 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); - } - } - } +// // 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 = &Updates[i * length_3block]; +// unsigned int *Updates_index_3block = &Updates_index[i * 3]; +// 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 = &Updates[n_of_3blocks * length_3block]; - unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; - 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 = &Updates[n_of_3blocks * length_3block]; - unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; - SM4(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); - } else { // remainder == 0 - // Nothing left to do. - } -} +// if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel +// double *Updates_2block = &Updates[n_of_3blocks * length_3block]; +// unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; +// 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 = &Updates[n_of_3blocks * length_3block]; +// unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; +// 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 @@ -177,6 +192,9 @@ void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double unsigned int length_1block = 1 * Dim; // Apply first 2*n_of_2blocks updates in n_of_2blocks blocks of 2 updates with Woodbury 2x2 kernel + double later_updates[Dim * N_updates]; + unsigned int later_index[N_updates]; + unsigned int later = 0; if (n_of_2blocks > 0) { for (unsigned int i = 0; i < n_of_2blocks; i++) { double *Updates_2block = &Updates[i * length_2block]; @@ -185,22 +203,29 @@ void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double 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 + #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); + #endif + unsigned int l = 0; + SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block, later_updates + (Dim * later), later_index + later, &l); + later = later + l; } } } - if (remainder == 1) { // Apply last remaining update with SM4 + if (remainder == 1) { // Apply last remaining update with SM2 double *Updates_1block = &Updates[n_of_2blocks * length_2block]; unsigned int *Updates_index_1block = &Updates_index[2 * n_of_2blocks]; - SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); - } else { // remainder == 0 - // Nothing left to do. + unsigned int l = 0; + SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates + (Dim * later), later_index + later, &l); + later = later + l; } + + if (later > 0) { + SM2(Slater_inv, Dim, later, later_updates, later_index); + } + } extern "C" { @@ -208,14 +233,14 @@ 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 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_Standard.cpp b/src/SM_Standard.cpp index 7423ce8..b2a40dd 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -3,6 +3,8 @@ #include "SM_Standard.hpp" #include "Helpers.hpp" +#define DEBUG1 + // Naïve Sherman Morrison void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { @@ -57,13 +59,13 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl; #endif - double C[Dim]; - double D[Dim]; - double later_updates[Dim * N_updates]; unsigned int later_index[N_updates]; unsigned int later = 0; + double C[Dim]; + double D[Dim]; + unsigned int l = 0; // For each update while (l < N_updates) { @@ -116,6 +118,66 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, } } +// Sherman Morrison, with J. Slagel splitting +// http://hdl.handle.net/10919/52966 +void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index, double *later_updates, unsigned int* later_index, unsigned int *later) { +#ifdef DEBUG1 + std::cerr << "Called SM2* with " << N_updates << " updates" << std::endl; +#endif + + double C[Dim]; + double D[Dim]; + + unsigned int l = 0; + // For each update + while (l < N_updates) { + // C = A^{-1} x U_l + for (unsigned int i = 0; i < Dim; i++) { + C[i] = 0; + for (unsigned int j = 0; j < Dim; j++) { + C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; + } + } + + // 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++) { + later_updates[*later * Dim + i] = Updates[l * Dim + i] / 2.0; + C[i] /= 2.0; + } + later_index[*later] = Updates_index[l]; + (*later)++; + + den = 1 + C[Updates_index[l] - 1]; + } + double iden = 1 / den; + + // D = v^T x A^{-1} + for (unsigned int j = 0; j < Dim; j++) { + D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; + } + + // A^{-1} = A^{-1} - C x D / den + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + double update = C[i] * D[j] * iden; + Slater_inv[i * Dim + j] -= update; + } + } + l += 1; + } +} + + // Sherman Morrison, leaving zero denominators for later void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp index f6573f2..b5b299a 100644 --- a/src/Woodbury.cpp +++ b/src/Woodbury.cpp @@ -10,6 +10,9 @@ #include "Woodbury.hpp" #include "Helpers.hpp" +// #define DEBUG1 +// #define DEBUG2 + // Woodbury 2x2 kernel bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, unsigned int *Updates_index) { @@ -18,7 +21,7 @@ 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 +#ifdef DEBUG1 std::cerr << "Called Woodbury 2x2 kernel" << std::endl; #endif @@ -43,7 +46,6 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, } } } - // matMul2(Updates, Slater_inv, C, 2, Dim, Dim); // Compute B = 1 + V * C double B[4]; @@ -53,24 +55,23 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, B[3] = C[row2 * 2 + 1]; B[0] += 1, B[3] += 1; - // Compute B^{-1} with explicit formula for 2x2 inversion - double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]); - double Binv[4]; - Binv[0] = idet * B[3]; - Binv[1] = -1.0 * idet * B[1]; - Binv[2] = -1.0 * idet * B[2]; - Binv[3] = idet * B[0]; - // Check if determinant of inverted matrix is not zero - double det = Binv[0] * Binv[3] - Binv[1] * Binv[2]; + double det = B[0] * B[3] - B[1] * B[2]; if (std::fabs(det) < threshold()) { + std::cerr << "Determinant too close to zero! No inverse found." << std::endl; #ifdef DEBUG1 - std::cerr << "Determinant too close to zero!" << std::endl; std::cerr << "Determinant = " << det << std::endl; #endif return false; } + // Compute B^{-1} with explicit formula for 2x2 inversion + double Binv[4], idet = 1.0 / det; + Binv[0] = idet * B[3]; + Binv[1] = -1.0 * idet * B[1]; + Binv[2] = -1.0 * idet * B[2]; + Binv[3] = idet * B[0]; + // Compute B^{-1} x D double tmp[2 * Dim]; matMul2(Binv, D, tmp, 2, 2, Dim); @@ -127,7 +128,6 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, } } } - // matMul2(Updates, Slater_inv, C, 2, Dim, Dim); #ifdef DEBUG2 showMatrix2(C, Dim, 3, "C = S_inv * U"); @@ -154,17 +154,25 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, showMatrix2(B, 3, 3, "B = 1 + V * C"); #endif - // Compute B^{-1} with explicit formula for 3x3 inversion - double Binv[9], det, idet; + // Check if determinant of B is not too close to zero + double det; det = B[0] * (B[4] * B[8] - B[5] * B[7]) - B[1] * (B[3] * B[8] - B[5] * B[6]) + B[2] * (B[3] * B[7] - B[4] * B[6]); - idet = 1.0 / det; - #ifdef DEBUG2 std::cerr << "Determinant of B = " << det << std::endl; #endif + if (std::fabs(det) < threshold()) { + // if (std::fabs(det) < 1000000) { + std::cerr << "Determinant too close to zero! No inverse found." << std::endl; +#ifdef DEBUG1 + std::cerr << "Determinant = " << det << std::endl; +#endif + return false; + } + // Compute B^{-1} with explicit formula for 3x3 inversion + double Binv[9], idet = 1.0 / det; Binv[0] = ( B[4] * B[8] - B[7] * B[5] ) * idet; Binv[1] = - ( B[1] * B[8] - B[7] * B[2] ) * idet; Binv[2] = ( B[1] * B[5] - B[4] * B[2] ) * idet; @@ -180,24 +188,6 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, showMatrix2(Binv, 3, 3, "Binv"); #endif - // Check if determinant of inverted matrix is not zero - // double det; - det = Binv[0] * (Binv[4] * Binv[8] - Binv[5] * Binv[7]) - - Binv[1] * (Binv[3] * Binv[8] - Binv[5] * Binv[6]) + - Binv[2] * (Binv[3] * Binv[7] - Binv[4] * Binv[6]); - -#ifdef DEBUG2 - std::cerr << "Determinant of Binv = " << det << std::endl; -#endif - - if (std::fabs(det) < threshold()) { -#ifdef DEBUG1 - std::cerr << "Determinant too close to zero!" << std::endl; - std::cerr << "Determinant = " << det << std::endl; -#endif - return false; - } - // Compute B^{-1} x D double tmp[3 * Dim]; matMul2(Binv, D, tmp, 3, 3, Dim); diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp index 9cb296b..c174f36 100644 --- a/tests/test_h5.cpp +++ b/tests/test_h5.cpp @@ -7,7 +7,7 @@ #include "Woodbury.hpp" #include "SMWB.hpp" -#define PERF +// #define PERF #ifdef PERF unsigned int repetition_number; @@ -101,10 +101,10 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { WB3(slater_inverse_nonpersistent, dim, u, col_update_index); } else if (version == "smwb1") { SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); - } else if (version == "smwb2") { - SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); - } else if (version == "smwb3") { - SMWB3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); + // } else if (version == "smwb2") { + // SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); + // } else if (version == "smwb3") { + // SMWB3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); } else if (version == "smwb4") { SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); #ifdef MKL @@ -138,10 +138,10 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { WB3(slater_inverse, dim, 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 == "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