mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-12-25 13:53:56 +01:00
- Passing break-down threshold as a function argument
- Renaming kernels to correspond with the ones in QMCkl - In the qmckl-version of the test program, chaning the way integer data is read from the HDF5 file.
This commit is contained in:
Francois Coppens
parent
6ce2055e59
commit
74bb333de1
7
Makefile
7
Makefile
@ -2,7 +2,7 @@
|
||||
ifeq ($(ENV),INTEL)
|
||||
CXX = icpx
|
||||
FC = ifort
|
||||
ARCH = -march=native
|
||||
ARCH = -xCORE-AVX2
|
||||
OPT = -O3
|
||||
DEBUG = -g -debug full
|
||||
else ifeq ($(ENV),LLVM)
|
||||
@ -14,8 +14,9 @@ else ifeq ($(ENV),LLVM)
|
||||
else ifeq ($(ENV),GNU)
|
||||
CXX = g++
|
||||
FC = gfortran
|
||||
ARCH = -mavx
|
||||
OPT = -O3
|
||||
# ARCH = -mavx
|
||||
ARCH =
|
||||
OPT = -O1
|
||||
DEBUG = -g
|
||||
else
|
||||
$(error No valid compiler environment set in $$ENV. \
|
||||
|
||||
@ -8,6 +8,9 @@
|
||||
#include <mkl_lapacke.h>
|
||||
#endif
|
||||
|
||||
#include <cstdio>
|
||||
|
||||
|
||||
// #define DEBUG
|
||||
#ifndef THRESHOLD
|
||||
#define THRESHOLD 1e-3
|
||||
@ -86,7 +89,8 @@ template <typename T> T *transpose(T *A, unsigned int M) {
|
||||
return B;
|
||||
}
|
||||
|
||||
template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
|
||||
template <typename T>
|
||||
void matMul(T *A, T *B, T *C, unsigned int M) {
|
||||
memset(C, 0, M * M * sizeof(T));
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
@ -263,7 +267,7 @@ template <typename T> T residual_frobenius2(T *A, unsigned int Dim) {
|
||||
res += delta * delta;
|
||||
}
|
||||
}
|
||||
return res;
|
||||
return sqrt(res);
|
||||
}
|
||||
|
||||
template <typename T> T residual2(T *A, unsigned int Dim) {
|
||||
|
||||
@ -1,11 +1,17 @@
|
||||
// Sherman-Morrison-Woodbury kernel 1
|
||||
// WB2, WB3, SM2 mixing scheme
|
||||
void SMWB1(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index);
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel 2
|
||||
// Sherman-Morrison-Woodbury kernel WB2s
|
||||
// WB2, SM2 mixing scheme
|
||||
void SMWB2(double *Slater_inv, const unsigned int Dim,
|
||||
void WB2s(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index);
|
||||
unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel WB3s
|
||||
// WB3, SM2 mixing scheme
|
||||
void WB3s(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel WB32s
|
||||
// WB3, WB2, SM2 mixing scheme
|
||||
void WB32s(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
@ -1,22 +1,22 @@
|
||||
// Naïve Sherman Morrison
|
||||
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
// Sherman Morrison, with J. Slagel splitting
|
||||
// http://hdl.handle.net/10919/52966
|
||||
void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
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);
|
||||
unsigned int *later, const double breakdown);
|
||||
|
||||
// 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);
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
// Sherman Morrison (SM3+SM2), leaving zero denominators for later (SM3), and
|
||||
// when none are left falling back on Splitting (SM2)
|
||||
void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index);
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// Woodbury 2x2 kernel
|
||||
bool WB2(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
const unsigned int *Updates_index);
|
||||
const unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
// Woodbury 3x3 kernel
|
||||
bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
const unsigned int *Updates_index);
|
||||
const unsigned int *Updates_index, const double breakdown);
|
||||
|
||||
2
qmckl
2
qmckl
@ -1 +1 @@
|
||||
Subproject commit 41be86fe594f561247986cd75dc4557e0bc5e6de
|
||||
Subproject commit a5e58c80d70978c12d24ff47b2362c9145af325c
|
||||
245
src/SMWB.cpp
245
src/SMWB.cpp
@ -3,88 +3,14 @@
|
||||
#include "SM_Standard.hpp"
|
||||
#include "Woodbury.hpp"
|
||||
|
||||
// #define DEBUG1
|
||||
#define DEBUG1
|
||||
// #define DEBUG2
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel 1
|
||||
// WB2, WB3, SM2 mixing scheme
|
||||
void SMWB1(double *Slater_inv, const unsigned int Dim,
|
||||
const 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_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;
|
||||
|
||||
// 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];
|
||||
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 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
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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 SM2
|
||||
#ifdef DEBUG2
|
||||
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2"
|
||||
<< std::endl;
|
||||
#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];
|
||||
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
|
||||
// Sherman-Morrison-Woodbury kernel WB2s
|
||||
// WB2, SM2 mixing scheme
|
||||
void SMWB2(double *Slater_inv, const unsigned int Dim,
|
||||
void WB2s(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index) {
|
||||
unsigned int *Updates_index, const double breakdown) {
|
||||
#ifdef DEBUG2
|
||||
std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates
|
||||
<< " updates" << std::endl;
|
||||
@ -95,7 +21,6 @@ void SMWB2(double *Slater_inv, const unsigned int Dim,
|
||||
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 2*n_of_2blocks updates in n_of_2blocks blocks of 2 updates with
|
||||
// Woodbury 2x2 kernel
|
||||
@ -107,10 +32,10 @@ void SMWB2(double *Slater_inv, const unsigned int Dim,
|
||||
double *Updates_2block = &Updates[i * length_2block];
|
||||
unsigned int *Updates_index_2block = &Updates_index[i * 2];
|
||||
bool ok;
|
||||
ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
|
||||
ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block, breakdown);
|
||||
if (!ok) { // Send the entire block to SM2
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2star"
|
||||
std::cerr << "Woodbury 2x2 block failed! Sending to SM w/ US"
|
||||
<< std::endl;
|
||||
#endif
|
||||
#ifdef DEBUG2
|
||||
@ -119,33 +44,169 @@ void SMWB2(double *Slater_inv, const unsigned int Dim,
|
||||
#endif
|
||||
unsigned int l = 0;
|
||||
SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block,
|
||||
later_updates + (Dim * later), later_index + later, &l);
|
||||
later_updates + (Dim * later), later_index + later, &l, breakdown);
|
||||
later = later + l;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (remainder == 1) { // Apply last remaining update with SM2
|
||||
if (remainder != 0) { // 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];
|
||||
unsigned int l = 0;
|
||||
SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block,
|
||||
later_updates + (Dim * later), later_index + later, &l);
|
||||
SM2star(Slater_inv, Dim, remainder, Updates_1block, Updates_index_1block,
|
||||
later_updates + (Dim * later), later_index + later, &l, breakdown);
|
||||
later = later + l;
|
||||
}
|
||||
|
||||
if (later > 0) {
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index);
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel WB3s
|
||||
// WB3, SM2 mixing scheme
|
||||
void WB3s(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index, const double breakdown) {
|
||||
#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;
|
||||
unsigned int length_3block = 3 * 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];
|
||||
unsigned int *Updates_index_3block = &Updates_index[i * 3];
|
||||
bool ok;
|
||||
ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block, breakdown);
|
||||
if (!ok) { // Send the entire block to SM2
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Woodbury 3x3 block failed! Sending to SM w/ US"
|
||||
<< std::endl;
|
||||
#endif
|
||||
#ifdef DEBUG2
|
||||
showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
|
||||
showMatrix2(Updates_index_3block, 1, 3, "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, breakdown);
|
||||
later = later + l;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (remainder != 0) { // 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];
|
||||
unsigned int l = 0;
|
||||
SM2star(Slater_inv, Dim, remainder, Updates_2block, Updates_index_2block,
|
||||
later_updates + (Dim * later), later_index + later, &l, breakdown);
|
||||
later = later + l;
|
||||
}
|
||||
|
||||
if (later > 0) {
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel WB32s
|
||||
// WB3, WB2, SM2 mixing scheme
|
||||
void WB32s(double *Slater_inv, const unsigned int Dim,
|
||||
const unsigned int N_updates, double *Updates,
|
||||
unsigned int *Updates_index, const double breakdown) {
|
||||
#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;
|
||||
unsigned int length_3block = 3 * 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];
|
||||
unsigned int *Updates_index_3block = &Updates_index[i * 3];
|
||||
bool ok;
|
||||
ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block, breakdown);
|
||||
if (!ok) { // Send the entire block to SM2
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Woodbury 3x3 block failed! Sending to SM w/ US"
|
||||
<< std::endl;
|
||||
#endif
|
||||
#ifdef DEBUG2
|
||||
showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
|
||||
showMatrix2(Updates_index_3block, 1, 3, "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, breakdown);
|
||||
later = later + l;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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, breakdown);
|
||||
if (!ok) { // Send the entire block to SM2
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Woodbury 2x2 block failed! Sending to SM w/ US"
|
||||
<< std::endl;
|
||||
#endif
|
||||
unsigned int l = 0;
|
||||
SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block,
|
||||
later_updates + (Dim * later), later_index + later, &l, breakdown);
|
||||
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];
|
||||
unsigned int l = 0;
|
||||
SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block,
|
||||
later_updates + (Dim * later), later_index + later, &l, breakdown);
|
||||
later = later + l;
|
||||
}
|
||||
|
||||
if (later > 0) {
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
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 WB2s_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
WB2s(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
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 WB3s_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
WB3s(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
void WB32s_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
WB32s(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
@ -6,8 +6,12 @@
|
||||
// #define DEBUG1
|
||||
|
||||
// Naïve Sherman Morrison
|
||||
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
void SM1(double *Slater_inv,
|
||||
unsigned int Dim,
|
||||
unsigned int N_updates,
|
||||
double *Updates,
|
||||
unsigned int *Updates_index,
|
||||
const double breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
@ -28,7 +32,7 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (std::fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
@ -51,12 +55,13 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
l += 1;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
// Sherman Morrison, with J. Slagel splitting
|
||||
// http://hdl.handle.net/10919/52966
|
||||
void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
@ -66,10 +71,10 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
unsigned int later = 0;
|
||||
|
||||
SM2star(Slater_inv, Dim, N_updates, Updates, Updates_index, later_updates,
|
||||
later_index, &later);
|
||||
later_index, &later, breakdown);
|
||||
|
||||
if (later > 0) {
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index);
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
@ -78,7 +83,7 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
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) {
|
||||
unsigned int *later, const double breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM2* with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
@ -99,7 +104,7 @@ void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (std::fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
@ -136,7 +141,7 @@ void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// 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) {
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM3 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
@ -161,7 +166,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (std::fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
@ -201,7 +206,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
}
|
||||
// If some have failed, make a recursive call
|
||||
else if (later > 0) {
|
||||
SM3(Slater_inv, Dim, later, later_updates, later_index);
|
||||
SM3(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
@ -209,7 +214,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
// Leave zero denominators for later (SM3), and when none are left then split
|
||||
// (SM2)
|
||||
void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
double *Updates, unsigned int *Updates_index, const double breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM4 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
@ -234,7 +239,7 @@ void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (std::fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
@ -266,32 +271,32 @@ void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// If all the updates have failed, fall back on splitting (SM2)
|
||||
if (later == N_updates) {
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index);
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
// If some have failed, make a recursive call
|
||||
else if (later > 0) {
|
||||
SM4(Slater_inv, Dim, later, later_updates, later_index);
|
||||
SM4(Slater_inv, Dim, later, later_updates, later_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
extern "C" {
|
||||
void SM1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index) {
|
||||
SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
|
||||
void SM2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index) {
|
||||
SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
|
||||
void SM3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index) {
|
||||
SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
|
||||
void SM4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index) {
|
||||
SM4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
double **linUpdates, unsigned int **Updates_index, const double breakdown) {
|
||||
SM4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
|
||||
}
|
||||
}
|
||||
|
||||
@ -15,7 +15,7 @@
|
||||
|
||||
// Woodbury 2x2 kernel
|
||||
bool WB2(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
const unsigned int *Updates_index) {
|
||||
const unsigned int *Updates_index, const double breakdown) {
|
||||
/*
|
||||
C := S^{-1} * U, dim x 2
|
||||
B := 1 + V * C, 2 x 2
|
||||
@ -48,7 +48,7 @@ bool WB2(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
|
||||
// Check if determinant of inverted matrix is not zero
|
||||
double det = B0 * B3 - B1 * B2;
|
||||
if (std::fabs(det) < threshold()) {
|
||||
if (std::fabs(det) < breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Determinant too close to zero! No inverse found."
|
||||
<< std::endl;
|
||||
@ -86,7 +86,7 @@ bool WB2(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
|
||||
// Woodbury 3x3 kernel
|
||||
bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
const unsigned int *Updates_index) {
|
||||
const unsigned int *Updates_index, const double breakdown) {
|
||||
/*
|
||||
C := S^{-1} * U, dim x 3
|
||||
B := 1 + V * C, 3 x 3
|
||||
@ -139,7 +139,7 @@ bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
#ifdef DEBUG2
|
||||
std::cerr << "Determinant of B = " << det << std::endl;
|
||||
#endif
|
||||
if (std::fabs(det) < threshold()) {
|
||||
if (std::fabs(det) < breakdown) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Determinant too close to zero! No inverse found."
|
||||
<< std::endl;
|
||||
@ -197,15 +197,15 @@ bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
|
||||
|
||||
extern "C" {
|
||||
bool WB2_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
unsigned int **Updates_index, const double breakdown) {
|
||||
bool ok;
|
||||
ok = WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
|
||||
ok = WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index, breakdown);
|
||||
return ok;
|
||||
}
|
||||
bool WB3_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
unsigned int **Updates_index, const double breakdown) {
|
||||
bool ok;
|
||||
ok = WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
|
||||
ok = WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index, breakdown);
|
||||
return ok;
|
||||
}
|
||||
}
|
||||
|
||||
@ -3,7 +3,6 @@
|
||||
|
||||
#include "Helpers.hpp"
|
||||
#include "qmckl.h"
|
||||
|
||||
#include <math.h>
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
@ -16,9 +15,9 @@ unsigned int repetition_number;
|
||||
|
||||
const H5std_string FILE_NAME("dataset.hdf5");
|
||||
|
||||
void read_int(H5::H5File file, std::string key, uint64_t *data) {
|
||||
void read_int(H5::H5File file, std::string key, unsigned int *data) {
|
||||
H5::DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, H5::PredType::STD_U64LE);
|
||||
ds.read(data, H5::PredType::STD_U32LE);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
@ -35,10 +34,12 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
std::string group = "cycle_" + std::to_string(cycle);
|
||||
|
||||
unsigned int col, i, j;
|
||||
unsigned int dim_32, nupdates_32;
|
||||
uint64_t dim, nupdates;
|
||||
|
||||
read_int(file, group + "/slater_matrix_dim", &dim);
|
||||
read_int(file, group + "/nupdates", &nupdates);
|
||||
read_int(file, group + "/slater_matrix_dim", &dim_32);
|
||||
read_int(file, group + "/nupdates", &nupdates_32);
|
||||
dim = dim_32; nupdates = nupdates_32;
|
||||
|
||||
double *slater_matrix = new double[dim * dim];
|
||||
read_double(file, group + "/slater_matrix", slater_matrix);
|
||||
@ -46,8 +47,12 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
double *slater_inverse = new double[dim * dim];
|
||||
read_double(file, group + "/slater_inverse", slater_inverse);
|
||||
|
||||
unsigned int *temp = new unsigned int[nupdates];
|
||||
uint64_t *col_update_index = new uint64_t[nupdates];
|
||||
read_int(file, group + "/col_update_index", col_update_index);
|
||||
read_int(file, group + "/col_update_index", temp);
|
||||
for (i = 0; i < nupdates; i++) {
|
||||
col_update_index[i] = temp[i];
|
||||
}
|
||||
|
||||
double *updates = new double[nupdates * dim];
|
||||
read_double(file, group + "/updates", updates);
|
||||
@ -74,7 +79,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
std::cout << "# of reps. = " << repetition_number << std::endl;
|
||||
double *slater_inverse_nonpersistent = new double[dim * dim];
|
||||
|
||||
if (version == "qmckl_sm1") {
|
||||
if (version == "sm1") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -82,7 +87,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
u, col_update_index, breakdown, slater_inverse_nonpersistent);
|
||||
}
|
||||
}
|
||||
else if (version == "qmckl_wb2") {
|
||||
else if (version == "wb2") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -90,7 +95,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
u, col_update_index, breakdown, slater_inverse_nonpersistent);
|
||||
}
|
||||
}
|
||||
else if (version == "qmckl_wb3") {
|
||||
else if (version == "wb3") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -98,7 +103,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
u, col_update_index, breakdown, slater_inverse_nonpersistent);
|
||||
}
|
||||
}
|
||||
else if (version == "qmckl_sm2") {
|
||||
else if (version == "sm2") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -106,7 +111,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
u, col_update_index, breakdown, slater_inverse_nonpersistent);
|
||||
}
|
||||
}
|
||||
else if (version == "qmckl_wb2s") {
|
||||
else if (version == "wb2s") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -114,7 +119,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
u, col_update_index, breakdown, slater_inverse_nonpersistent);
|
||||
}
|
||||
}
|
||||
else if (version == "qmckl_wb3s") {
|
||||
else if (version == "wb3s") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -122,7 +127,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
u, col_update_index, breakdown, slater_inverse_nonpersistent);
|
||||
}
|
||||
}
|
||||
else if (version == "qmckl_wb32s") {
|
||||
else if (version == "wb32s") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
@ -138,21 +143,21 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
dim * dim * sizeof(double));
|
||||
delete[] slater_inverse_nonpersistent;
|
||||
#else // No performance measurements repetition
|
||||
if (version == "qmckl_sm1") {
|
||||
if (version == "sm1") {
|
||||
qmckl_context context;
|
||||
context = qmckl_context_create();
|
||||
qmckl_exit_code rc;
|
||||
rc = qmckl_sherman_morrison_c(context, dim, nupdates,
|
||||
u, col_update_index, breakdown, slater_inverse);
|
||||
}
|
||||
else if (version == "qmckl_wb2") {
|
||||
else if (version == "wb2") {
|
||||
qmckl_context context;
|
||||
context = qmckl_context_create();
|
||||
qmckl_exit_code rc;
|
||||
rc = qmckl_woodbury_2_c(context, dim,
|
||||
u, col_update_index, breakdown, slater_inverse);
|
||||
}
|
||||
else if (version == "qmckl_wb3") {
|
||||
else if (version == "wb3") {
|
||||
qmckl_context context;
|
||||
context = qmckl_context_create();
|
||||
qmckl_exit_code rc;
|
||||
@ -166,8 +171,68 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
|
||||
#endif // PERF
|
||||
delete[] u, col_update_index;
|
||||
|
||||
showMatrix(slater_matrix, dim, "Slater Matrix");
|
||||
showMatrix(slater_inverse, dim, "Slater Inverse");
|
||||
double *res = new double[dim * dim]{0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
{
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
for (unsigned int k = 0; k < dim; k++) {
|
||||
res[i * dim + j] += slater_matrix[i * dim + k] * slater_inverse[k * dim + j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//matMul2(slater_matrix, slater_inverse, res, dim_32, dim_32, dim_32);
|
||||
//
|
||||
//
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
printf("[");
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
if (slater_matrix[i * dim + j] >= 0) {
|
||||
printf(" %17.10e,", slater_matrix[i * dim + j]);
|
||||
} else {
|
||||
printf(" %17.10e,", slater_matrix[i * dim + j]);
|
||||
}
|
||||
}
|
||||
printf(" ],\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
printf("[");
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
if (slater_inverse[i * dim + j] >= 0) {
|
||||
printf(" %17.10e,", slater_inverse[i * dim + j]);
|
||||
} else {
|
||||
printf(" %17.10e,", slater_inverse[i * dim + j]);
|
||||
}
|
||||
}
|
||||
printf(" ],\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
printf("[");
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
if (res[i * dim + j] >= 0) {
|
||||
printf(" %17.10e,", res[i * dim + j]);
|
||||
} else {
|
||||
printf(" %17.10e,", res[i * dim + j]);
|
||||
}
|
||||
}
|
||||
printf(" ],\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
//
|
||||
//
|
||||
bool ok = is_identity(res, dim, tolerance);
|
||||
double res_max = residual_max(res, dim);
|
||||
double res2 = residual_frobenius2(res, dim);
|
||||
|
||||
@ -13,31 +13,29 @@
|
||||
unsigned int repetition_number;
|
||||
#endif
|
||||
|
||||
using namespace H5;
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
const H5std_string FILE_NAME("dataset.hdf5");
|
||||
|
||||
void read_int(H5File file, std::string key, unsigned int *data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::STD_U32LE);
|
||||
void read_int(H5::H5File file, std::string key, unsigned int *data) {
|
||||
H5::DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, H5::PredType::STD_U32LE);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
void read_double(H5File file, std::string key, double *data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::IEEE_F64LE);
|
||||
void read_double(H5::H5File file, std::string key, double *data) {
|
||||
H5::DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, H5::PredType::IEEE_F64LE);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown, double tolerance) {
|
||||
|
||||
/* Read the data */
|
||||
|
||||
std::string group = "cycle_" + std::to_string(cycle);
|
||||
|
||||
unsigned int dim, nupdates, col, i, j;
|
||||
|
||||
|
||||
read_int(file, group + "/slater_matrix_dim", &dim);
|
||||
read_int(file, group + "/nupdates", &nupdates);
|
||||
|
||||
@ -56,9 +54,6 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
double *u = new double[nupdates * dim];
|
||||
|
||||
/* Test */
|
||||
#ifdef DEBUG2
|
||||
showMatrix(slater_inverse, dim, "OLD Inverse");
|
||||
#endif
|
||||
|
||||
// Transform replacement updates in 'updates[]' into additive updates in 'u[]'
|
||||
for (j = 0; j < nupdates; j++) {
|
||||
@ -69,54 +64,79 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
slater_matrix[i * dim + (col - 1)] = updates[i + j * dim];
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(slater_matrix, dim, "OLD Slater");
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(u, dim, "Updates");
|
||||
#endif
|
||||
delete[] updates;
|
||||
|
||||
#ifdef PERF
|
||||
std::cout << "# of reps. = " << repetition_number << std::endl;
|
||||
double *slater_inverse_nonpersistent = new double[dim * dim];
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
if (version == "maponia3") {
|
||||
MaponiA3(slater_inverse_nonpersistent, dim, nupdates, u,
|
||||
col_update_index);
|
||||
} else if (version == "maponia3s") {
|
||||
MaponiA3S(slater_inverse_nonpersistent, dim, nupdates, u,
|
||||
col_update_index);
|
||||
} else if (version == "sm1") {
|
||||
SM1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm2") {
|
||||
SM2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm3") {
|
||||
SM3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm4") {
|
||||
SM4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "wb2") {
|
||||
WB2(slater_inverse_nonpersistent, dim, u, col_update_index);
|
||||
} else if (version == "wb3") {
|
||||
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);
|
||||
#ifdef MKL
|
||||
} else if (version == "lapack") {
|
||||
memcpy(slater_inverse_nonpersistent, slater_matrix,
|
||||
dim * dim * sizeof(double));
|
||||
inverse(slater_inverse_nonpersistent, dim);
|
||||
#endif // MKL
|
||||
} else {
|
||||
std::cerr << "Unknown version " << version << std::endl;
|
||||
exit(1);
|
||||
|
||||
if (version == "sm1") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
SM1(slater_inverse_nonpersistent, dim, nupdates,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
else if (version == "wb2") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
WB2(slater_inverse_nonpersistent, dim,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
else if (version == "wb3") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
WB3(slater_inverse_nonpersistent, dim,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
else if (version == "sm2") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
SM2(slater_inverse_nonpersistent, dim, nupdates,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
else if (version == "wb2s") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
WB2s(slater_inverse_nonpersistent, dim, nupdates,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
else if (version == "wb3s") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
WB3s(slater_inverse_nonpersistent, dim, nupdates,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
else if (version == "wb32s") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
memcpy(slater_inverse_nonpersistent, slater_inverse,
|
||||
dim * dim * sizeof(double));
|
||||
WB32s(slater_inverse_nonpersistent, dim, nupdates,
|
||||
u, col_update_index, breakdown);
|
||||
}
|
||||
}
|
||||
#ifdef MKL
|
||||
else if (version == "lapack") {
|
||||
memcpy(slater_inverse_nonpersistent, slater_matrix,
|
||||
dim * dim * sizeof(double));
|
||||
inverse(slater_inverse_nonpersistent, dim);
|
||||
}
|
||||
#endif // MKL
|
||||
else {
|
||||
std::cerr << "Unknown version " << version << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
std::memcpy(slater_inverse, slater_inverse_nonpersistent,
|
||||
dim * dim * sizeof(double));
|
||||
delete[] slater_inverse_nonpersistent;
|
||||
@ -137,10 +157,12 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
WB2(slater_inverse, dim, u, col_update_index);
|
||||
} else if (version == "wb3") {
|
||||
WB3(slater_inverse, dim, u, col_update_index);
|
||||
} else if (version == "smwb1") {
|
||||
SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "smwb4") {
|
||||
SMWB4(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "wb2s") {
|
||||
WB2s(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "wb3s") {
|
||||
WB3s(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "wb32s") {
|
||||
WB32s(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
#ifdef MKL
|
||||
} else if (version == "lapack") {
|
||||
memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double));
|
||||
@ -151,17 +173,70 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
exit(1);
|
||||
}
|
||||
#endif // PERF
|
||||
delete[] u, col_update_index;
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(slater_matrix, dim, "NEW Slater");
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(slater_inverse, dim, "NEW Inverse");
|
||||
#endif
|
||||
|
||||
showMatrix(slater_matrix, dim, "Slater Matrix");
|
||||
showMatrix(slater_inverse, dim, "Slater Inverse");
|
||||
double *res = new double[dim * dim]{0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
{
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
for (unsigned int k = 0; k < dim; k++) {
|
||||
res[i * dim + j] += slater_matrix[i * dim + k] * slater_inverse[k * dim + j];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//matMul2(slater_matrix, slater_inverse, res, dim, dim, dim);
|
||||
//
|
||||
//
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
printf("[");
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
if (slater_matrix[i * dim + j] >= 0) {
|
||||
printf(" %17.10e,", slater_matrix[i * dim + j]);
|
||||
} else {
|
||||
printf(" %17.10e,", slater_matrix[i * dim + j]);
|
||||
}
|
||||
}
|
||||
printf(" ],\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
printf("[");
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
if (slater_inverse[i * dim + j] >= 0) {
|
||||
printf(" %17.10e,", slater_inverse[i * dim + j]);
|
||||
} else {
|
||||
printf(" %17.10e,", slater_inverse[i * dim + j]);
|
||||
}
|
||||
}
|
||||
printf(" ],\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
for (unsigned int i = 0; i < dim; i++) {
|
||||
printf("[");
|
||||
for (unsigned int j = 0; j < dim; j++) {
|
||||
if (res[i * dim + j] >= 0) {
|
||||
printf(" %17.10e,", res[i * dim + j]);
|
||||
} else {
|
||||
printf(" %17.10e,", res[i * dim + j]);
|
||||
}
|
||||
}
|
||||
printf(" ],\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
//
|
||||
//
|
||||
bool ok = is_identity(res, dim, tolerance);
|
||||
double res_max = residual_max(res, dim);
|
||||
double res2 = residual_frobenius2(res, dim);
|
||||
@ -173,25 +248,25 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
showMatrix(res, dim, "Result");
|
||||
#endif
|
||||
|
||||
delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse;
|
||||
delete[] res, slater_matrix, slater_inverse;
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
#ifdef PERF
|
||||
if (argc != 6) {
|
||||
if (argc != 7) {
|
||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||
std::cerr
|
||||
<< "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance> <number of reps.>"
|
||||
<< "usage: test_h5 <version> <start cycle> <stop cycle> <break-down threshold> <tolerance> <number of reps.>"
|
||||
<< std::endl;
|
||||
return 1;
|
||||
}
|
||||
#else
|
||||
if (argc != 5) {
|
||||
if (argc != 6) {
|
||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||
std::cerr
|
||||
<< "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>"
|
||||
<< "usage: test_h5 <version> <start cycle> <stop cycle> <break-down threshold> <tolerance>"
|
||||
<< std::endl;
|
||||
return 1;
|
||||
}
|
||||
@ -200,16 +275,17 @@ int main(int argc, char **argv) {
|
||||
std::string version(argv[1]);
|
||||
int start_cycle = std::stoi(argv[2]);
|
||||
int stop_cycle = std::stoi(argv[3]);
|
||||
double tolerance = std::stod(argv[4]);
|
||||
H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
||||
double breakdown = std::stod(argv[4]);
|
||||
double tolerance = std::stod(argv[5]);
|
||||
H5::H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
||||
|
||||
#ifdef PERF
|
||||
repetition_number = std::stoi(argv[5]);
|
||||
repetition_number = std::stoi(argv[6]);
|
||||
#endif
|
||||
|
||||
bool ok;
|
||||
for (int cycle = start_cycle; cycle < stop_cycle + 1; cycle++) {
|
||||
ok = test_cycle(file, cycle, version, tolerance);
|
||||
ok = test_cycle(file, cycle, version, breakdown, tolerance);
|
||||
if (ok) {
|
||||
std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;
|
||||
} else {
|
||||
|
||||
Loading…
Reference in New Issue
Block a user