From aa6895a4aadb5c4ac10694fd49895dfabf1fff80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 3 Feb 2021 12:13:09 +0100 Subject: [PATCH] Put helper functions in separate header file from the Sherman-Morrison algo. --- Helpers.hpp | 127 ++++++++++++++++++ SM-Coppens-1.m => Octave-tests/SM-Coppens-1.m | 0 SMMaponiA3.m => Octave-tests/SMMaponiA3.m | 0 SMMaponiA4.m => Octave-tests/SMMaponiA4.m | 0 SM-MaponiA3.cpp | 5 + SM-MaponiA3.hpp | 126 +---------------- main.cpp | 27 ++-- 7 files changed, 145 insertions(+), 140 deletions(-) create mode 100644 Helpers.hpp rename SM-Coppens-1.m => Octave-tests/SM-Coppens-1.m (100%) rename SMMaponiA3.m => Octave-tests/SMMaponiA3.m (100%) rename SMMaponiA4.m => Octave-tests/SMMaponiA4.m (100%) diff --git a/Helpers.hpp b/Helpers.hpp new file mode 100644 index 0000000..a81e38a --- /dev/null +++ b/Helpers.hpp @@ -0,0 +1,127 @@ +// Helpers.hpp +// Some usefull helper functions to support the Maponi algorithm. +#include +#include +#include +using namespace std; + +template +unsigned int getMaxIndex(T *vector, unsigned int size) { + unsigned int i; + unsigned int max = vector[0]; + unsigned int maxi = 0; + for (i = 1; i < size; i++) { + if (vector[i] > max) { + max = vector[i]; + maxi = i; + } + } + return maxi; +} + +template +void showScalar(T scalar, string name) { + cout << name << " = " << scalar << endl << endl; +} + +template +void showVector(T* vector, unsigned int size, string name) { + cout << name << " = " << endl; + for (unsigned int i = 0; i < size; i++) { + cout << "[ " << vector[i] << " ]" << endl; + } + cout << endl; +} + +template +void showMatrix(T** matrix, unsigned int size, string name) { + cout << name << " = " << endl; + for (unsigned int i = 0; i < size; i++) { + cout << "[ "; + for (unsigned int j = 0; j < size; j++) { + cout << matrix[i][j] << " "; + } + cout << " ]" << endl; + } + cout << endl; +} + +template +void showMatrixT(T** matrix, unsigned int size, string name) { + cout << name << " = " << endl; + for (unsigned int i = 0; i < size; i++) { + cout << "[ "; + for (unsigned int j = 0; j < size; j++) { + cout << matrix[j][i] << " "; + } + cout << " ]" << endl; + } + cout << endl; +} + +template +T** matMul(T** A, T** B, unsigned int size) { + T** C = new T*[size]; + for (unsigned int i = 0; i < size; i++) { + C[i] = new T[size]; + } + for (unsigned int i = 0; i < size; i++) { + for (unsigned int j = 0; j < size; j++) { + for (unsigned int k = 0; k < size; k++) { + C[i][j] += A[i][k] * B[k][j]; + } + } + } + showMatrix(C, size, "C"); + return C; +} + +template +T1** outProd(T1* vec1, T2* vec2, unsigned int size) { + T1** C = new T1*[size]; + for (unsigned int i = 0; i < size; i++) { + C[i] = new T1[size]; + } + for (unsigned int i = 0; i < size; i++) { + for (unsigned int j = 0; j < size; j++) { + C[i][j] = vec1[i+1] * vec2[j]; + } + } + 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; + 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; + } + delete [] temp; +} \ No newline at end of file diff --git a/SM-Coppens-1.m b/Octave-tests/SM-Coppens-1.m similarity index 100% rename from SM-Coppens-1.m rename to Octave-tests/SM-Coppens-1.m diff --git a/SMMaponiA3.m b/Octave-tests/SMMaponiA3.m similarity index 100% rename from SMMaponiA3.m rename to Octave-tests/SMMaponiA3.m diff --git a/SMMaponiA4.m b/Octave-tests/SMMaponiA4.m similarity index 100% rename from SMMaponiA4.m rename to Octave-tests/SMMaponiA4.m diff --git a/SM-MaponiA3.cpp b/SM-MaponiA3.cpp index abae3a7..207682a 100644 --- a/SM-MaponiA3.cpp +++ b/SM-MaponiA3.cpp @@ -1,5 +1,8 @@ // SM-MaponiA3.cpp +// Algorithm 3 from P. Maponi, +// p. 283, doi:10.1016/j.laa.2006.07.007 #include "SM-MaponiA3.hpp" +#include "Helpers.hpp" void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index) { unsigned int k, l, lbar, i, j, tmp, M = *Dim; @@ -83,7 +86,9 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns Al[i][j] = Id[i][j] - U[i][j] / beta; } } + showMatrix(Slater_inv, M, "Slater_inv"); Slater_inv = matMul(Al, Slater_inv, M); + showMatrix(Slater_inv, M, "Slater_inv"); } delete [] p, breakdown; diff --git a/SM-MaponiA3.hpp b/SM-MaponiA3.hpp index 7a501ec..9b5a48f 100644 --- a/SM-MaponiA3.hpp +++ b/SM-MaponiA3.hpp @@ -1,127 +1,3 @@ // SM-MaponiA3.hpp -#include -#include -#include -using namespace std; -void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index); - -template -unsigned int getMaxIndex(T *vector, unsigned int size) { - unsigned int i; - unsigned int max = vector[0]; - unsigned int maxi = 0; - for (i = 1; i < size; i++) { - if (vector[i] > max) { - max = vector[i]; - maxi = i; - } - } - return maxi; -} - -template -void showScalar(T scalar, string name) { - cout << name << " = " << scalar << endl << endl; -} - -template -void showVector(T* vector, unsigned int size, string name) { - cout << name << " = " << endl; - for (unsigned int i = 0; i < size; i++) { - cout << "[ " << vector[i] << " ]" << endl; - } - cout << endl; -} - -template -void showMatrix(T** matrix, unsigned int size, string name) { - cout << name << " = " << endl; - for (unsigned int i = 0; i < size; i++) { - cout << "[ "; - for (unsigned int j = 0; j < size; j++) { - cout << matrix[i][j] << " "; - } - cout << " ]" << endl; - } - cout << endl; -} - -template -void showMatrixT(T** matrix, unsigned int size, string name) { - cout << name << " = " << endl; - for (unsigned int i = 0; i < size; i++) { - cout << "[ "; - for (unsigned int j = 0; j < size; j++) { - cout << matrix[j][i] << " "; - } - cout << " ]" << endl; - } - cout << endl; -} - -template -T** matMul(T** A, T** B, unsigned int size) { - T** C = new T*[size]; - for (unsigned int i = 0; i < size; i++) { - C[i] = new T[size]; - } - for (unsigned int i = 0; i < size; i++) { - for (unsigned int j = 0; j < size; j++) { - for (unsigned int k = 0; k < size; k++) { - C[i][j] += A[i][k] * B[k][j]; - } - } - } - return C; -} - -template -T1** outProd(T1* vec1, T2* vec2, unsigned int size) { - T1** C = new T1*[size]; - for (unsigned int i = 0; i < size; i++) { - C[i] = new T1[size]; - } - for (unsigned int i = 0; i < size; i++) { - for (unsigned int j = 0; j < size; j++) { - C[i][j] = vec1[i+1] * vec2[j]; - } - } - 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; - 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; - } - delete [] temp; -} \ No newline at end of file +void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index); \ No newline at end of file diff --git a/main.cpp b/main.cpp index eec6795..7853531 100644 --- a/main.cpp +++ b/main.cpp @@ -1,7 +1,6 @@ -// Algorithm 3 from P. Maponi, -// p. 283, doi:10.1016/j.laa.2006.07.007 - +// main.cpp #include "SM-MaponiA3.hpp" +#include "Helpers.hpp" #include #include @@ -12,7 +11,7 @@ int main() { unsigned int M = 3; // Dimension of the Slater-matrix unsigned int i, j; // Indices for iterators - // Declare and allocate all vectors and matrices + // Declare, allocate all vectors and matrices and fill them with zeros unsigned int *Ar_index = new unsigned int[M]; int **A = new int*[M]; // The matrix to be inverted int **A0 = new int*[M]; // A diagonal matrix with the digonal elements of A @@ -24,8 +23,7 @@ int main() { Ar[i] = new int[M]; A0_inv[i] = new double[M]; } - - // Initialize all matrices with zeros + // Fill with zeros for (i = 0; i < M; i++) { for (j = 0; j < M; j++) { A0[i][j] = 0; @@ -34,11 +32,10 @@ int main() { } } - // Initialize A with M=3 and Eq. (17) from paper + // Initialize A with M=3 and fill acc. to Eq. (17) from paper A[0][0] = 1; A[0][1] = 1; A[0][2] = -1; A[1][0] = 1; A[1][1] = 1; A[1][2] = 0; A[2][0] = -1; A[2][1] = 0; A[2][2] = -1; - // // Fill A with random numbers from [-randRange,randRange] // // and check if A and A0 are invertable // do { @@ -51,11 +48,13 @@ int main() { // A0[i][i] = A[i][i]; // } // } while (matDet(A, M) == 0 || matDet(A0, M) == 0); - showMatrix(A, M, "A"); - // Init the update matrix Ar, A0_inv and Ar_index + // 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][i] = A[i][i]; A0_inv[i][i] = 1.0/A[i][i]; Ar_index[i] = i; for (j = 0; j < M; j++) { @@ -63,6 +62,7 @@ int main() { } } + // Define pointers dim and n_updates to use in Sherman-Morrison(...) function call unsigned int *dim = new unsigned int(M); unsigned int *n_updates = new unsigned int(M); Sherman_Morrison(A0, A0_inv, dim, n_updates, Ar, Ar_index); @@ -71,13 +71,10 @@ int main() { // Deallocate all vectors and matrices for (i = 0; i < M; i++) { - delete [] A[i]; - delete [] A0[i]; - delete [] A0_inv[i]; - delete [] Ar[i]; + delete [] A[i], A0[i], A0_inv[i], Ar[i]; } delete [] A, A0, A0_inv, Ar, Ar_index; delete dim, n_updates; return 0; -} +} \ No newline at end of file