diff --git a/Helpers.hpp b/Helpers.hpp index 7bc5efd..6dab844 100644 --- a/Helpers.hpp +++ b/Helpers.hpp @@ -7,9 +7,9 @@ using namespace std; template unsigned int getMaxIndex(T *vector, unsigned int size) { - unsigned int i; - unsigned int max = vector[0]; - unsigned int maxi = 0; + unsigned int i = 0; + unsigned int maxi = i; + unsigned int max = vector[maxi]; for (i = 1; i < size; i++) { if (vector[i] > max) { max = vector[i]; @@ -64,7 +64,7 @@ T *transpose(T *A, unsigned int M) { template T *matMul(T *A, T *B, unsigned int M) { - T *C = new T[M*M]; + T *C = new T[M*M] {0}; for (unsigned int i = 0; i < M; i++) { for (unsigned int j = 0; j < M; j++) { for (unsigned int k = 0; k < M; k++) { @@ -75,20 +75,6 @@ T *matMul(T *A, T *B, unsigned int M) { return C; } -template -T *matMul2(T *A, T *B, unsigned int M) { - T *C = new T[M*M]; - for (unsigned int i = 0; i < M; i++) { - for (unsigned int j = 0; j < M; j++) { - for (unsigned int k = 0; k < M; k++) { - C[i*M+j] += A[i*M+k] * B[k*M+j]; - } - } - } - showMatrix(C,M,"C"); - return C; -} - template T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) { T1 *C = new T1[M*M]; diff --git a/QMCChem_dataset_test.f90 b/QMCChem_dataset_test.f90 index 508c0f3..e4c4986 100644 --- a/QMCChem_dataset_test.f90 +++ b/QMCChem_dataset_test.f90 @@ -9,7 +9,7 @@ program QMCChem_dataset_test real(c_double), dimension(:,:), allocatable :: Updates real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t - call Read_dataset("datasets/update_cycle_8169.dat", & + call Read_dataset("datasets/update_cycle_13.dat", & cycle_id, & dim, & n_updates, & diff --git a/SM_MaponiA3.cpp b/SM_MaponiA3.cpp index 187ed65..6a1e6e7 100644 --- a/SM_MaponiA3.cpp +++ b/SM_MaponiA3.cpp @@ -4,15 +4,17 @@ #include "SM_MaponiA3.hpp" #include "Helpers.hpp" -void MaponiA3(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) { unsigned int k, l, lbar, i, j, tmp, component; - unsigned int *p = new unsigned int[N_updates + 1]; + unsigned int *p = new unsigned int[N_updates + 1] {0}; double alpha, beta; - double *breakdown = new double[N_updates + 1]; + double *breakdown = new double[N_updates + 1] {0}; double *Al = new double[Dim * Dim]; - p[0] = 0; + + // Populate update-order vector for (i = 0; i < N_updates; i++) { p[i + 1] = i + 1; } @@ -43,8 +45,8 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]); } lbar = getMaxIndex(breakdown, N_updates + 1); - // Reset breakdown back to 0 for next round to avoid case where its - // first element is always the largest + // Reset breakdown back to 0 for next round to avoid case where + // its first element is always the largest for (i = 0; i < N_updates + 1; i++) { breakdown[i] = 0; } @@ -60,14 +62,16 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, for (k = l + 1; k < N_updates + 1; k++) { alpha = ylk[l - 1][p[k]][component] / beta; for (i = 1; i < Dim + 1; i++) { - ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - alpha * ylk[l - 1][p[l]][i]; + ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] + - alpha * ylk[l - 1][p[l]][i]; } } } // Keep the memory location of the passed array 'Slater_inv' before - // 'Slater_inv' gets reassigned by 'matMul(...)' in the next line, by creating - // a new pointer 'copy' that points to whereever 'Slater_inv' points to now. + // 'Slater_inv' gets reassigned by 'matMul(...)' in the next line, + // by creating a new pointer 'copy' that points to whereever + // 'Slater_inv' points to now. double *copy = Slater_inv; // Construct A-inverse from A0-inverse and the ylk @@ -77,13 +81,15 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, beta = 1 + ylk[l][p[k]][component]; for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { - Al[i * Dim + j] = (i == j) - (j == component - 1) * ylk[l][p[k]][i + 1] / beta; + Al[i*Dim + j] = (i == j) - (j == component-1) + * ylk[l][p[k]][i + 1] / beta; } } Slater_inv = matMul(Al, Slater_inv, Dim); } - // Assign the new values of 'Slater_inv' to the old values in 'copy[][]' + // Assign the new values of 'Slater_inv' to the old values + // in 'copy[][]' for (i = 0; i < Dim; i++) { for (j = 0; j < Dim; j++) { copy[i * Dim + j] = Slater_inv[i * Dim + j]; diff --git a/tests/test.cpp b/tests/test.cpp index 199bdef..01ae2b5 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -48,11 +48,11 @@ int test_cycle(H5File file, int cycle) { /* Test */ #ifdef DEBUG - showMatrix(slater_matrix, dim, "Slater"); + showMatrix(slater_matrix, dim, "OLD Slater"); #endif #ifdef DEBUG - showMatrix(slater_inverse, dim, "Inverse"); + showMatrix(slater_inverse, dim, "OLD Inverse"); #endif for (j = 0; j < nupdates; j++) { @@ -65,21 +65,22 @@ int test_cycle(H5File file, int cycle) { MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index); #ifdef DEBUG - showMatrix(slater_matrix, dim, "Slater"); + showMatrix(slater_matrix, dim, "NEW Slater"); #endif #ifdef DEBUG - showMatrix(slater_inverse, dim, "Inverse"); + showMatrix(slater_inverse, dim, "NEW Inverse"); #endif - double * res = matMul2(slater_matrix, slater_inverse, dim); - bool ok = is_identity(res, dim, 1.0e-8); + double * res = matMul(slater_matrix, slater_inverse, dim); + bool ok = is_identity(res, dim, 0.5e-4); #ifdef DEBUG showMatrix(res, dim, "Result"); #endif - delete [] res, updates, col_update_index, slater_matrix, slater_inverse; + delete [] res, updates, col_update_index, + slater_matrix, slater_inverse; return ok; } @@ -95,10 +96,12 @@ int main(int argc, char **argv) { bool ok = test_cycle(file, cycle); if (ok) { - std::cerr << "ok -- 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; + std::cerr << "failed -- cycle " << std::to_string(cycle) + << std::endl; } return ok; }