From f096db6411eb55c7d6069b768516963e47693f66 Mon Sep 17 00:00:00 2001 From: Pablo Oliveira Date: Mon, 15 Mar 2021 14:01:38 +0100 Subject: [PATCH] New implementations for Sherman Morisson - SM2: Slagel splitting implementation (http://hdl.handle.net/10919/52966) - SM3: Keep close to zero denominators for later - Update tests to show the squared residual --- Makefile | 8 +-- src/SM_Standard.cpp | 141 +++++++++++++++++++++++++++++++++++-- tests/test_external_h5.cpp | 116 ------------------------------ tests/test_internal_h5.cpp | 21 +++++- 4 files changed, 158 insertions(+), 128 deletions(-) delete mode 100644 tests/test_external_h5.cpp diff --git a/Makefile b/Makefile index 1249960..1de1eb4 100644 --- a/Makefile +++ b/Makefile @@ -1,11 +1,11 @@ ## Compilers -ARCH = -CXX = icpc -FC = ifort +ARCH = +CXX = clang++ +FC = flang-7 H5CXX = h5c++ ## Compiler flags -CXXFLAGS = -O0 -debug full -traceback +CXXFLAGS = -O0 FFLAGS = $(CXXFLAGS) H5CXXFLAGS = $(CXXFLAGS) -fPIC FLIBS = -lstdc++ diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp index ae1e7ff..2d1633f 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -4,15 +4,18 @@ #include "Helpers.hpp" #include -void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + +// Naïve Sherman Morrison +void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { + std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl; double C[Dim]; double D[Dim]; + unsigned int l = 0; // For each update - for (unsigned int l = 0; l < N_updates; l++) { - + while (l < N_updates) { // C = A^{-1} x U_l for (unsigned int i = 0; i < Dim; i++) { C[i] = 0; @@ -23,8 +26,8 @@ void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // Denominator double den = 1 + C[Updates_index[l] - 1]; - if (den < 1e-6) { - std::cerr << "Breakdown condition triggered" << std::endl; + if (fabs(den) < 1e-6) { + std::cerr << "Breakdown condition triggered at " << l << std::endl; } double iden = 1 / den; @@ -40,5 +43,133 @@ void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates, Slater_inv[i * Dim + j] -= update; } } + + 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) { + + std::cerr << "Called SM2 with updates " << N_updates << std::endl; + double C[Dim]; + double D[Dim]; + + double later_updates[Dim*N_updates]; + unsigned int later_index[N_updates]; + unsigned int later = 0; + + 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 (fabs(den) < 1e-6) { + std::cerr << "Breakdown condition triggered at " << l << std::endl; + + // U_l = U_l / 2 (do the split) + for (unsigned int j = 0; j < Dim; j++) { + later_updates[later*Dim+j] = Updates[l*Dim+j]/2.0; + C[j] /= 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) { + SM2(Slater_inv, Dim, later, later_updates, later_index); + } +} + +// 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) { + + std::cerr << "Called SM2 with updates " << N_updates << std::endl; + double C[Dim]; + double D[Dim]; + + double later_updates[Dim*N_updates]; + unsigned int later_index[N_updates]; + unsigned int later = 0; + + 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 (fabs(den) < 1e-6) { + std::cerr << "Breakdown condition triggered at " << l << std::endl; + + for (unsigned int j = 0; j < Dim; j++) { + later_updates[later*Dim+j] = Updates[l*Dim+j]; + } + later_index[later] = Updates_index[l]; + later++; + l += 1; + continue; + } + 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) { + SM3(Slater_inv, Dim, later, later_updates, later_index); + } +} + + + +void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index) { + SM2(Slater_inv, Dim, N_updates, Updates, Updates_index); +} diff --git a/tests/test_external_h5.cpp b/tests/test_external_h5.cpp deleted file mode 100644 index 44c4a5d..0000000 --- a/tests/test_external_h5.cpp +++ /dev/null @@ -1,116 +0,0 @@ -#include -#include -#include "hdf5/serial/hdf5.h" -#include "hdf5/serial/H5Cpp.h" - -#include "SM_MaponiA3.hpp" -#include "SM_Standard.hpp" -#include "Helpers.hpp" - -using namespace H5; -//#define DEBUG - -const H5std_string FILE_NAME( "datasets.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) { - - /* 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); - slater_inverse = transpose(slater_inverse, dim); - - 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 DEBUG - showMatrix(slater_matrix, dim, "OLD Slater"); -#endif - -#ifdef DEBUG - showMatrix(slater_inverse, dim, "OLD Inverse"); -#endif - - 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]; - } - } - - //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index); - SM(slater_inverse, dim, nupdates, u, col_update_index); - -#ifdef DEBUG - showMatrix(slater_matrix, dim, "NEW Slater"); -#endif - -#ifdef DEBUG - 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, 0.5e-4); - -#ifdef DEBUG - showMatrix(res, dim, "Result"); -#endif - - delete [] res, updates, u, col_update_index, - slater_matrix, slater_inverse; - - return ok; -} - -int main(int argc, char **argv) { - if (argc != 2) { - std::cerr << "Execute from within 'datasets/'" << std::endl; - std::cerr << "usage: test_external_h5 " << std::endl; - return 1; - } - int cycle = std::stoi(argv[1]); - H5File file(FILE_NAME, H5F_ACC_RDONLY); - - bool ok = test_cycle(file, cycle); - - if (ok) { - std::cerr << "ok -- cycle " << std::to_string(cycle) - << std::endl; - } - else { - std::cerr << "failed -- cycle " << std::to_string(cycle) - << std::endl; - } - return ok; -} - diff --git a/tests/test_internal_h5.cpp b/tests/test_internal_h5.cpp index 5d55804..edafc8a 100644 --- a/tests/test_internal_h5.cpp +++ b/tests/test_internal_h5.cpp @@ -12,6 +12,18 @@ using namespace H5; const H5std_string FILE_NAME( "datasets.hdf5" ); + +double residual2(double * A, unsigned int Dim) { + double res = 0.0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + double delta = (A[i * Dim + j] - (i == j)); + res += delta*delta; + } + } + return res; +} + void read_int(H5File file, std::string key, unsigned int * data) { DataSet ds = file.openDataSet(key); ds.read(data, PredType::STD_U32LE); @@ -66,8 +78,8 @@ int test_cycle(H5File file, int cycle) { } } - MaponiA3(slater_inverse, dim, nupdates, u, col_update_index); - //SM(slater_inverse, dim, nupdates, updates, col_update_index); + //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index); + SM(slater_inverse, dim, nupdates, u, col_update_index); #ifdef DEBUG showMatrix(slater_matrix, dim, "NEW Slater"); @@ -79,7 +91,10 @@ int test_cycle(H5File file, int cycle) { double * res = new double[dim*dim] {0}; matMul(slater_matrix, slater_inverse, res, dim); - bool ok = is_identity(res, dim, 0.5e-4); + bool ok = is_identity(res, dim, 1e-4); + + double res2 = residual2(res, dim); + std::cerr << "Residual = " << res2 << std::endl; #ifdef DEBUG showMatrix(res, dim, "Result");