mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-11-03 20:54:08 +01:00
- Added various Sherman-Morrison-Woodbury kernels
- Separated debug information into 2 different debug levels.
This commit is contained in:
Francois Coppens
parent
4b34fd7e17
commit
837e160f17
2
Makefile
2
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++
|
||||
|
||||
@ -7,6 +7,7 @@
|
||||
#ifdef MKL
|
||||
#include <mkl_lapacke.h>
|
||||
#endif
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
#ifndef THRESHOLD
|
||||
@ -275,3 +276,38 @@ template <typename T> T residual2(T * A, unsigned int Dim) {
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
// Computes the condition number of A using a previously computed B=A^{-1}
|
||||
template <typename T> 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 <typename T> 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
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -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
|
||||
#endif
|
||||
|
||||
215
src/SMWB.cpp
215
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);
|
||||
}
|
||||
}
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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++) {
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
## Call: $ ./run_test.sh <{sm1|sm2|sm3|maponia3}> <start cycle> <stop cycle>
|
||||
## Example: ./run_test.sh sm2 1 500
|
||||
|
||||
TEST=test_h5
|
||||
TEST=test_h5.O2
|
||||
ALGO=$1
|
||||
START=$2
|
||||
STOP=$3
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user