diff --git a/include/SM_Helpers.hpp b/include/SM_Helpers.hpp index 049df9c..e8c308e 100644 --- a/include/SM_Helpers.hpp +++ b/include/SM_Helpers.hpp @@ -1,127 +1,122 @@ // SM_Helpers.hpp // Some usefull helper functions to support the Maponi algorithm. +#include +#include #include #include -#include -#include void Switch(unsigned int *p, unsigned int l, unsigned int lbar); -void selectLargestDenominator(unsigned int l, unsigned int N_updates, - unsigned int *Updates_index, unsigned int *p, - double ***ylk); +void selectLargestDenominator(unsigned int l, unsigned int N_updates, + unsigned int *Updates_index, unsigned int *p, + double ***ylk); -template -void showScalar(T scalar, std::string name) { - std::cout << name << " = " << scalar << std::endl << std::endl; +template void showScalar(T scalar, std::string name) { + std::cout << name << " = " << scalar << std::endl << std::endl; } -template +template void showVector(T *vector, unsigned int size, std::string name) { - std::cout << name << " = " << std::endl; - for (unsigned int i = 0; i < size; i++) { - std::cout << "[ " << vector[i] << " ]" << std::endl; - } - std::cout << std::endl; + std::cout << name << " = " << std::endl; + for (unsigned int i = 0; i < size; i++) { + std::cout << "[ " << vector[i] << " ]" << std::endl; + } + std::cout << std::endl; } -template +template void showMatrix(T *matrix, unsigned int M, std::string name) { - std::cout.precision(17); - std::cout << name << " = [" << std::endl; - for (unsigned int i = 0; i < M; i++) { - std::cout << "["; - for (unsigned int j = 0; j < M; j++) { - if (matrix[i*M + j] >= 0) { - std::cout << " " << matrix[i*M + j] << ","; - } - else { - std::cout << " " << matrix[i*M + j] << "," ; - } - } - std::cout << " ]," << std::endl; + std::cout.precision(17); + std::cout << name << " = [" << std::endl; + for (unsigned int i = 0; i < M; i++) { + std::cout << "["; + for (unsigned int j = 0; j < M; j++) { + if (matrix[i * M + j] >= 0) { + std::cout << " " << matrix[i * M + j] << ","; + } else { + std::cout << " " << matrix[i * M + j] << ","; + } } - std::cout << "]" << std::endl; - std::cout << std::endl; + std::cout << " ]," << std::endl; + } + std::cout << "]" << std::endl; + std::cout << std::endl; } -template -T *transpose(T *A, unsigned int M) { - T *B = new T[M*M]; - for (unsigned int i = 0; i < M; i++) { - for (unsigned int j = 0; j < M; j++) { - B[i*M + j] = A[i + j*M]; - } +template T *transpose(T *A, unsigned int M) { + T *B = new T[M * M]; + for (unsigned int i = 0; i < M; i++) { + for (unsigned int j = 0; j < M; j++) { + B[i * M + j] = A[i + j * M]; } - return B; + } + return B; } -template -void matMul(T *A, T *B, T *C, unsigned int M) { - memset(C, 0, M*M*sizeof(T)); - 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]; - } - } +template void matMul(T *A, T *B, T *C, unsigned int M) { + memset(C, 0, M * M * sizeof(T)); + 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]; + } } + } } -template +template T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) { - T1 *C = new T1[M*M]; - for (unsigned int i = 0; i < M; i++) { - for (unsigned int j = 0; j < M; j++) { - C[i*M+j] = vec1[i+1] * vec2[j]; - } + T1 *C = new T1[M * M]; + for (unsigned int i = 0; i < M; i++) { + for (unsigned int j = 0; j < M; j++) { + C[i * M + j] = vec1[i + 1] * vec2[j]; } - return C; + } + return C; } -template -T matDet(T **A, unsigned int M) { - int det = 0, p, h, k, i, j; - T **temp = new T*[M]; - for (int i = 0; i < M; i++) temp[i] = new T[M]; - if(M == 1) { - return A[0][0]; - } - else if(M == 2) { - det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]); - return det; - } - else { - for(p = 0; p < M; p++) { - h = 0; +template T matDet(T **A, unsigned int M) { + int det = 0, p, h, k, i, j; + T **temp = new T *[M]; + for (int i = 0; i < M; i++) + temp[i] = new T[M]; + if (M == 1) { + return A[0][0]; + } else if (M == 2) { + det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]); + return det; + } else { + for (p = 0; p < M; p++) { + h = 0; + k = 0; + for (i = 1; i < M; i++) { + for (j = 0; j < M; j++) { + if (j == p) { + continue; + } + temp[h][k] = A[i][j]; + k++; + if (k == M - 1) { + h++; k = 0; - for(i = 1; i < M; i++) { - for( j = 0; j < M; j++) { - if(j == p) { - continue; - } - temp[h][k] = A[i][j]; - k++; - if(k == M-1) { - h++; - k = 0; - } - } - } - det = det + A[0][p] * pow(-1, p) * matDet(temp, M-1); + } } - return det; + } + det = det + A[0][p] * pow(-1, p) * matDet(temp, M - 1); } - delete [] temp; + return det; + } + delete[] temp; } -template -bool is_identity(T *A, unsigned int M, double tolerance) { - for (unsigned int i = 0; i < M; i++) { - for (unsigned int j = 0; j < M; j++) { - if (i==j && fabs(A[i*M+j]-1) > tolerance) return false; - if (i!=j && fabs(A[i*M+j]) > tolerance) return false; - } - } - return true; +template bool is_identity(T *A, unsigned int M, double tolerance) { + for (unsigned int i = 0; i < M; i++) { + for (unsigned int j = 0; j < M; j++) { + if (i == j && fabs(A[i * M + j] - 1) > tolerance) + return false; + if (i != j && fabs(A[i * M + j]) > tolerance) + return false; + } + } + return true; } diff --git a/include/SM_MaponiA3.hpp b/include/SM_MaponiA3.hpp index 8d8867a..c4c566b 100644 --- a/include/SM_MaponiA3.hpp +++ b/include/SM_MaponiA3.hpp @@ -1,3 +1,2 @@ -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); diff --git a/include/SM_MaponiA3S.hpp b/include/SM_MaponiA3S.hpp index 86775fc..9085b7e 100644 --- a/include/SM_MaponiA3S.hpp +++ b/include/SM_MaponiA3S.hpp @@ -1,3 +1,2 @@ -void MaponiA3S(double *Slater_inv, unsigned int Dim, - unsigned int N_updates, double *Updates, - unsigned int *Updates_index); +void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index); diff --git a/include/SM_Standard.hpp b/include/SM_Standard.hpp index 0332ec7..5541517 100644 --- a/include/SM_Standard.hpp +++ b/include/SM_Standard.hpp @@ -1,12 +1,12 @@ // Naïve Sherman Morrison void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index); + 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); + 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); + double *Updates, unsigned int *Updates_index); diff --git a/src/SM_Helpers.cpp b/src/SM_Helpers.cpp index 185c367..4e07f3d 100644 --- a/src/SM_Helpers.cpp +++ b/src/SM_Helpers.cpp @@ -1,15 +1,15 @@ #include "SM_Helpers.hpp" void Switch(unsigned int *p, unsigned int l, unsigned int lbar) { - unsigned int tmp = p[l+1]; - p[l+1] = p[lbar]; + unsigned int tmp = p[l + 1]; + p[l + 1] = p[lbar]; p[lbar] = tmp; } -void selectLargestDenominator(unsigned int l, unsigned int N_updates, - unsigned int *Updates_index, unsigned int *p, - double ***ylk) { - unsigned int lbar = l+1, max =0; +void selectLargestDenominator(unsigned int l, unsigned int N_updates, + unsigned int *Updates_index, unsigned int *p, + double ***ylk) { + unsigned int lbar = l + 1, max = 0; unsigned int index = 0, component = 0; unsigned int tmp = 0; double breakdown = 0; @@ -17,11 +17,12 @@ void selectLargestDenominator(unsigned int l, unsigned int N_updates, index = p[j]; component = Updates_index[index - 1]; breakdown = abs(1 + ylk[l][index][component]); - #ifdef DEBUG +#ifdef DEBUG std::cout << "Inside selectLargestDenominator()" << std::endl; - std::cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << std::endl; + std::cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" + << component << "]) = " << breakdown << std::endl; std::cout << std::endl; - #endif +#endif if (breakdown > max) { max = breakdown; lbar = j; diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index 1196ed5..3495bec 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -6,19 +6,18 @@ // #define DEBUG -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) { /* DECLARE AND INITIALISE ARRAYS */ unsigned int k, l, i, j, component; - unsigned int *p = new unsigned int[N_updates + 1] {0}; + 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 *next = new double[Dim * Dim]{0}; double *last = Slater_inv, *tmp; // Populate update-order vector @@ -31,101 +30,99 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, 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}; + 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 +#ifdef DEBUG std::cout << "Compute y0k: " << std::endl; std::cout << "ylk[0][" << k << "][:]"; std::cout << std::endl; - #endif +#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)]; + ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] * + Updates[(k - 1) * Dim + (j - 1)]; } } - #ifdef DEBUG +#ifdef DEBUG showVector(ylk[0][k], Dim, ""); - #endif +#endif } // Calculate the {y_{l,k}} from the {y_{0,k}} for (l = 0; l < N_updates; l++) { - #ifdef DEBUG +#ifdef DEBUG std::cout << "In outer compute-ylk-loop: l = " << l << std::endl; std::cout << std::endl; - #endif +#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; + 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 << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component + << "] = " << beta << std::endl; std::cout << std::endl; - #endif +#endif if (fabs(beta) < 1e-3) { std::cerr << "Break-down occured." << std::endl; } - // Compute intermediate update to Slater_inv - #ifdef DEBUG +// 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 << "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 +#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] / beta; + Al[i * Dim + j] = + (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta; } } matMul(Al, last, next, Dim); tmp = next; next = last; last = tmp; - #ifdef DEBUG +#ifdef DEBUG showMatrix(last, Dim, "last"); - #endif +#endif // For given l != 0 compute the next {y_{l,k}} - for (k = l+2; k < N_updates+1; k++) { + for (k = l + 2; k < N_updates + 1; k++) { alpha = ylk[l][p[k]][component] / beta; - #ifdef DEBUG +#ifdef DEBUG std::cout << "Inside k-loop: k = " << k << std::endl; - std::cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << std::endl; + std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl; std::cout << std::endl; - #endif +#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]; + 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)); - - + 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]; @@ -136,11 +133,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, } 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); +} } diff --git a/src/SM_MaponiA3S.cpp b/src/SM_MaponiA3S.cpp index 249f018..24a50c2 100644 --- a/src/SM_MaponiA3S.cpp +++ b/src/SM_MaponiA3S.cpp @@ -6,19 +6,17 @@ #define DEBUG -void MaponiA3S(double *Slater_inv, unsigned int Dim, - unsigned int N_updates, double *Updates, - unsigned int *Updates_index) { - +void MaponiA3S(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}; + 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 *next = new double[Dim * Dim]{0}; double *last = Slater_inv, *tmp; // Populate update-order vector @@ -31,101 +29,99 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, 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}; + 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 +#ifdef DEBUG std::cout << "Compute y0k: " << std::endl; std::cout << "ylk[0][" << k << "][:]"; std::cout << std::endl; - #endif +#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)]; + ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] * + Updates[(k - 1) * Dim + (j - 1)]; } } - #ifdef DEBUG +#ifdef DEBUG showVector(ylk[0][k], Dim, ""); - #endif +#endif } // Calculate the {y_{l,k}} from the {y_{0,k}} for (l = 0; l < N_updates; l++) { - #ifdef DEBUG +#ifdef DEBUG std::cout << "In outer compute-ylk-loop: l = " << l << std::endl; std::cout << std::endl; - #endif +#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; + 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 << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component + << "] = " << beta << std::endl; std::cout << std::endl; - #endif +#endif if (fabs(beta) < 1e-3) { std::cerr << "Break-down occured." << std::endl; } - // Compute intermediate update to Slater_inv - #ifdef DEBUG +// 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 << "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 +#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] / beta; + Al[i * Dim + j] = + (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta; } } matMul(Al, last, next, Dim); tmp = next; next = last; last = tmp; - #ifdef DEBUG +#ifdef DEBUG showMatrix(last, Dim, "last"); - #endif +#endif // For given l != 0 compute the next {y_{l,k}} - for (k = l+2; k < N_updates+1; k++) { + for (k = l + 2; k < N_updates + 1; k++) { alpha = ylk[l][p[k]][component] / beta; - #ifdef DEBUG +#ifdef DEBUG std::cout << "Inside k-loop: k = " << k << std::endl; - std::cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << std::endl; + std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl; std::cout << std::endl; - #endif +#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]; + 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)); - - + 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]; @@ -136,11 +132,9 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, } 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 2354450..b3294ec 100644 --- a/src/SM_Standard.cpp +++ b/src/SM_Standard.cpp @@ -11,46 +11,38 @@ double threshold() { // Naïve Sherman Morrison void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index) -{ - + 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 - while (l < N_updates) - { + while (l < N_updates) { // C = A^{-1} x U_l - for (unsigned int i = 0; i < Dim; i++) - { + for (unsigned int i = 0; i < Dim; i++) { C[i] = 0; - for (unsigned int j = 0; j < Dim; j++) - { + 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) < threshold()) - { - std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; + if (fabs(den) < threshold()) { + std::cerr << "Breakdown condition triggered at " << Updates_index[l] + << std::endl; } double iden = 1 / den; // D = v^T x A^{-1} - for (unsigned int j = 0; j < Dim; j++) - { + 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++) - { + 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; } @@ -63,9 +55,7 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, // 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) -{ - + double *Updates, unsigned int *Updates_index) { std::cerr << "Called SM2 with updates " << N_updates << std::endl; double C[Dim]; double D[Dim]; @@ -76,27 +66,23 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, unsigned int l = 0; // For each update - while (l < N_updates) - { + while (l < N_updates) { // C = A^{-1} x U_l - for (unsigned int i = 0; i < Dim; i++) - { + for (unsigned int i = 0; i < Dim; i++) { C[i] = 0; - for (unsigned int j = 0; j < Dim; j++) - { + 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) < threshold()) - { - std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; + if (fabs(den) < threshold()) { + std::cerr << "Breakdown condition triggered at " << Updates_index[l] + << std::endl; // U_l = U_l / 2 (do the split) - for (unsigned int j = 0; j < Dim; j++) - { + for (unsigned int j = 0; j < Dim; j++) { later_updates[later * Dim + j] = Updates[l * Dim + j] / 2.0; C[j] /= 2.0; } @@ -108,16 +94,13 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double iden = 1 / den; // D = v^T x A^{-1} - for (unsigned int j = 0; j < Dim; j++) - { + 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++) - { + 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; } @@ -125,17 +108,14 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, l += 1; } - if (later > 0) - { + 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) -{ - + double *Updates, unsigned int *Updates_index) { std::cerr << "Called SM3 with updates " << N_updates << std::endl; double C[Dim]; double D[Dim]; @@ -146,26 +126,22 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, unsigned int l = 0; // For each update - while (l < N_updates) - { + while (l < N_updates) { // C = A^{-1} x U_l - for (unsigned int i = 0; i < Dim; i++) - { + for (unsigned int i = 0; i < Dim; i++) { C[i] = 0; - for (unsigned int j = 0; j < Dim; j++) - { + 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) < threshold()) - { - std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; + if (fabs(den) < threshold()) { + std::cerr << "Breakdown condition triggered at " << Updates_index[l] + << std::endl; - for (unsigned int j = 0; j < Dim; j++) - { + for (unsigned int j = 0; j < Dim; j++) { later_updates[later * Dim + j] = Updates[l * Dim + j]; } later_index[later] = Updates_index[l]; @@ -176,16 +152,13 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double iden = 1 / den; // D = v^T x A^{-1} - for (unsigned int j = 0; j < Dim; j++) - { + 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++) - { + 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; } @@ -193,8 +166,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, l += 1; } - if (later > 0) - { + if (later > 0) { SM3(Slater_inv, Dim, later, later_updates, later_index); } }