TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2024-12-25 22:03:57 +01:00
Code Issues Projects Releases Wiki Activity

Merge pull request #51 from fmgjcoppens:qmckl_integration

Browse Source 
qmckl_integration
This commit is contained in:
François Coppens 2021-07-19 08:17:01 +02:00 committed by GitHub
commit 5ca2ac4d5c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 49 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Makefile
View File

@ -1,6 +1,6 @@
## Compilers, compiler flags & external libs ## Compilers, compiler flags & external libs
ifeq ($(ENV),INTEL) ifeq ($(ENV),INTEL)
CXX = icpc CXX = icpx
FC = ifort FC = ifort
ARCH = -march=native ARCH = -march=native
OPT = -O3 OPT = -O3
@ -15,7 +15,7 @@ else ifeq ($(ENV),GNU)
CXX = g++ CXX = g++
FC = gfortran FC = gfortran
ARCH = -mavx ARCH = -mavx
OPT = -O0 OPT = -O3
DEBUG = -g DEBUG = -g
else else
$(error No valid compiler environment set in $$ENV. \ $(error No valid compiler environment set in $$ENV. \

14
include/SMWB.hpp
View File

@ -4,18 +4,8 @@ void SMWB1(double *Slater_inv, const unsigned int Dim,
const unsigned int N_updates, double *Updates, const unsigned int N_updates, double *Updates,
unsigned int *Updates_index); unsigned int *Updates_index);
// // Sherman-Morrison-Woodbury kernel 2 // 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 // WB2, SM2 mixing scheme
void SMWB4(double *Slater_inv, const unsigned int Dim, void SMWB2(double *Slater_inv, const unsigned int Dim,
const unsigned int N_updates, double *Updates, const unsigned int N_updates, double *Updates,
unsigned int *Updates_index); unsigned int *Updates_index);

8
smvars.sh
View File

@ -27,20 +27,23 @@ export SMROOT
case $ENV in case $ENV in
intel) intel)
echo "* SM build environment set to 'intel'" echo "* SM build environment set to 'intel'"
export HDF5_CXX=icpc export HDF5_CXX=icpx
export HDF5_CXXLINKER=icpc export HDF5_CXXLINKER=icpx
export HDF5_CLINKER=icpx
export ENV=INTEL export ENV=INTEL
;; ;;
llvm) llvm)
echo "* SM build environment set to 'llvm'" echo "* SM build environment set to 'llvm'"
export HDF5_CXX=clang++ export HDF5_CXX=clang++
export HDF5_CXXLINKER=clang++ export HDF5_CXXLINKER=clang++
export HDF5_CLINKER=clang++
export ENV=LLVM export ENV=LLVM
;; ;;
vfc) vfc)
echo "* SM build environment set to 'vfc'" echo "* SM build environment set to 'vfc'"
export HDF5_CXX=clang++ export HDF5_CXX=clang++
export HDF5_CXXLINKER=clang++ export HDF5_CXXLINKER=clang++
export HDF5_CLINKER=clang++
export ENV=LLVM export ENV=LLVM
export VFC_BACKENDS="libinterflop_ieee.so --count-op" export VFC_BACKENDS="libinterflop_ieee.so --count-op"
;; ;;
@ -48,6 +51,7 @@ case $ENV in
echo "* SM build environment set to 'gnu'" echo "* SM build environment set to 'gnu'"
export HDF5_CXX=g++ export HDF5_CXX=g++
export HDF5_CXXLINKER=g++ export HDF5_CXXLINKER=g++
export HDF5_CLINKER=g++
export ENV=GNU export ENV=GNU
;; ;;
*) *)

8
src/SMWB.cpp
View File

@ -80,9 +80,9 @@ void SMWB1(double *Slater_inv, const unsigned int Dim,
} }
} }
// Sherman-Morrison-Woodbury kernel 4 // Sherman-Morrison-Woodbury kernel 2
// WB2, SM2 mixing scheme // WB2, SM2 mixing scheme
void SMWB4(double *Slater_inv, const unsigned int Dim, void SMWB2(double *Slater_inv, const unsigned int Dim,
const unsigned int N_updates, double *Updates, const unsigned int N_updates, double *Updates,
unsigned int *Updates_index) { unsigned int *Updates_index) {
#ifdef DEBUG2 #ifdef DEBUG2
@ -144,8 +144,8 @@ void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
double **linUpdates, unsigned int **Updates_index) { double **linUpdates, unsigned int **Updates_index) {
SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
} }
void SMWB4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, void SMWB2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
double **linUpdates, unsigned int **Updates_index) { double **linUpdates, unsigned int **Updates_index) {
SMWB4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); SMWB2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
} }
} }

57
src/SM_Standard.cpp
View File

@ -65,55 +65,8 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int later_index[N_updates]; unsigned int later_index[N_updates];
unsigned int later = 0; unsigned int later = 0;
double C[Dim]; SM2star(Slater_inv, Dim, N_updates, Updates, Updates_index, later_updates,
double D[Dim]; later_index, &later);
unsigned int l = 0;
// For each update
while (l < N_updates) {
// C = A^{-1} x U_l
for (unsigned int i = 0; i < Dim; i++) {
C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) {
C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
}
}
// Denominator
double den = 1 + C[Updates_index[l] - 1];
if (std::fabs(den) < threshold()) {
#ifdef DEBUG1
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
<< std::endl;
std::cerr << "Denominator = " << den << std::endl;
#endif
// U_l = U_l / 2 (do the split)
for (unsigned int i = 0; i < Dim; i++) {
later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
C[i] /= 2.0;
}
later_index[later] = Updates_index[l];
later++;
den = 1 + C[Updates_index[l] - 1];
}
double iden = 1 / den;
// D = v^T x A^{-1}
for (unsigned int j = 0; j < Dim; j++) {
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
}
// A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update;
}
}
l += 1;
}
if (later > 0) { if (later > 0) {
SM2(Slater_inv, Dim, later, later_updates, later_index); SM2(Slater_inv, Dim, later, later_updates, later_index);
@ -136,7 +89,7 @@ void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) { while (l < N_updates) {
// C = A^{-1} x U_l // C = S^{-1} x U_l
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++) {
C[i] = 0; C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++) {
@ -165,12 +118,12 @@ void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
} }
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // D = v^T x S^{-1}
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++) {
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // S^{-1} = S^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++) {
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;

22
src/Woodbury.cpp
View File

@ -134,8 +134,8 @@ bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
// Check if determinant of B is not too close to zero // Check if determinant of B is not too close to zero
double det; double det;
det = B0 * (B4 * B8 - B5 * B7) - det = B0 * (B4 * B8 - B5 * B7) - B1 * (B3 * B8 - B5 * B6) +
B1 * (B3 * B8 - B5 * B6) + B2 * (B3 * B7 - B4 * B6); B2 * (B3 * B7 - B4 * B6);
#ifdef DEBUG2 #ifdef DEBUG2
std::cerr << "Determinant of B = " << det << std::endl; std::cerr << "Determinant of B = " << det << std::endl;
#endif #endif
@ -150,15 +150,15 @@ bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
// Compute B^{-1} with explicit formula for 3x3 inversion // Compute B^{-1} with explicit formula for 3x3 inversion
double Binv[9], idet = 1.0 / det; double Binv[9], idet = 1.0 / det;
Binv[0] = (B4 * B8 - B7 * B5) * idet; Binv[0] = (B4 * B8 - B7 * B5) * idet;
Binv[1] = -(B1 * B8 - B7 * B2) * idet; Binv[1] = -(B1 * B8 - B7 * B2) * idet;
Binv[2] = (B1 * B5 - B4 * B2) * idet; Binv[2] = (B1 * B5 - B4 * B2) * idet;
Binv[3] = -(B3 * B8 - B6 * B5) * idet; Binv[3] = -(B3 * B8 - B6 * B5) * idet;
Binv[4] = (B0 * B8 - B6 * B2) * idet; Binv[4] = (B0 * B8 - B6 * B2) * idet;
Binv[5] = -(B0 * B5 - B3 * B2) * idet; Binv[5] = -(B0 * B5 - B3 * B2) * idet;
Binv[6] = (B3 * B7 - B6 * B4) * idet; Binv[6] = (B3 * B7 - B6 * B4) * idet;
Binv[7] = -(B0 * B7 - B6 * B1) * idet; Binv[7] = -(B0 * B7 - B6 * B1) * idet;
Binv[8] = (B0 * B4 - B3 * B1) * idet; Binv[8] = (B0 * B4 - B3 * B1) * idet;
#ifdef DEBUG2 #ifdef DEBUG2
std::cerr << "Conditioning number of B = " << condition1(B, Binv, 3) std::cerr << "Conditioning number of B = " << condition1(B, Binv, 3)
@ -198,10 +198,14 @@ bool WB3(double *Slater_inv, const unsigned int Dim, double *Updates,
extern "C" { extern "C" {
bool WB2_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates, bool WB2_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
unsigned int **Updates_index) { unsigned int **Updates_index) {
WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index); bool ok;
ok = WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
return ok;
} }
bool WB3_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates, bool WB3_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
unsigned int **Updates_index) { unsigned int **Updates_index) {
WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index); bool ok;
ok = WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
return ok;
} }
} }

18
tests/fnu_test_h5.cpp
View File

@ -10,7 +10,7 @@
#include <vector> #include <vector>
#define PERF #define PERF
// #define STATUS #define STATUS
// #define RESIDUAL // #define RESIDUAL
#ifdef PERF #ifdef PERF
@ -126,11 +126,11 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
dim * dim * sizeof(double)); dim * dim * sizeof(double));
SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
} }
} else if (version == "smwb4") { } else if (version == "smwb2") {
for (unsigned int i = 0; i < repetition_number; i++) { for (unsigned int i = 0; i < repetition_number; i++) {
std::memcpy(slater_inverse_nonpersistent, slater_inverse, std::memcpy(slater_inverse_nonpersistent, slater_inverse,
dim * dim * sizeof(double)); dim * dim * sizeof(double));
SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
} }
#ifdef MKL #ifdef MKL
} else if (version == "lapack") { } else if (version == "lapack") {
@ -166,12 +166,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
WB3(slater_inverse, dim, u, col_update_index); WB3(slater_inverse, dim, u, col_update_index);
} else if (version == "smwb1") { } else if (version == "smwb1") {
SMWB1(slater_inverse, dim, nupdates, u, col_update_index); SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
// } else if (version == "smwb2") { } else if (version == "smwb2") {
// SMWB2(slater_inverse, dim, nupdates, u, col_update_index); SMWB2(slater_inverse, dim, nupdates, u, col_update_index);
// } else if (version == "smwb3") {
// SMWB3(slater_inverse, dim, nupdates, u, col_update_index);
} else if (version == "smwb4") {
SMWB4(slater_inverse, dim, nupdates, u, col_update_index);
#ifdef MKL #ifdef MKL
} else if (version == "lapack") { } else if (version == "lapack") {
memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double)); memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double));
@ -213,14 +209,14 @@ int main(int argc, char **argv) {
if (argc != 5) { if (argc != 5) {
std::cerr << "Execute from within 'datasets/'" << std::endl; std::cerr << "Execute from within 'datasets/'" << std::endl;
std::cerr std::cerr
<< "usage: test_h5 <version> <cycle file> <tolerance> <number of reps.>" << "usage: fnu_test_h5 <version> <cycle file> <tolerance> <number of reps.>"
<< std::endl; << std::endl;
return 1; return 1;
} }
#else #else
if (argc != 4) { if (argc != 4) {
std::cerr << "Execute from within 'datasets/'" << std::endl; std::cerr << "Execute from within 'datasets/'" << std::endl;
std::cerr << "usage: test_h5 <version> <cycle file> <tolerance>" std::cerr << "usage: fnu_test_h5 <version> <cycle file> <tolerance>"
<< std::endl; << std::endl;
return 1; return 1;
} }

23
tests/test_h5.cpp
View File

@ -104,14 +104,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
WB3(slater_inverse_nonpersistent, dim, u, col_update_index); WB3(slater_inverse_nonpersistent, dim, u, col_update_index);
} else if (version == "smwb1") { } else if (version == "smwb1") {
SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index); SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
// } else if (version == "smwb2") { } else if (version == "smwb2") {
// SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
// col_update_index);
// } else if (version == "smwb3") {
// SMWB3(slater_inverse_nonpersistent, dim, nupdates, u,
// col_update_index);
} else if (version == "smwb4") {
SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
#ifdef MKL #ifdef MKL
} else if (version == "lapack") { } else if (version == "lapack") {
memcpy(slater_inverse_nonpersistent, slater_matrix, memcpy(slater_inverse_nonpersistent, slater_matrix,
@ -145,10 +139,6 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
WB3(slater_inverse, dim, u, col_update_index); WB3(slater_inverse, dim, u, col_update_index);
} else if (version == "smwb1") { } else if (version == "smwb1") {
SMWB1(slater_inverse, dim, nupdates, u, col_update_index); SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
// } else if (version == "smwb2") {
// 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") { } else if (version == "smwb4") {
SMWB4(slater_inverse, dim, nupdates, u, col_update_index); SMWB4(slater_inverse, dim, nupdates, u, col_update_index);
#ifdef MKL #ifdef MKL
@ -191,15 +181,22 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
int main(int argc, char **argv) { int main(int argc, char **argv) {
#ifdef PERF #ifdef PERF
if (argc != 6) { if (argc != 6) {
std::cerr << "Execute from within 'datasets/'" << std::endl;
std::cerr
<< "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance> <number of reps.>"
<< std::endl;
return 1;
}
#else #else
if (argc != 5) { if (argc != 5) {
#endif
std::cerr << "Execute from within 'datasets/'" << std::endl; std::cerr << "Execute from within 'datasets/'" << std::endl;
std::cerr std::cerr
<< "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>" << "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>"
<< std::endl; << std::endl;
return 1; return 1;
} }
#endif
std::string version(argv[1]); std::string version(argv[1]);
int start_cycle = std::stoi(argv[2]); int start_cycle = std::stoi(argv[2]);
int stop_cycle = std::stoi(argv[3]); int stop_cycle = std::stoi(argv[3]);