TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2024-12-26 14:23:47 +01:00
Code Issues Projects Releases Wiki Activity

- Moved check on determinants in Woodbury kernels before inversion of B that was there by mistake.

Browse Source 
- 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.
This commit is contained in:
François Coppens 2021-07-07 12:06:31 +02:00
parent 6f282f329c
commit 67f5379bea
6 changed files with 255 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
include/SMWB.hpp
View File

@ -3,15 +3,15 @@
void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index); double *Updates, unsigned int *Updates_index);
// Sherman-Morrison-Woodbury kernel 2 // // Sherman-Morrison-Woodbury kernel 2
// WB2, WB3, SM3 mixing scheme // // WB2, WB3, SM3 mixing scheme
void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index); // double *Updates, unsigned int *Updates_index);
// Sherman-Morrison-Woodbury kernel 3 // // Sherman-Morrison-Woodbury kernel 3
// WB2, WB3, SM4 mixing scheme // // WB2, WB3, SM4 mixing scheme
void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index); // double *Updates, unsigned int *Updates_index);
// Sherman-Morrison-Woodbury kernel 4 // Sherman-Morrison-Woodbury kernel 4
// WB2, SM2 mixing scheme // WB2, SM2 mixing scheme

3
include/SM_Standard.hpp
View File

@ -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, void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index); 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 // Sherman Morrison, leaving zero denominators for later
void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index); double *Updates, unsigned int *Updates_index);

265
src/SMWB.cpp
View File

@ -3,6 +3,9 @@
#include "Woodbury.hpp" #include "Woodbury.hpp"
#include "Helpers.hpp" #include "Helpers.hpp"
// #define DEBUG1
// #define DEBUG2
// Sherman-Morrison-Woodbury kernel 1 // Sherman-Morrison-Woodbury kernel 1
// WB2, WB3, SM2 mixing scheme // WB2, WB3, SM2 mixing scheme
void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { 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; 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 // 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) { if (n_of_3blocks > 0) {
for (unsigned int i = 0; i < n_of_3blocks; i++) { for (unsigned int i = 0; i < n_of_3blocks; i++) {
double *Updates_3block = &Updates[i * length_3block]; 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; bool ok;
ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block); ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
if (!ok) { // Send the entire block to SM2 if (!ok) { // Send the entire block to SM2
#ifdef DEBUG2 #ifdef DEBUG2
std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM2" << std::endl; std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM2" << std::endl;
showMatrix2(Updates_3block, 3, Dim, "Updates_3block"); showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block"); showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
#endif #endif
SM2(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block); 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; bool ok;
ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
if (!ok) { // Send the entire block to SM2 if (!ok) { // Send the entire block to SM2
#ifdef DEBUG2 #ifdef DEBUG2
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl; std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl;
#endif #endif
SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); 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 } else if (remainder == 1) { // Apply last remaining update with SM2
double *Updates_1block = &Updates[n_of_3blocks * length_3block]; double *Updates_1block = &Updates[n_of_3blocks * length_3block];
unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); unsigned int l = 0;
} else { // remainder == 0 SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates+(Dim * later), later_index + later, &l);
// Nothing left to do. later = later + l;
} }
if (later > 0) {
SM2(Slater_inv, Dim, later, later_updates, later_index);
}
} }
// Sherman-Morrison-Woodbury kernel 2 // // Sherman-Morrison-Woodbury kernel 2
// WB2, WB3, SM3 mixing scheme // // WB2, WB3, SM3 mixing scheme
void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { // void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
#ifdef DEBUG2 // #ifdef DEBUG2
std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; // std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
#endif // #endif
unsigned int n_of_3blocks = N_updates / 3; // unsigned int n_of_3blocks = N_updates / 3;
unsigned int remainder = N_updates % 3; // unsigned int remainder = N_updates % 3;
unsigned int length_3block = 3 * Dim; // unsigned int length_3block = 3 * Dim;
unsigned int length_2block = 2 * Dim; // unsigned int length_2block = 2 * Dim;
unsigned int length_1block = 1 * Dim; // unsigned int length_1block = 1 * Dim;
#ifdef DEBUG2 // #ifdef DEBUG2
showMatrix2(Updates_index, 1, N_updates, "Updates_index"); // showMatrix2(Updates_index, 1, N_updates, "Updates_index");
showMatrix2(Updates, N_updates, Dim, "Updates"); // showMatrix2(Updates, N_updates, Dim, "Updates");
#endif // #endif
// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel // // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
if (n_of_3blocks > 0) { // if (n_of_3blocks > 0) {
for (unsigned int i = 0; i < n_of_3blocks; i++) { // for (unsigned int i = 0; i < n_of_3blocks; i++) {
double *Updates_3block = &Updates[i * length_3block]; // double *Updates_3block = &Updates[i * length_3block];
unsigned int *Updates_index_3block = &Updates_index[i * 3]; // unsigned int *Updates_index_3block = &Updates_index[i * 3];
bool ok; // bool ok;
ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block); // ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
if (!ok) { // Send the entire block to SM3 // if (!ok) { // Send the entire block to SM3
#ifdef DEBUG2 // #ifdef DEBUG2
std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM3" << std::endl; // std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM3" << std::endl;
showMatrix2(Updates_3block, 3, Dim, "Updates_3block"); // showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block"); // showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
#endif // #endif
SM3(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block); // SM3(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
} // }
} // }
} // }
if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel // if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
double *Updates_2block = &Updates[n_of_3blocks * length_3block]; // double *Updates_2block = &Updates[n_of_3blocks * length_3block];
unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; // unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
bool ok; // bool ok;
ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); // ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
if (!ok) { // Send the entire block to SM3 // if (!ok) { // Send the entire block to SM3
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM3" << std::endl; // std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM3" << std::endl;
SM3(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); // SM3(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
} // }
} else if (remainder == 1) { // Apply last remaining update with SM3 // } else if (remainder == 1) { // Apply last remaining update with SM3
double *Updates_1block = &Updates[n_of_3blocks * length_3block]; // double *Updates_1block = &Updates[n_of_3blocks * length_3block];
unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; // unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
SM3(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); // SM3(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
} else { // remainder == 0 // } else { // remainder == 0
// Nothing left to do. // // Nothing left to do.
} // }
} // }
// Sherman-Morrison-Woodbury kernel 3 // // Sherman-Morrison-Woodbury kernel 3
// WB2, WB3, SM4 mixing scheme // // WB2, WB3, SM4 mixing scheme
void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { // 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; // std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
unsigned int n_of_3blocks = N_updates / 3; // unsigned int n_of_3blocks = N_updates / 3;
unsigned int remainder = N_updates % 3; // unsigned int remainder = N_updates % 3;
unsigned int length_3block = 3 * Dim; // unsigned int length_3block = 3 * Dim;
unsigned int length_2block = 2 * Dim; // unsigned int length_2block = 2 * Dim;
unsigned int length_1block = 1 * Dim; // unsigned int length_1block = 1 * Dim;
#ifdef DEBUG2 // #ifdef DEBUG2
showMatrix2(Updates_index, 1, N_updates, "Updates_index"); // showMatrix2(Updates_index, 1, N_updates, "Updates_index");
showMatrix2(Updates, N_updates, Dim, "Updates"); // showMatrix2(Updates, N_updates, Dim, "Updates");
#endif // #endif
// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel // // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
if (n_of_3blocks > 0) { // if (n_of_3blocks > 0) {
for (unsigned int i = 0; i < n_of_3blocks; i++) { // for (unsigned int i = 0; i < n_of_3blocks; i++) {
double *Updates_3block = &Updates[i * length_3block]; // double *Updates_3block = &Updates[i * length_3block];
unsigned int *Updates_index_3block = &Updates_index[i * 3]; // unsigned int *Updates_index_3block = &Updates_index[i * 3];
bool ok; // bool ok;
ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block); // ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
if (!ok) { // Send the entire block to SM4 // if (!ok) { // Send the entire block to SM4
std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM4" << std::endl; // std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM4" << std::endl;
#ifdef DEBUG2 // #ifdef DEBUG2
showMatrix2(Updates_3block, 3, Dim, "Updates_3block"); // showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block"); // showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
#endif // #endif
SM4(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block); // SM4(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
} // }
} // }
} // }
if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel // if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
double *Updates_2block = &Updates[n_of_3blocks * length_3block]; // double *Updates_2block = &Updates[n_of_3blocks * length_3block];
unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; // unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
bool ok; // bool ok;
ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block); // ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
if (!ok) { // Send the entire block to SM4 // if (!ok) { // Send the entire block to SM4
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM4" << std::endl; // std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM4" << std::endl;
SM4(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); // SM4(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
} // }
} else if (remainder == 1) { // Apply last remaining update with SM4 // } else if (remainder == 1) { // Apply last remaining update with SM4
double *Updates_1block = &Updates[n_of_3blocks * length_3block]; // double *Updates_1block = &Updates[n_of_3blocks * length_3block];
unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; // unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
SM4(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); // SM4(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
} else { // remainder == 0 // } else { // remainder == 0
// Nothing left to do. // // Nothing left to do.
} // }
} // }
// Sherman-Morrison-Woodbury kernel 4 // Sherman-Morrison-Woodbury kernel 4
// WB2, SM2 mixing scheme // 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; 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 // 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) { if (n_of_2blocks > 0) {
for (unsigned int i = 0; i < n_of_2blocks; i++) { for (unsigned int i = 0; i < n_of_2blocks; i++) {
double *Updates_2block = &Updates[i * length_2block]; 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); ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
if (!ok) { // Send the entire block to SM2 if (!ok) { // Send the entire block to SM2
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl; 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_2block, 2, Dim, "Updates_2block");
showMatrix2(Updates_index_2block, 1, 2, "Updates_index_2block"); showMatrix2(Updates_index_2block, 1, 2, "Updates_index_2block");
#endif #endif
SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block); 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]; double *Updates_1block = &Updates[n_of_2blocks * length_2block];
unsigned int *Updates_index_1block = &Updates_index[2 * n_of_2blocks]; unsigned int *Updates_index_1block = &Updates_index[2 * n_of_2blocks];
SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block); unsigned int l = 0;
} else { // remainder == 0 SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates + (Dim * later), later_index + later, &l);
// Nothing left to do. later = later + l;
} }
if (later > 0) {
SM2(Slater_inv, Dim, later, later_updates, later_index);
}
} }
extern "C" { 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) { double **linUpdates, unsigned int **Updates_index) {
SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
} }
void SMWB2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, // void SMWB2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
double **linUpdates, unsigned int **Updates_index) { // double **linUpdates, unsigned int **Updates_index) {
SMWB2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); // SMWB2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
} // }
void SMWB3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, // void SMWB3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
double **linUpdates, unsigned int **Updates_index) { // double **linUpdates, unsigned int **Updates_index) {
SMWB3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); // SMWB3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
} // }
void SMWB4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, void SMWB4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
double **linUpdates, unsigned int **Updates_index) { double **linUpdates, unsigned int **Updates_index) {
SMWB4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); SMWB4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);

68
src/SM_Standard.cpp
View File

@ -3,6 +3,8 @@
#include "SM_Standard.hpp" #include "SM_Standard.hpp"
#include "Helpers.hpp" #include "Helpers.hpp"
#define DEBUG1
// Naïve Sherman Morrison // Naïve Sherman Morrison
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index) { 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; std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl;
#endif #endif
double C[Dim];
double D[Dim];
double later_updates[Dim * N_updates]; double later_updates[Dim * N_updates];
unsigned int later_index[N_updates]; unsigned int later_index[N_updates];
unsigned int later = 0; unsigned int later = 0;
double C[Dim];
double D[Dim];
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) { 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 // Sherman Morrison, leaving zero denominators for later
void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index) { double *Updates, unsigned int *Updates_index) {

60
src/Woodbury.cpp
View File

@ -10,6 +10,9 @@
#include "Woodbury.hpp" #include "Woodbury.hpp"
#include "Helpers.hpp" #include "Helpers.hpp"
// #define DEBUG1
// #define DEBUG2
// Woodbury 2x2 kernel // Woodbury 2x2 kernel
bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
unsigned int *Updates_index) { 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 B := 1 + V * C, 2 x 2
D := V * S^{-1}, 2 x dim D := V * S^{-1}, 2 x dim
*/ */
#ifdef DEBUG1 #ifdef DEBUG1
std::cerr << "Called Woodbury 2x2 kernel" << std::endl; std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
#endif #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 // Compute B = 1 + V * C
double B[4]; double B[4];
@ -53,24 +55,23 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
B[3] = C[row2 * 2 + 1]; B[3] = C[row2 * 2 + 1];
B[0] += 1, B[3] += 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 // 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()) { if (std::fabs(det) < threshold()) {
std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
#ifdef DEBUG1 #ifdef DEBUG1
std::cerr << "Determinant too close to zero!" << std::endl;
std::cerr << "Determinant = " << det << std::endl; std::cerr << "Determinant = " << det << std::endl;
#endif #endif
return false; 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 // Compute B^{-1} x D
double tmp[2 * Dim]; double tmp[2 * Dim];
matMul2(Binv, D, tmp, 2, 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 #ifdef DEBUG2
showMatrix2(C, Dim, 3, "C = S_inv * U"); 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"); showMatrix2(B, 3, 3, "B = 1 + V * C");
#endif #endif
// Compute B^{-1} with explicit formula for 3x3 inversion // Check if determinant of B is not too close to zero
double Binv[9], det, idet; double det;
det = B[0] * (B[4] * B[8] - B[5] * B[7]) - det = B[0] * (B[4] * B[8] - B[5] * B[7]) -
B[1] * (B[3] * B[8] - B[5] * B[6]) + B[1] * (B[3] * B[8] - B[5] * B[6]) +
B[2] * (B[3] * B[7] - B[4] * B[6]); B[2] * (B[3] * B[7] - B[4] * B[6]);
idet = 1.0 / det;
#ifdef DEBUG2 #ifdef DEBUG2
std::cerr << "Determinant of B = " << det << std::endl; std::cerr << "Determinant of B = " << det << std::endl;
#endif #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[0] = ( B[4] * B[8] - B[7] * B[5] ) * idet;
Binv[1] = - ( B[1] * B[8] - B[7] * B[2] ) * idet; Binv[1] = - ( B[1] * B[8] - B[7] * B[2] ) * idet;
Binv[2] = ( B[1] * B[5] - B[4] * 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"); showMatrix2(Binv, 3, 3, "Binv");
#endif #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 // Compute B^{-1} x D
double tmp[3 * Dim]; double tmp[3 * Dim];
matMul2(Binv, D, tmp, 3, 3, Dim); matMul2(Binv, D, tmp, 3, 3, Dim);

18
tests/test_h5.cpp
View File

@ -7,7 +7,7 @@
#include "Woodbury.hpp" #include "Woodbury.hpp"
#include "SMWB.hpp" #include "SMWB.hpp"
#define PERF // #define PERF
#ifdef PERF #ifdef PERF
unsigned int repetition_number; 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); WB3(slater_inverse_nonpersistent, dim, u, col_update_index);
} else if (version == "smwb1") { } else if (version == "smwb1") {
SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
} else if (version == "smwb2") { // } else if (version == "smwb2") {
SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); // SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
} else if (version == "smwb3") { // } else if (version == "smwb3") {
SMWB3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); // SMWB3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
} else if (version == "smwb4") { } else if (version == "smwb4") {
SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
#ifdef MKL #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); WB3(slater_inverse, dim, u, col_update_index);
} else if (version == "smwb1") { } else if (version == "smwb1") {
SMWB1(slater_inverse, dim, nupdates, u, col_update_index); SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
} else if (version == "smwb2") { // } else if (version == "smwb2") {
SMWB2(slater_inverse, dim, nupdates, u, col_update_index); // SMWB2(slater_inverse, dim, nupdates, u, col_update_index);
} else if (version == "smwb3") { // } else if (version == "smwb3") {
SMWB3(slater_inverse, dim, nupdates, u, col_update_index); // SMWB3(slater_inverse, dim, nupdates, u, col_update_index);
} else if (version == "smwb4") { } else if (version == "smwb4") {
SMWB4(slater_inverse, dim, nupdates, u, col_update_index); SMWB4(slater_inverse, dim, nupdates, u, col_update_index);
#ifdef MKL #ifdef MKL