mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-11-03 20:54:08 +01:00
Merge pull request #46 from fmgjcoppens/bug/woodbury
Solved bug in Woodbury and SM/WB mix kernels
This commit is contained in:
François Coppens
GitHub
commit
61df6e2827
@ -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
|
||||
|
||||
@ -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);
|
||||
|
||||
265
src/SMWB.cpp
265
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);
|
||||
|
||||
@ -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) {
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user