diff --git a/include/Helpers.hpp b/include/Helpers.hpp index 015c2ee..2384db1 100644 --- a/include/Helpers.hpp +++ b/include/Helpers.hpp @@ -56,6 +56,45 @@ void showMatrix(T *matrix, unsigned int M, std::string name) { std::cout << std::endl; } +template +void showMatrix2(T *matrix, unsigned int M, unsigned int N, 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 < N; j++) { + if (matrix[i * N + j] >= 0) { + std::cout << " " << matrix[i * N + j] << ","; + } else { + std::cout << " " << matrix[i * N + j] << ","; + } + } + std::cout << " ]," << std::endl; + } + std::cout << "]" << std::endl; + std::cout << std::endl; +} + + +template +void showMatrixNS(T *matrix, unsigned int M, unsigned int N, 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 < N; j++) { + if (matrix[i * N + j] >= 0) { + std::cout << " " << matrix[i * N + j] << ","; + } else { + std::cout << " " << matrix[i * N + j] << ","; + } + } + 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++) { @@ -77,6 +116,18 @@ template void matMul(T *A, T *B, T *C, unsigned int M) { } } +template +static inline void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) { + for(unsigned int i = 0; i < M; i++) { + for(unsigned int j = 0; j < P; j++) { + C[i * P + j] = 0; + for(unsigned int k = 0; k < N; k++) { + C[i * P + j] += A[i * N + k] * B[k * P + j]; + } + } + } +} + template T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) { T1 *C = new T1[M * M]; diff --git a/include/Woodbury.hpp b/include/Woodbury.hpp index 001d912..c03e509 100644 --- a/include/Woodbury.hpp +++ b/include/Woodbury.hpp @@ -1,7 +1,7 @@ // Woodbury 2x2 kernel -void WB2(double *Slater_inv, unsigned int Dim, +bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, unsigned int *Updates_index); // Woodbury 3x3 kernel -void WB3(double *Slater_inv, unsigned int Dim, +bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, unsigned int *Updates_index); diff --git a/smvars.sh b/smvars.sh index 6e846ee..2e9e868 100644 --- a/smvars.sh +++ b/smvars.sh @@ -3,6 +3,10 @@ unset THRESHOLD unset MKL ENV=$1 +echo +echo "Sherman-Morrison-Woodbury parameters" +echo "------------------------------------" + ## Set Sherman-Morrison root dir PWD=$(pwd) SRCDIR=$(dirname $BASH_SOURCE) diff --git a/src/SMWB.cpp b/src/SMWB.cpp index 0590ff5..670920b 100644 --- a/src/SMWB.cpp +++ b/src/SMWB.cpp @@ -5,12 +5,15 @@ // Sherman-Morrison-Woodbury kernel 1 // WB2, WB3, SM2 mixing scheme 1 -void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index) { +void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) { std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl; - WB3(Slater_inv, Dim, Updates, Updates_index); - WB2(Slater_inv, Dim, Updates, Updates_index); - SM2(Slater_inv, Dim, N_updates, Updates, Updates_index); + + bool ok; + ok = WB2(Slater_inv, Dim, Updates, Updates_index); + if (!ok) { + std::cerr << "Woodbury kernel failed!" << std::endl; + SM2(Slater_inv, Dim, N_updates, Updates, Updates_index); + } } extern "C" { diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp index 21b9d9d..980464b 100644 --- a/src/Woodbury.cpp +++ b/src/Woodbury.cpp @@ -5,65 +5,126 @@ // (S + U * V)^{-1} = S^{-1} - S^{-1} * U * B^{-1} * V * S^{-1} // B := 1 + V * C, 2 x 2 // C := S^{-1} * U, dim x 2 +// D := V * S^{-1}, 2 x dim +// +// All matrices are stored in row-major order + #include "Woodbury.hpp" #include "Helpers.hpp" -// Woodbury 2x2 kernel: -void WB2(double *Slater_inv, unsigned int Dim, - double *Updates, unsigned int *Updates_index) { +// Woodbury 2x2 kernel: +bool WB2(double *Slater_inv, unsigned int Dim, double *Updates, + unsigned int *Updates_index) { std::cerr << "Called Woodbury 2x2 kernel" << std::endl; // Construct V from Updates_index - unsigned int *V = new unsigned int[2 * Dim]{0}; - for (unsigned int i = 0; i < Dim; i++) { - if (i == Updates_index[0] - 1) V[i] = 1; - if (Dim + i == Updates_index[1] - 1) V[Dim + i] = 1; - } + unsigned int V[2 * Dim]{0}; // 2 x Dim matrix stored in row-major order + V[Updates_index[0] - 1] = 1; + V[Dim + Updates_index[1] - 1] = 1; - // Compute B from U, V and Slater_inv + // Compute C + double C[2 * Dim]{0}; + matMul2(Slater_inv, Updates, C, Dim, Dim, 2); + // Compute B double B[4]{0}; - double Binv[4]{0}; - double C[4]{0}; - double D[2*Dim]{0}; - // Compute C - for (unsigned int i = 0; i < Dim; i++) { - for (unsigned int j = 0; j < 2; j++) { - for (unsigned int k = 0; k < Dim; k++) { - C[i * Dim + j] += Slater_inv[i * Dim + k] * Updates[k * Dim + j]; - } - } - } - // Compute B - for (unsigned int i = 0; i < 2; i++) { - for (unsigned int j = 0; j < 2; j++) { - for (unsigned int k = 0; k < Dim; k++) { - B[i * Dim + j] += (i == j) + V[i * Dim + k] * C[k * Dim + j]; - } - } - } + matMul2(V, C, B, 2, Dim, 2); + // Compute 1 + B + B[0] += 1; + B[3] += 1; - // Invert B with explicit formula for 2x2 inversion + // Invert 1 + B with explicit formula for 2x2 inversion double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]); + double Binv[4]{0}; Binv[0] = idet * B[3]; Binv[1] = -1.0 * idet * B[1]; Binv[2] = -1.0 * idet * B[2]; Binv[3] = idet * B[0]; + // Check if determinant of inverted matrix is not zero + double det = B[0] * B[3] - B[1] * B[2]; + if (std::fabs(det) < threshold()) { + std::cerr << "Determinant approached 0!" << std::endl; + return false; + } + // Compute (S + U * V)^{-1} with Woobury identity - + double D[2 * Dim]{0}; + matMul2(V, Slater_inv, D, 2, Dim, Dim); + double tmp[2 * Dim]{0}; + matMul2(Binv, D, tmp, 2, 2, Dim); + double tmp2[Dim * Dim]{0}; + matMul2(C, tmp, tmp2, Dim, 2, Dim); + for (unsigned int i = 0; i < Dim * Dim; i++) { + Slater_inv[i] -= tmp2[i]; + } + return true; } // Woodbury 3x3 kernel -void WB3(double *Slater_inv, unsigned int Dim, - double *Updates, unsigned int *Updates_index) { - std::cerr << "Called Woodbury 3x3 kernel" << std::endl; +bool WB3(double *Slater_inv, unsigned int Dim, double *Updates, + unsigned int *Updates_index) { + std::cerr << "Called Woodbury 3x3 kernel" << std::endl; + + // Construct V from Updates_index + unsigned int V[3 * Dim]{0}; // 2 x Dim matrix stored in row-major order + V[Updates_index[0] - 1] = 1; + V[Dim + Updates_index[1] - 1] = 1; + V[2 * Dim + Updates_index[2] - 1] = 1; - // Construct V from Updates_index + // Compute C + double C[3 * Dim]{0}; + matMul2(Slater_inv, Updates, C, Dim, Dim, 3); + // Compute B + double B[9]{0}; + matMul2(V, C, B, 3, Dim, 3); + // Compute 1 + B + B[0] += 1; + B[4] += 1; + B[8] += 1; - // Compute B from U, V and Slater_inv + double Binv[9]; + Binv[0] = B[4] * B[8] - B[5] * B[7]; + Binv[3] = B[5] * B[6] - B[3] * B[8]; + Binv[6] = B[3] * B[7] - B[4] * B[6]; - // Invert B with explicit formula for 3x3 inversion + Binv[1] = B[2] * B[7] - B[1] * B[8]; + Binv[4] = B[0] * B[8] - B[2] * B[6]; + Binv[7] = B[1] * B[6] - B[0] * B[7]; - // Compute (S + U * V)^{-1} with Woobury identity - - } + Binv[2] = B[1] * B[5] - B[2] * B[4]; + Binv[5] = B[2] * B[3] - B[0] * B[5]; + Binv[8] = B[0] * B[4] - B[1] * B[3]; + + // Check if determinant of inverted matrix is not zero + // If so, exigt and return false. + double det; + det = B[0] * (B[4] * B[8] - B[5] * B[7]) - + B[1] * (B[3] * B[8] - B[5] * B[6]) + B[2] * (B[3] * B[7] - B[4] * B[6]); + if (std::fabs(det) < threshold()) { + std::cerr << "Determinant approached 0!" << std::endl; + return false; + } + + // Compute (S + U * V)^{-1} with Woobury identity + double D[3 * Dim]{0}; + matMul2(V, Slater_inv, D, 3, Dim, Dim); + double tmp[3 * Dim]{0}; + matMul2(Binv, D, tmp, 3, 3, Dim); + double tmp2[Dim * Dim]{0}; + matMul2(C, tmp, tmp2, Dim, 3, Dim); + for (unsigned int i = 0; i < Dim * Dim; i++) { + Slater_inv[i] -= tmp2[i]; + } + return true; +} + +extern "C" { +bool WB2_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates, + unsigned int **Updates_index) { + WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index); +} +bool WB3_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates, + unsigned int **Updates_index) { + WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index); +} +} diff --git a/todo.txt b/todo.txt deleted file mode 100644 index 5a579c0..0000000 --- a/todo.txt +++ /dev/null @@ -1,10 +0,0 @@ --=={TODO}==- - -* Looking at the eigenvalues of the intermediate matrix after a problematic update is -* plotting the ratio of the moduli of the maximum and the minimum eigenvalues of - -* Use valgrind to find the problem with MaponiA3S -* Contact Claudia if she would be interested in sharing datasets from Champ to test with SM - -* Make a representative selection of update cycles for presentatoin puroses -