diff --git a/include/SM_Helpers.hpp b/include/SM_Helpers.hpp index 73748b9..f9a7458 100644 --- a/include/SM_Helpers.hpp +++ b/include/SM_Helpers.hpp @@ -128,7 +128,7 @@ template bool is_identity(T *A, unsigned int M, double tolerance) { return true; } -template T norm_max(T * A, unsigned int Dim) { +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++) { @@ -142,37 +142,37 @@ template T norm_max(T * A, unsigned int Dim) { return res; } -template T norm_frobenius2(T * A, unsigned int Dim) { +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; + 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) { +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); - res += delta*delta; + 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/src/SM_Maponi.cpp b/src/SM_Maponi.cpp index 7b8263e..7ea5607 100644 --- a/src/SM_Maponi.cpp +++ b/src/SM_Maponi.cpp @@ -80,7 +80,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; } - double ibeta = 1.0 / beta; + double ibeta = 1.0 / beta; // Compute intermediate update to Slater_inv #ifdef DEBUG @@ -184,7 +184,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, } } #ifdef DEBUG - showVector(ylk[0][k], Dim+1, ""); + showVector(ylk[0][k], Dim + 1, ""); #endif } @@ -237,7 +237,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta; } } - matMul(Al, last, next, Dim); tmp = next; next = last; last = tmp; + matMul(Al, last, next, Dim); + tmp = next; + next = last; + last = tmp; #ifdef DEBUG showMatrix(last, Dim, "last"); #endif @@ -275,17 +278,17 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, } 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); - } +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, - unsigned int **Updates_index) { - MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); - } +void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim, + unsigned int *N_updates, double **linUpdates, + unsigned int **Updates_index) { + MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} } diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp index 3585b7e..cb4cc2e 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -173,7 +173,8 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, } // Sherman Morrison, mix between SM3 + SM2 -// Leave zero denominators for later (SM3), and when none are left then split (SM2) +// 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 " << N_updates << " updates" << std::endl; @@ -237,27 +238,23 @@ void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, } extern "C" { - void SM1_f(double **linSlater_inv, unsigned int *Dim, - unsigned int *N_updates, double **linUpdates, - unsigned int **Updates_index) { - SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); - } - - void SM2_f(double **linSlater_inv, unsigned int *Dim, - unsigned int *N_updates, double **linUpdates, - unsigned int **Updates_index) { - SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); - } - - void SM3_f(double **linSlater_inv, unsigned int *Dim, - unsigned int *N_updates, double **linUpdates, - 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) { - SM4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); - } +void SM1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, + double **linUpdates, unsigned int **Updates_index) { + SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} + +void SM2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, + double **linUpdates, unsigned int **Updates_index) { + SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} + +void SM3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, + double **linUpdates, 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) { + SM4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index); +} } diff --git a/tests/cMaponiA3_test_3x3_3.cpp b/tests/cMaponiA3_test_3x3_3.cpp index fc278bc..bcd363a 100644 --- a/tests/cMaponiA3_test_3x3_3.cpp +++ b/tests/cMaponiA3_test_3x3_3.cpp @@ -1,53 +1,61 @@ // main.cpp -#include "SM_Maponi.hpp" #include "SM_Helpers.hpp" +#include "SM_Maponi.hpp" int main() { - unsigned int M = 3; // Dimension of the Slater-matrix - unsigned int i, j; // Indices for iterators + unsigned int M = 3; // Dimension of the Slater-matrix + unsigned int i, j; // Indices for iterators - // Declare, allocate all vectors and matrices and fill them with zeros - unsigned int *Ar_index = new unsigned int[M]; - double *A = new double[M*M]; // The matrix to be inverted - double *A0 = new double[M*M]; // A diagonal matrix with the digonal elements of A - double *Ar = new double[M*M]; // The update matrix - double *A0_inv = new double[M*M]; // The inverse + // Declare, allocate all vectors and matrices and fill them with zeros + unsigned int *Ar_index = new unsigned int[M]; + double *A = new double[M * M]; // The matrix to be inverted + double *A0 = + new double[M * M]; // A diagonal matrix with the digonal elements of A + double *Ar = new double[M * M]; // The update matrix + double *A0_inv = new double[M * M]; // The inverse - // Fill with zeros - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - A0[i*M+j] = 0; - Ar[i*M+j] = 0; - A0_inv[i*M+j] = 0; - } + // Fill with zeros + for (i = 0; i < M; i++) { + for (j = 0; j < M; j++) { + A0[i * M + j] = 0; + Ar[i * M + j] = 0; + A0_inv[i * M + j] = 0; } + } - // Initialize A with M=3 and fill acc. to Eq. (17) from paper - A[0] = 1; A[3] = 1; A[6] = -1; - A[1] = 1; A[4] = 1; A[7] = 0; - A[2] = -1; A[5] = 0; A[8] = -1; + // Initialize A with M=3 and fill acc. to Eq. (17) from paper + A[0] = 1; + A[3] = 1; + A[6] = -1; + A[1] = 1; + A[4] = 1; + A[7] = 0; + A[2] = -1; + A[5] = 0; + A[8] = -1; - showMatrix(A, M, "A"); + showMatrix(A, M, "A"); - // Initialize the diagonal matrix A0, - // the inverse of A0_inv of diagonal matrix A0_inv - // and the update matrix Ar - for (i = 0; i < M; i++) { - A0[i*M + i] = A[i*M + i]; - A0_inv[i*M + i] = 1.0/A[i*M + i]; - Ar_index[i] = i+1; // ! First column needs to start with 1 ! - for (j = 0; j < M; j++) { - Ar[i*M + j] = A[i*M + j] - A0[i*M + j]; - } + // Initialize the diagonal matrix A0, + // the inverse of A0_inv of diagonal matrix A0_inv + // and the update matrix Ar + for (i = 0; i < M; i++) { + A0[i * M + i] = A[i * M + i]; + A0_inv[i * M + i] = 1.0 / A[i * M + i]; + Ar_index[i] = i + 1; // ! First column needs to start with 1 ! + for (j = 0; j < M; j++) { + Ar[i * M + j] = A[i * M + j] - A0[i * M + j]; } + } - // Define pointers dim and n_updates to use in Sherman-Morrison(...) function call - MaponiA3(A0_inv, M, M, Ar, Ar_index); - showMatrix(A0_inv, M, "A0_inv"); + // Define pointers dim and n_updates to use in Sherman-Morrison(...) function + // call + MaponiA3(A0_inv, M, M, Ar, Ar_index); + showMatrix(A0_inv, M, "A0_inv"); - // Deallocate all vectors and matrices - delete [] A, A0, A0_inv, Ar, Ar_index; + // Deallocate all vectors and matrices + delete[] A, A0, A0_inv, Ar, Ar_index; - return 0; + return 0; } diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp index cd4a70a..6fe6bea 100644 --- a/tests/test_h5.cpp +++ b/tests/test_h5.cpp @@ -1,22 +1,22 @@ -#include "hdf5/serial/hdf5.h" #include "hdf5/serial/H5Cpp.h" +#include "hdf5/serial/hdf5.h" +#include "SM_Helpers.hpp" #include "SM_Maponi.hpp" #include "SM_Standard.hpp" -#include "SM_Helpers.hpp" using namespace H5; // #define DEBUG -const H5std_string FILE_NAME( "dataset.hdf5" ); +const H5std_string FILE_NAME("dataset.hdf5"); -void read_int(H5File file, std::string key, unsigned int * data) { +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) { +void read_double(H5File file, std::string key, double *data) { DataSet ds = file.openDataSet(key); ds.read(data, PredType::IEEE_F64LE); ds.close(); @@ -32,19 +32,19 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { read_int(file, group + "/slater_matrix_dim", &dim); read_int(file, group + "/nupdates", &nupdates); - double * slater_matrix = new double[dim*dim]; + double *slater_matrix = new double[dim * dim]; read_double(file, group + "/slater_matrix", slater_matrix); - double * slater_inverse = new double[dim*dim]; + double *slater_inverse = new double[dim * dim]; read_double(file, group + "/slater_inverse", slater_inverse); - unsigned int * col_update_index = new unsigned int[nupdates]; + 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]; + double *updates = new double[nupdates * dim]; read_double(file, group + "/updates", updates); - double * u = new double[nupdates*dim]; + double *u = new double[nupdates * dim]; /* Test */ #ifdef DEBUG @@ -55,8 +55,9 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { 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]; + u[i + j * dim] = + updates[i + j * dim] - slater_matrix[i * dim + (col - 1)]; + slater_matrix[i * dim + (col - 1)] = updates[i + j * dim]; } } @@ -93,20 +94,20 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { showMatrix(slater_inverse, dim, "NEW Inverse"); #endif - double * res = new double[dim*dim] {0}; + 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); - std::cout << "Residual = " << version << " " << cycle << " " << res_max << " " << res2 << std::endl; + std::cout << "Residual = " << version << " " << cycle << " " << res_max << " " + << res2 << std::endl; #ifdef DEBUG showMatrix(res, dim, "Result"); #endif - delete [] res, updates, u, col_update_index, - slater_matrix, slater_inverse; + delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse; return ok; } @@ -114,7 +115,9 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) { int main(int argc, char **argv) { if (argc != 5) { std::cerr << "Execute from within 'datasets/'" << std::endl; - std::cerr << "usage: test_h5 " << std::endl; + std::cerr + << "usage: test_h5 " + << std::endl; return 1; } std::string version(argv[1]); @@ -124,15 +127,12 @@ int main(int argc, char **argv) { H5File file(FILE_NAME, H5F_ACC_RDONLY); bool ok; - for (int cycle = start_cycle; cycle < stop_cycle+1; cycle++) { + for (int cycle = start_cycle; cycle < stop_cycle + 1; cycle++) { ok = test_cycle(file, cycle, version, tolerance); if (ok) { - std::cerr << "ok -- cycle " << std::to_string(cycle) - << std::endl; - } - else { - std::cerr << "failed -- cycle " << std::to_string(cycle) - << std::endl; + std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl; + } else { + std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl; } } diff --git a/tests/vfc_test_h5.cpp b/tests/vfc_test_h5.cpp index 23cae1d..91109cf 100644 --- a/tests/vfc_test_h5.cpp +++ b/tests/vfc_test_h5.cpp @@ -3,31 +3,30 @@ // cycles in a CSV file, instead of accepting a start and an end cycle (which // makes it easier to select the exact cycles we are interested in with vfc_ci). -#include #include +#include -#include #include #include +#include - +#include "SM_Helpers.hpp" #include "SM_Maponi.hpp" #include "SM_Standard.hpp" -#include "SM_Helpers.hpp" #include "vfc_probe.h" using namespace H5; // #define DEBUG -const H5std_string FILE_NAME( "datasets/ci_dataset.hdf5" ); +const H5std_string FILE_NAME("datasets/ci_dataset.hdf5"); -void read_int(H5File file, std::string key, unsigned int * data) { +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) { +void read_double(H5File file, std::string key, double *data) { DataSet ds = file.openDataSet(key); ds.read(data, PredType::IEEE_F64LE); ds.close(); @@ -42,14 +41,15 @@ std::vector get_cycles_list(std::string path) { std::string cycle_str; std::vector cycles_list = {}; - while(string_stream >> cycle_str) { + while (string_stream >> cycle_str) { cycles_list.push_back(std::stoi(cycle_str)); } return cycles_list; } -int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) { +int test_cycle(H5File file, int cycle, std::string version, + vfc_probes *probes) { /* Read the data */ @@ -59,11 +59,12 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) // being zero-padded. This is used when calling vfc_put_probe later on. std::string zero_padded_group = std::to_string(cycle); zero_padded_group = "cycle_" + - std::string(5 - zero_padded_group.length(), '0') + zero_padded_group; + std::string(5 - zero_padded_group.length(), '0') + + zero_padded_group; - try{ + try { file.openGroup(group); - } catch(H5::Exception& e){ + } catch (H5::Exception &e) { std::cerr << "group " << group << "not found" << std::endl; return 0; } @@ -72,21 +73,20 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) read_int(file, group + "/slater_matrix_dim", &dim); read_int(file, group + "/nupdates", &nupdates); - - double * slater_matrix = new double[dim*dim]; + double *slater_matrix = new double[dim * dim]; read_double(file, group + "/slater_matrix", slater_matrix); - double * slater_inverse = new double[dim*dim]; + double *slater_inverse = new double[dim * dim]; read_double(file, group + "/slater_inverse", slater_inverse); - //slater_inverse = transpose(slater_inverse, dim); + // slater_inverse = transpose(slater_inverse, dim); - unsigned int * col_update_index = new unsigned int[nupdates]; + 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]; + double *updates = new double[nupdates * dim]; read_double(file, group + "/updates", updates); - double * u = new double[nupdates*dim]; + double *u = new double[nupdates * dim]; /* Test */ #ifdef DEBUG @@ -100,8 +100,9 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) 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]; + u[i + j * dim] = + updates[i + j * dim] - slater_matrix[i * dim + (col - 1)]; + slater_matrix[i * dim + (col - 1)] = updates[i + j * dim]; } } @@ -128,7 +129,7 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) showMatrix(slater_inverse, dim, "NEW Inverse"); #endif - double * res = new double[dim*dim] {0}; + double *res = new double[dim * dim]{0}; matMul(slater_matrix, slater_inverse, res, dim); bool ok = is_identity(res, dim, 1e-3); @@ -139,11 +140,11 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) showMatrix(res, dim, "Result"); #endif - vfc_put_probe(probes, &(zero_padded_group)[0], &("res_max_" + version)[0], res_max); + vfc_put_probe(probes, &(zero_padded_group)[0], &("res_max_" + version)[0], + res_max); vfc_put_probe(probes, &(zero_padded_group)[0], &("res2_" + version)[0], res2); - delete [] res, updates, u, col_update_index, - slater_matrix, slater_inverse; + delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse; return ok; } @@ -165,12 +166,9 @@ int main(int argc, char **argv) { for (int i = 0; i < cycles_list.size(); i++) { ok = test_cycle(file, cycles_list[i], version, &probes); if (ok) { - std::cout << "ok -- cycle " << std::to_string(i) - << std::endl; - } - else { - std::cerr << "failed -- cycle " << std::to_string(i) - << std::endl; + std::cout << "ok -- cycle " << std::to_string(i) << std::endl; + } else { + std::cerr << "failed -- cycle " << std::to_string(i) << std::endl; } }