diff --git a/include/SM_Standard.hpp b/include/SM_Standard.hpp index e40feda..0332ec7 100644 --- a/include/SM_Standard.hpp +++ b/include/SM_Standard.hpp @@ -1,3 +1,12 @@ -void SM(double *Slater_inv, unsigned int Dim, - unsigned int N_updates, double *Updates, - unsigned int *Updates_index); +// Naïve Sherman Morrison +void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); + +// 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); + +// 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); diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp index 2d1633f..6a7a9d8 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -166,10 +166,3 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, 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_internal_h5.cpp b/tests/test_internal_h5.cpp index edafc8a..4fb2e42 100644 --- a/tests/test_internal_h5.cpp +++ b/tests/test_internal_h5.cpp @@ -36,7 +36,7 @@ void read_double(H5File file, std::string key, double * data) { ds.close(); } -int test_cycle(H5File file, int cycle) { +int test_cycle(H5File file, int cycle, std::string version) { /* Read the data */ @@ -78,8 +78,19 @@ int test_cycle(H5File file, int cycle) { } } - //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index); - SM(slater_inverse, dim, nupdates, u, col_update_index); + if (version == "maponia3") { + MaponiA3(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 { + std::cerr << "Unknown version " << version << std::endl; + exit(1); + } + #ifdef DEBUG showMatrix(slater_matrix, dim, "NEW Slater"); @@ -91,10 +102,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, 1e-4); + bool ok = is_identity(res, dim, 1e-3); double res2 = residual2(res, dim); - std::cerr << "Residual = " << res2 << std::endl; + std::cout << "Residual = " << version << " " << cycle << " " << res2 << std::endl; #ifdef DEBUG showMatrix(res, dim, "Result"); @@ -107,18 +118,19 @@ int test_cycle(H5File file, int cycle) { } int main(int argc, char **argv) { - if (argc != 3) { + if (argc != 4) { std::cerr << "Execute from within 'datasets/'" << std::endl; - std::cerr << "usage: test_internal_h5 " << std::endl; + std::cerr << "usage: test_internal_h5 " << std::endl; return 1; } - int start_cycle = std::stoi(argv[1]); - int stop_cycle = std::stoi(argv[2]); + std::string version(argv[1]); + int start_cycle = std::stoi(argv[2]); + int stop_cycle = std::stoi(argv[3]); H5File file(FILE_NAME, H5F_ACC_RDONLY); bool ok; for (int cycle = start_cycle; cycle < stop_cycle+1; cycle++) { - ok = test_cycle(file, cycle); + ok = test_cycle(file, cycle, version); if (ok) { std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;