From adf893ed6f359fb8595749833c7e8ed77fa23a86 Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Fri, 7 May 2021 17:11:04 +0200 Subject: [PATCH] - Various fixes - Merged MaponiA3 and MaponiA3S into one file. --- Makefile | 7 +- Makefile.verificarlo | 7 +- include/SM_Helpers.hpp | 50 +++++++ include/SM_Maponi.hpp | 8 ++ include/SM_MaponiA3.hpp | 2 - include/SM_MaponiA3S.hpp | 2 - smvars.sh | 1 + src/{SM_MaponiA3S.cpp => SM_Maponi.cpp} | 163 ++++++++++++++++++++--- src/SM_MaponiA3.cpp | 142 -------------------- src/SM_Standard.cpp | 7 +- src/{SM_mod.f90 => finterface_mod.f90} | 0 src/{Helpers_mod.f90 => helpers_mod.f90} | 0 tests/cMaponiA3_test_3x3_3.cpp | 2 +- tests/run_test.sh | 3 +- tests/test_h5.cpp | 32 +---- tests/vfc_test_h5.cpp | 28 +--- 16 files changed, 220 insertions(+), 234 deletions(-) create mode 100644 include/SM_Maponi.hpp delete mode 100644 include/SM_MaponiA3.hpp delete mode 100644 include/SM_MaponiA3S.hpp rename src/{SM_MaponiA3S.cpp => SM_Maponi.cpp} (51%) delete mode 100644 src/SM_MaponiA3.cpp rename src/{SM_mod.f90 => finterface_mod.f90} (100%) rename src/{Helpers_mod.f90 => helpers_mod.f90} (100%) diff --git a/Makefile b/Makefile index 439b5e1..f144ad7 100644 --- a/Makefile +++ b/Makefile @@ -29,13 +29,12 @@ FLIBS = -lstdc++ ## Includes and dependencies INCLUDE = -I $(INC_DIR)/ -DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \ - $(OBJ_DIR)/SM_MaponiA3S.o \ +DEPS_CXX = $(OBJ_DIR)/SM_Maponi.o \ $(OBJ_DIR)/SM_Standard.o \ $(OBJ_DIR)/SM_Helpers.o DEPS_F = $(DEPS_CXX) \ - $(OBJ_DIR)/SM_mod.o \ - $(OBJ_DIR)/Helpers_mod.o + $(OBJ_DIR)/finterface_mod.o \ + $(OBJ_DIR)/helpers_mod.o ## Directory structure SRC_DIR := src diff --git a/Makefile.verificarlo b/Makefile.verificarlo index 7e5c347..ec35224 100644 --- a/Makefile.verificarlo +++ b/Makefile.verificarlo @@ -16,13 +16,12 @@ CXXFLAGS = -O0 -g $(H5FLAGS) FFLAGS = $(CXXFLAGS) INCLUDE = -I $(INC_DIR)/ -DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \ - $(OBJ_DIR)/SM_MaponiA3S.o \ +DEPS_CXX = $(OBJ_DIR)/SM_Maponi.o \ $(OBJ_DIR)/SM_Standard.o \ $(OBJ_DIR)/SM_Helpers.o DEPS_F = $(DEPS_CXX) \ - $(OBJ_DIR)/SM_mod.o \ - $(OBJ_DIR)/Helpers_mod.o + $(OBJ_DIR)/finterface_mod.o \ + $(OBJ_DIR)/helpers_mod.o FLIBS = -lstdc++ SRC_DIR := src diff --git a/include/SM_Helpers.hpp b/include/SM_Helpers.hpp index ce86d02..73748b9 100644 --- a/include/SM_Helpers.hpp +++ b/include/SM_Helpers.hpp @@ -127,3 +127,53 @@ template bool is_identity(T *A, unsigned int M, double tolerance) { } return true; } + +template T norm_max(T * A, unsigned int Dim) { + T res = 0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + T delta = A[i * Dim + j]; + delta = fabs(delta); + if (delta > res) { + res = delta; + } + } + } + return res; +} + +template T norm_frobenius2(T * A, unsigned int Dim) { + T res = 0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + T delta = A[i * Dim + j]; + res += delta*delta; + } + } + return res; +} + +template T residual_max(T * A, unsigned int Dim) { + T res = 0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + T delta = A[i * Dim + j] - (i == j); + delta = fabs(delta); + if (delta > res) { + res = delta; + } + } + } + return res; +} + +template T residual_frobenius2(T * A, unsigned int Dim) { + T res = 0; + for (unsigned int i = 0; i < Dim; i++) { + for (unsigned int j = 0; j < Dim; j++) { + T delta = A[i * Dim + j] - (i == j); + res += delta*delta; + } + } + return res; +} diff --git a/include/SM_Maponi.hpp b/include/SM_Maponi.hpp new file mode 100644 index 0000000..873cd45 --- /dev/null +++ b/include/SM_Maponi.hpp @@ -0,0 +1,8 @@ +// P. Maponi Algorithm 3 +void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); + +// P. Maponi Algorithm 3 with J. Slagel splitting +// http://hdl.handle.net/10919/52966 +void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); diff --git a/include/SM_MaponiA3.hpp b/include/SM_MaponiA3.hpp deleted file mode 100644 index c4c566b..0000000 --- a/include/SM_MaponiA3.hpp +++ /dev/null @@ -1,2 +0,0 @@ -void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index); diff --git a/include/SM_MaponiA3S.hpp b/include/SM_MaponiA3S.hpp deleted file mode 100644 index 9085b7e..0000000 --- a/include/SM_MaponiA3S.hpp +++ /dev/null @@ -1,2 +0,0 @@ -void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index); diff --git a/smvars.sh b/smvars.sh index caf11d4..db22e5a 100644 --- a/smvars.sh +++ b/smvars.sh @@ -45,6 +45,7 @@ case $ENV in *) echo "Unknown environment descriptor given." echo "Usage: source smvars.sh {intel | llvm | vfc | gnu}" + echo "Usage: source smvars.sh {intel | llvm | vfc | gnu} [threshold]" return 1 ;; esac diff --git a/src/SM_MaponiA3S.cpp b/src/SM_Maponi.cpp similarity index 51% rename from src/SM_MaponiA3S.cpp rename to src/SM_Maponi.cpp index 0b8c4f3..7b8263e 100644 --- a/src/SM_MaponiA3S.cpp +++ b/src/SM_Maponi.cpp @@ -1,17 +1,146 @@ -// SM-MaponiA3_f.cpp +// SM_Maponi.cpp // Algorithm 3 from P. Maponi, // p. 283, doi:10.1016/j.laa.2006.07.007 -#include "SM_MaponiA3S.hpp" +#include "SM_Maponi.hpp" #include "SM_Helpers.hpp" // #define DEBUG -void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index) { +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; + + /* + DECLARE AND INITIALISE ARRAYS + */ + unsigned int k, l, i, j, component; + unsigned int *p = new unsigned int[N_updates + 1]{0}; + double alpha, beta; + double *Al = new double[Dim * Dim]; + double *next = new double[Dim * Dim]{0}; + double *last = Slater_inv, *tmp; + + // Populate update-order vector + for (i = 0; i < N_updates; i++) { + p[i + 1] = i + 1; + } + + // Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim] + double ***ylk = new double **[N_updates]; + for (l = 0; l < N_updates; l++) { + ylk[l] = new double *[N_updates + 1]; + for (k = 0; k < N_updates + 1; k++) { + ylk[l][k] = new double[Dim + 1]{0}; + } + } + + /* + START ALGORITHM + */ + + // Calculate the {y_{0,k}} + for (k = 1; k < N_updates + 1; k++) { +#ifdef DEBUG + std::cerr << "Compute y0k: " << std::endl; + std::cerr << "ylk[0][" << k << "][:]"; + std::cerr << std::endl; +#endif + for (i = 1; i < Dim + 1; i++) { + for (j = 1; j < Dim + 1; j++) { + ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] * + Updates[(k - 1) * Dim + (j - 1)]; + } + } +#ifdef DEBUG + showVector(ylk[0][k], Dim, ""); +#endif + } + + // Calculate the {y_{l,k}} from the {y_{0,k}} + for (l = 0; l < N_updates; l++) { +#ifdef DEBUG + std::cerr << "In outer compute-ylk-loop: l = " << l << std::endl; + std::cerr << std::endl; +#endif + + // For given l select intermediate update with largest break-down val + selectLargestDenominator(l, N_updates, Updates_index, p, ylk); + + // Select component and comp. bd-condition. + component = Updates_index[p[l + 1] - 1]; + beta = 1 + ylk[l][p[l + 1]][component]; +#ifdef DEBUG + std::cerr << "p[l+1] = " << p[l + 1] << std::endl; + std::cerr << "component = " << component << std::endl; + std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component + << "] = " << beta << std::endl; + std::cerr << std::endl; +#endif + if (fabs(beta) < threshold()) { + std::cerr << "Breakdown condition triggered at " << Updates_index[l] + << std::endl; + } + double ibeta = 1.0 / beta; + +// Compute intermediate update to Slater_inv +#ifdef DEBUG + std::cerr << "Compute intermediate update to Slater_inv" << std::endl; + std::cerr << "component = " << component << std::endl; + std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component + << "]" << std::endl; + std::cerr << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]" + << std::endl; + std::cerr << std::endl; +#endif + for (i = 0; i < Dim; i++) { + for (j = 0; j < Dim; j++) { + Al[i * Dim + j] = + (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta; + } + } + matMul(Al, last, next, Dim); + tmp = next; + next = last; + last = tmp; +#ifdef DEBUG + showMatrix(last, Dim, "last"); +#endif + + // For given l != 0 compute the next {y_{l,k}} + for (k = l + 2; k < N_updates + 1; k++) { + alpha = ylk[l][p[k]][component] * ibeta; +#ifdef DEBUG + std::cerr << "Inside k-loop: k = " << k << std::endl; + std::cerr << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl; + std::cerr << std::endl; +#endif + for (i = 1; i < Dim + 1; i++) { + ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i]; + } + } + } + memcpy(Slater_inv, last, Dim * Dim * sizeof(double)); + + /* + CLEANUP MEMORY + */ + + for (l = 0; l < N_updates; l++) { + for (k = 0; k < N_updates + 1; k++) { + delete[] ylk[l][k]; + } + delete[] ylk[l]; + } + delete[] Al, next, p, ylk; +} + +void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index) { + std::cerr << "Called MaponiA3S with " << N_updates << " updates" << std::endl; + /* DECLARE AND INITIALISE ARRAYS */ - showVector(Updates_index, Dim, "Updates_index"); unsigned int k, l, i, j, component; unsigned int *p = new unsigned int[N_updates + 1]{0}; double alpha, beta; @@ -40,7 +169,6 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, /* START ALGORITHM */ - showMatrix(Slater_inv, Dim, "S^{-1}"); // Calculate the y_{0,k} = S_0^{-1} u_k for (k = 1; k < N_updates + 1; k++) { @@ -81,27 +209,17 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, std::cerr << std::endl; #endif if (fabs(beta) < threshold()) { -#ifdef DEBUG - std::cerr << "Breakdown condition triggered at l = " << l << " and component "<< component + std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; -#endif for (unsigned int i = 1; i < Dim + 1; i++) { - std::cout << "Let's split..." << std::endl; ylk[l][p[l + 1]][i] *= 0.5; - // later_updates[later * Dim + i - 1] = ylk[l][p[l + 1]][i]; // Updates[l * Dim + i] * 0.5; later_updates[later * Dim + i - 1] = Updates[l * Dim + i - 1] * 0.5; } later_index[later] = Updates_index[p[l + 1] - 1]; - std::cout << "Later = " << later << std::endl; later++; - std::cout << "Later = " << later << std::endl; beta = 1 + ylk[l][p[l + 1]][component]; - std::cout << "New beta is: " << beta << std::endl; } double ibeta = 1.0 / beta; - showVector(ylk[l][p[l + 1]], Dim+1, ""); - showVector(later_updates+(Dim*0), Dim, "later_updateds"); - // Compute intermediate update to Slater_inv #ifdef DEBUG @@ -135,7 +253,6 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, for (i = 1; i < Dim + 1; i++) { ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i]; } - showVector(ylk[l + 1][p[k]], Dim+1, ""); } } memcpy(Slater_inv, last, Dim * Dim * sizeof(double)); @@ -153,12 +270,18 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, delete[] Al, next, p, ylk; if (later > 0) { - std::cerr << "Entering recursive loop with " << later << " updates" << std::endl; - std::cout << later_index[0] << std::endl; MaponiA3S(Slater_inv, Dim, later, later_updates, later_index); } } +extern "C" { + void MaponiA3_f(double **linSlater_inv, unsigned int *Dim, + unsigned int *N_updates, double **linUpdates, + unsigned int **Updates_index) { + MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); + } +} + extern "C" { void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, double **linUpdates, diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp deleted file mode 100644 index f16f625..0000000 --- a/src/SM_MaponiA3.cpp +++ /dev/null @@ -1,142 +0,0 @@ -// SM-MaponiA3_f.cpp -// Algorithm 3 from P. Maponi, -// p. 283, doi:10.1016/j.laa.2006.07.007 -#include "SM_MaponiA3.hpp" -#include "SM_Helpers.hpp" - -// #define DEBUG - -void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index) { - - /* - DECLARE AND INITIALISE ARRAYS - */ - - unsigned int k, l, i, j, component; - unsigned int *p = new unsigned int[N_updates + 1]{0}; - double alpha, beta; - double *Al = new double[Dim * Dim]; - double *next = new double[Dim * Dim]{0}; - double *last = Slater_inv, *tmp; - - // Populate update-order vector - for (i = 0; i < N_updates; i++) { - p[i + 1] = i + 1; - } - - // Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim] - double ***ylk = new double **[N_updates]; - for (l = 0; l < N_updates; l++) { - ylk[l] = new double *[N_updates + 1]; - for (k = 0; k < N_updates + 1; k++) { - ylk[l][k] = new double[Dim + 1]{0}; - } - } - - /* - START ALGORITHM - */ - - // Calculate the {y_{0,k}} - for (k = 1; k < N_updates + 1; k++) { -#ifdef DEBUG - std::cout << "Compute y0k: " << std::endl; - std::cout << "ylk[0][" << k << "][:]"; - std::cout << std::endl; -#endif - for (i = 1; i < Dim + 1; i++) { - for (j = 1; j < Dim + 1; j++) { - ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] * - Updates[(k - 1) * Dim + (j - 1)]; - } - } -#ifdef DEBUG - showVector(ylk[0][k], Dim, ""); -#endif - } - - // Calculate the {y_{l,k}} from the {y_{0,k}} - for (l = 0; l < N_updates; l++) { -#ifdef DEBUG - std::cout << "In outer compute-ylk-loop: l = " << l << std::endl; - std::cout << std::endl; -#endif - - // For given l select intermediate update with largest break-down val - selectLargestDenominator(l, N_updates, Updates_index, p, ylk); - - // Select component and comp. bd-condition. - component = Updates_index[p[l + 1] - 1]; - beta = 1 + ylk[l][p[l + 1]][component]; -#ifdef DEBUG - std::cout << "p[l+1] = " << p[l + 1] << std::endl; - std::cout << "component = " << component << std::endl; - std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component - << "] = " << beta << std::endl; - std::cout << std::endl; -#endif - if (fabs(beta) < threshold()) { - std::cerr << "Break-down occured." << std::endl; - } - double ibeta = 1.0 / beta; - -// Compute intermediate update to Slater_inv -#ifdef DEBUG - std::cout << "Compute intermediate update to Slater_inv" << std::endl; - std::cout << "component = " << component << std::endl; - std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component - << "]" << std::endl; - std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]" - << std::endl; - std::cout << std::endl; -#endif - for (i = 0; i < Dim; i++) { - for (j = 0; j < Dim; j++) { - Al[i * Dim + j] = - (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta; - } - } - matMul(Al, last, next, Dim); - tmp = next; - next = last; - last = tmp; -#ifdef DEBUG - showMatrix(last, Dim, "last"); -#endif - - // For given l != 0 compute the next {y_{l,k}} - for (k = l + 2; k < N_updates + 1; k++) { - alpha = ylk[l][p[k]][component] * ibeta; -#ifdef DEBUG - std::cout << "Inside k-loop: k = " << k << std::endl; - std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl; - std::cout << std::endl; -#endif - for (i = 1; i < Dim + 1; i++) { - ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i]; - } - } - } - memcpy(Slater_inv, last, Dim * Dim * sizeof(double)); - - /* - CLEANUP MEMORY - */ - - for (l = 0; l < N_updates; l++) { - for (k = 0; k < N_updates + 1; k++) { - delete[] ylk[l][k]; - } - delete[] ylk[l]; - } - delete[] Al, next, p, ylk; -} - -extern "C" { - void MaponiA3_f(double **linSlater_inv, unsigned int *Dim, - unsigned int *N_updates, double **linUpdates, - unsigned int **Updates_index) { - MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); - } -} diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp index 5314387..3585b7e 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -50,7 +50,7 @@ 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) { - std::cerr << "Called SM2 with updates " << N_updates << std::endl; + std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl; double C[Dim]; double D[Dim]; @@ -110,7 +110,7 @@ void SM2(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) { - std::cerr << "Called SM3 with updates " << N_updates << std::endl; + std::cerr << "Called SM3 with " << N_updates << " updates" << std::endl; double C[Dim]; double D[Dim]; @@ -176,7 +176,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) { - std::cerr << "Called SM4 with updates " << N_updates << std::endl; + std::cerr << "Called SM4 with " << N_updates << " updates" << std::endl; double C[Dim]; double D[Dim]; @@ -254,6 +254,7 @@ extern "C" { unsigned int **Updates_index) { SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); } + void SM4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, double **linUpdates, unsigned int **Updates_index) { diff --git a/src/SM_mod.f90 b/src/finterface_mod.f90 similarity index 100% rename from src/SM_mod.f90 rename to src/finterface_mod.f90 diff --git a/src/Helpers_mod.f90 b/src/helpers_mod.f90 similarity index 100% rename from src/Helpers_mod.f90 rename to src/helpers_mod.f90 diff --git a/tests/cMaponiA3_test_3x3_3.cpp b/tests/cMaponiA3_test_3x3_3.cpp index 20f7d45..fc278bc 100644 --- a/tests/cMaponiA3_test_3x3_3.cpp +++ b/tests/cMaponiA3_test_3x3_3.cpp @@ -1,5 +1,5 @@ // main.cpp -#include "SM_MaponiA3.hpp" +#include "SM_Maponi.hpp" #include "SM_Helpers.hpp" int main() { diff --git a/tests/run_test.sh b/tests/run_test.sh index c23676d..a4b81a9 100755 --- a/tests/run_test.sh +++ b/tests/run_test.sh @@ -12,7 +12,8 @@ ALGO=$1 START=$2 STOP=$3 BLOCKSIZE=329 -TOLERANCE=1e-3 +#TOLERANCE=1e-3 +TOLERANCE=$4 BLOCKS=$((STOP / BLOCKSIZE)) LEFT=$((STOP % BLOCKSIZE)) diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp index 005eca8..cd4a70a 100644 --- a/tests/test_h5.cpp +++ b/tests/test_h5.cpp @@ -1,39 +1,15 @@ #include "hdf5/serial/hdf5.h" #include "hdf5/serial/H5Cpp.h" -#include "SM_MaponiA3.hpp" -#include "SM_MaponiA3S.hpp" +#include "SM_Maponi.hpp" #include "SM_Standard.hpp" #include "SM_Helpers.hpp" using namespace H5; -#define DEBUG +// #define DEBUG const H5std_string FILE_NAME( "dataset.hdf5" ); -double residual_max(double * A, unsigned int Dim) { - double max = 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)); - delta = abs(delta); - if (delta > max) max = delta; - } - } - return max; -} - -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); @@ -92,8 +68,6 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { showMatrix(u, dim, "Updates"); #endif - - if (version == "maponia3") { MaponiA3(slater_inverse, dim, nupdates, u, col_update_index); } else if (version == "maponia3s") { @@ -124,7 +98,7 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { bool ok = is_identity(res, dim, tolerance); double res_max = residual_max(res, dim); - double res2 = residual2(res, dim); + double res2 = residual_frobenius2(res, dim); std::cout << "Residual = " << version << " " << cycle << " " << res_max << " " << res2 << std::endl; #ifdef DEBUG diff --git a/tests/vfc_test_h5.cpp b/tests/vfc_test_h5.cpp index 7b63b5d..23cae1d 100644 --- a/tests/vfc_test_h5.cpp +++ b/tests/vfc_test_h5.cpp @@ -11,7 +11,7 @@ #include -#include "SM_MaponiA3.hpp" +#include "SM_Maponi.hpp" #include "SM_Standard.hpp" #include "SM_Helpers.hpp" #include "vfc_probe.h" @@ -21,29 +21,6 @@ using namespace H5; const H5std_string FILE_NAME( "datasets/ci_dataset.hdf5" ); -double residual_max(double * A, unsigned int Dim) { - double max = 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)); - delta = abs(delta); - if (delta > max) max = delta; - } - } - return max; -} - -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); @@ -56,7 +33,6 @@ void read_double(H5File file, std::string key, double * data) { ds.close(); } - /* Return a vector containing all cycles to execute by reading a data file */ std::vector get_cycles_list(std::string path) { std::ifstream file_stream(path); @@ -157,7 +133,7 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) bool ok = is_identity(res, dim, 1e-3); double res_max = residual_max(res, dim); - double res2 = residual2(res, dim); + double res2 = residual_frobenius2(res, dim); #ifdef DEBUG showMatrix(res, dim, "Result");