mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-11-04 05:03:59 +01:00
- Added tests/fnu_test_h5.cpp that takes its input cycle numbers from a file instead of from the command line.
- Suppressed all debug output so as to not polute performance measurements.
This commit is contained in:
Francois Coppens
parent
99fcaa60b6
commit
43b996dad9
12
Makefile
12
Makefile
@ -2,14 +2,14 @@
|
||||
ifeq ($(ENV),INTEL)
|
||||
CXX = icpc
|
||||
FC = ifort
|
||||
ARCH = -xCORE-AVX2
|
||||
OPT = -O2
|
||||
ARCH = -march=native
|
||||
OPT = -O3
|
||||
DEBUG = -g -debug full
|
||||
else ifeq ($(ENV),LLVM)
|
||||
CXX = clang++
|
||||
FC = flang
|
||||
ARCH =
|
||||
OPT = -O0
|
||||
ARCH = -march=native
|
||||
OPT = -O3
|
||||
DEBUG = -g
|
||||
else ifeq ($(ENV),GNU)
|
||||
CXX = g++
|
||||
@ -58,6 +58,7 @@ BIN_DIR := bin
|
||||
|
||||
EXEC := $(BIN_DIR)/cMaponiA3_test_3x3_3 \
|
||||
$(BIN_DIR)/test_h5 \
|
||||
$(BIN_DIR)/fnu_test_h5 \
|
||||
$(BIN_DIR)/fMaponiA3_test_3x3_3 \
|
||||
$(BIN_DIR)/fMaponiA3_test_4x4_2 \
|
||||
$(BIN_DIR)/QMCChem_dataset_test
|
||||
@ -111,6 +112,9 @@ $(BIN_DIR)/cMaponiA3_test_3x3_3: $(OBJ_DIR)/cMaponiA3_test_3x3_3.o $(DEPS_CXX) |
|
||||
$(BIN_DIR)/test_h5: $(OBJ_DIR)/test_h5.o $(DEPS_CXX) | $(BIN_DIR)
|
||||
$(H5CXX) $(H5LFLAGS) -o $@ $^
|
||||
|
||||
$(BIN_DIR)/fnu_test_h5: $(OBJ_DIR)/fnu_test_h5.o $(DEPS_CXX) | $(BIN_DIR)
|
||||
$(H5CXX) $(H5LFLAGS) -o $@ $^
|
||||
|
||||
$(BIN_DIR)/fMaponiA3_test_3x3_3: $(DEPS_F) $(OBJ_DIR)/fMaponiA3_test_3x3_3.o | $(BIN_DIR)
|
||||
$(FC) $(LFLAGS) $(FLIBS) -o $@ $^
|
||||
|
||||
|
||||
@ -9,10 +9,10 @@
|
||||
#endif
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
#ifndef THRESHOLD
|
||||
#define THRESHOLD 1e-3
|
||||
#endif
|
||||
|
||||
double threshold();
|
||||
|
||||
void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
|
||||
|
||||
118
src/SMWB.cpp
118
src/SMWB.cpp
@ -71,112 +71,6 @@ 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 = &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.
|
||||
// }
|
||||
// }
|
||||
|
||||
// // 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 = &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.
|
||||
// }
|
||||
// }
|
||||
|
||||
// Sherman-Morrison-Woodbury kernel 4
|
||||
// WB2, SM2 mixing scheme
|
||||
void SMWB4(double *Slater_inv, const unsigned int Dim, const unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
|
||||
@ -202,7 +96,9 @@ void SMWB4(double *Slater_inv, const unsigned int Dim, const unsigned int N_upda
|
||||
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 DEBUG1
|
||||
std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2star" << std::endl;
|
||||
#endif
|
||||
#ifdef DEBUG2
|
||||
showMatrix2(Updates_2block, 2, Dim, "Updates_2block");
|
||||
showMatrix2(Updates_index_2block, 1, 2, "Updates_index_2block");
|
||||
@ -233,14 +129,6 @@ 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);
|
||||
|
||||
@ -4,8 +4,6 @@
|
||||
#include "SM_Maponi.hpp"
|
||||
#include "Helpers.hpp"
|
||||
|
||||
// #define DEBUG
|
||||
|
||||
void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
std::cerr << "Called MaponiA3 with " << N_updates << " updates" << std::endl;
|
||||
@ -77,8 +75,10 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
std::cerr << std::endl;
|
||||
#endif
|
||||
if (std::fabs(beta) < threshold()) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
#endif
|
||||
}
|
||||
double ibeta = 1.0 / beta;
|
||||
|
||||
@ -209,8 +209,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
std::cerr << std::endl;
|
||||
#endif
|
||||
if (std::fabs(beta) < threshold()) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
#endif
|
||||
for (unsigned int i = 1; i < Dim + 1; i++) {
|
||||
ylk[l][p[l + 1]][i] *= 0.5;
|
||||
later_updates[later * Dim + i - 1] = Updates[l * Dim + i - 1] * 0.5;
|
||||
|
||||
@ -3,14 +3,14 @@
|
||||
#include "SM_Standard.hpp"
|
||||
#include "Helpers.hpp"
|
||||
|
||||
#define DEBUG1
|
||||
// #define DEBUG1
|
||||
|
||||
// Naïve Sherman Morrison
|
||||
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
#ifdef DEBUG1
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
@ -29,8 +29,10 @@ 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()) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
#endif
|
||||
}
|
||||
double iden = 1 / den;
|
||||
|
||||
@ -55,9 +57,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
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
double later_updates[Dim * N_updates];
|
||||
unsigned int later_index[N_updates];
|
||||
@ -80,11 +82,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
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
std::cerr << "Denominator = " << den << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
// U_l = U_l / 2 (do the split)
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
@ -122,9 +124,9 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
// 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
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM2* with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
@ -143,11 +145,11 @@ 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()) {
|
||||
#ifdef DEBUG1
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
std::cerr << "Denominator = " << den << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
// U_l = U_l / 2 (do the split)
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
@ -181,9 +183,9 @@ 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) {
|
||||
#ifdef DEBUG1
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM3 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
@ -206,9 +208,10 @@ 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()) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
|
||||
#endif
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
later_updates[later * Dim + j] = Updates[l * Dim + j];
|
||||
}
|
||||
@ -236,8 +239,10 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// If all the updates have failed, exit early with an error
|
||||
if (later == N_updates) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "SM3 cannot invert matrix." << std::endl;
|
||||
showMatrix(Slater_inv, Dim, "Slater_inverse");
|
||||
#endif
|
||||
return;
|
||||
}
|
||||
// If some have failed, make a recursive call
|
||||
@ -251,9 +256,9 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
// (SM2)
|
||||
void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
#ifdef DEBUG1
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Called SM4 with " << N_updates << " updates" << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
@ -276,9 +281,10 @@ 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()) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
|
||||
#endif
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
later_updates[later * Dim + j] = Updates[l * Dim + j];
|
||||
}
|
||||
|
||||
@ -51,10 +51,10 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
|
||||
// Check if determinant of inverted matrix is not zero
|
||||
double det = B[0] * B[3] - B[1] * B[2];
|
||||
if (std::fabs(det) < threshold()) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Determinant = " << det << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
return false;
|
||||
}
|
||||
|
||||
@ -160,11 +160,10 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
|
||||
std::cerr << "Determinant of B = " << det << std::endl;
|
||||
#endif
|
||||
if (std::fabs(det) < threshold()) {
|
||||
// if (std::fabs(det) < 1000000) {
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "Determinant = " << det << std::endl;
|
||||
#endif
|
||||
#endif
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
251
tests/fnu_test_h5.cpp
Normal file
251
tests/fnu_test_h5.cpp
Normal file
@ -0,0 +1,251 @@
|
||||
#include "hdf5/serial/H5Cpp.h"
|
||||
#include "hdf5/serial/hdf5.h"
|
||||
|
||||
#include "Helpers.hpp"
|
||||
#include "SM_Maponi.hpp"
|
||||
#include "SM_Standard.hpp"
|
||||
#include "Woodbury.hpp"
|
||||
#include "SMWB.hpp"
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
|
||||
#define PERF
|
||||
// #define STATUS
|
||||
// #define RESIDUAL
|
||||
|
||||
#ifdef PERF
|
||||
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);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
void read_double(H5File file, std::string key, double *data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::IEEE_F64LE);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
int test_cycle(H5File file, int cycle, std::string version, 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);
|
||||
|
||||
double *slater_matrix = new double[dim * dim];
|
||||
read_double(file, group + "/slater_matrix", slater_matrix);
|
||||
|
||||
double *slater_inverse = new double[dim * dim];
|
||||
read_double(file, group + "/slater_inverse", slater_inverse);
|
||||
|
||||
unsigned int *col_update_index = new unsigned int[nupdates];
|
||||
read_int(file, group + "/col_update_index", col_update_index);
|
||||
|
||||
double *updates = new double[nupdates * dim];
|
||||
read_double(file, group + "/updates", updates);
|
||||
|
||||
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++) {
|
||||
for (i = 0; i < dim; i++) {
|
||||
col = col_update_index[j];
|
||||
u[i + j * dim] =
|
||||
updates[i + j * dim] - slater_matrix[i * dim + (col - 1)];
|
||||
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
|
||||
|
||||
#ifdef PERF
|
||||
#ifdef DEBUG1
|
||||
std::cerr << "# of reps. = " << repetition_number << std::endl;
|
||||
#endif
|
||||
double *slater_inverse_nonpersistent = new double[dim * dim];
|
||||
if (version == "sm1") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
SM1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
}
|
||||
} else if (version == "sm2") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
SM2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
}
|
||||
} else if (version == "sm3") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
SM3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
}
|
||||
} else if (version == "sm4") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
SM4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
}
|
||||
} else if (version == "wb2") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
WB2(slater_inverse_nonpersistent, dim, u, col_update_index);
|
||||
}
|
||||
} else if (version == "wb3") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
WB3(slater_inverse_nonpersistent, dim, u, col_update_index);
|
||||
}
|
||||
} else if (version == "smwb1") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
}
|
||||
} else if (version == "smwb4") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
|
||||
SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
|
||||
}
|
||||
#ifdef MKL
|
||||
} else if (version == "lapack") {
|
||||
for (unsigned int i = 0; i < repetition_number; i++) {
|
||||
std::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;
|
||||
#else
|
||||
if (version == "maponia3") {
|
||||
MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "maponia3s") {
|
||||
MaponiA3S(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm1") {
|
||||
SM1(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm2") {
|
||||
SM2(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm3") {
|
||||
SM3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "sm4") {
|
||||
SM4(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
} else if (version == "wb2") {
|
||||
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 == "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));
|
||||
inverse(slater_inverse, dim);
|
||||
#endif // MKL
|
||||
} else {
|
||||
std::cerr << "Unknown version " << version << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
#endif // PERF
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(slater_matrix, dim, "NEW Slater");
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(slater_inverse, dim, "NEW Inverse");
|
||||
#endif
|
||||
|
||||
double *res = new double[dim * dim]{0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
bool ok = is_identity(res, dim, tolerance);
|
||||
double res_max = residual_max(res, dim);
|
||||
double res2 = residual_frobenius2(res, dim);
|
||||
|
||||
#ifdef RESIDUAL
|
||||
std::cout << "Residual = " << version << " " << cycle << " " << res_max << " "
|
||||
<< res2 << std::endl;
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG2
|
||||
showMatrix(res, dim, "Result");
|
||||
#endif
|
||||
|
||||
delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse;
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
#ifdef PERF
|
||||
if (argc != 5) {
|
||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||
std::cerr
|
||||
<< "usage: test_h5 <version> <cycle file> <tolerance> <number of reps.>"
|
||||
<< std::endl;
|
||||
return 1;
|
||||
}
|
||||
#else
|
||||
if (argc != 4) {
|
||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||
std::cerr
|
||||
<< "usage: test_h5 <version> <cycle file> <tolerance>"
|
||||
<< std::endl;
|
||||
return 1;
|
||||
}
|
||||
#endif
|
||||
std::string version(argv[1]);
|
||||
std::string cyclefile_name(argv[2]);
|
||||
std::ifstream cyclefile(cyclefile_name);
|
||||
std::vector<int> cycles;
|
||||
unsigned int cycle;
|
||||
while (cyclefile >> cycle) cycles.push_back(cycle);
|
||||
double tolerance = std::stod(argv[3]);
|
||||
H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
||||
|
||||
#ifdef PERF
|
||||
repetition_number = std::stoi(argv[4]);
|
||||
#endif
|
||||
|
||||
bool ok;
|
||||
for (auto & cycle : cycles) {
|
||||
ok = test_cycle(file, cycle, version, tolerance);
|
||||
#ifdef STATUS
|
||||
if (ok) {
|
||||
std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;
|
||||
} else {
|
||||
std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
return ok;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user