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 1/2] 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 From d03c9c665b9e8fbdf86759f32435dc8e61441857 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Wed, 3 Feb 2021 14:49:13 +0100 Subject: [PATCH 2/2] Update of A0_inv fixed. matMul(...) breaks link between A0_inv and Slater_inv. Fixed by creating a new pointer that points to where Slater_inv points to before it gets reassigned by matMul(...) and then copy the updated elements back to the passed array pointed to by the new pointer. --- Helpers.hpp | 19 +++++++++--------- SM-MaponiA3.cpp | 51 ++++++++++++++++++++++++++----------------------- main.cpp | 1 - 3 files changed, 36 insertions(+), 35 deletions(-) diff --git a/Helpers.hpp b/Helpers.hpp index a81e38a..646ba6e 100644 --- a/Helpers.hpp +++ b/Helpers.hpp @@ -25,7 +25,7 @@ void showScalar(T scalar, string name) { } template -void showVector(T* vector, unsigned int size, string name) { +void showVector(T *vector, unsigned int size, string name) { cout << name << " = " << endl; for (unsigned int i = 0; i < size; i++) { cout << "[ " << vector[i] << " ]" << endl; @@ -34,7 +34,7 @@ void showVector(T* vector, unsigned int size, string name) { } template -void showMatrix(T** matrix, unsigned int size, string name) { +void showMatrix(T **matrix, unsigned int size, string name) { cout << name << " = " << endl; for (unsigned int i = 0; i < size; i++) { cout << "[ "; @@ -47,7 +47,7 @@ void showMatrix(T** matrix, unsigned int size, string name) { } template -void showMatrixT(T** matrix, unsigned int size, string name) { +void showMatrixT(T **matrix, unsigned int size, string name) { cout << name << " = " << endl; for (unsigned int i = 0; i < size; i++) { cout << "[ "; @@ -60,8 +60,8 @@ void showMatrixT(T** matrix, unsigned int size, string name) { } template -T** matMul(T** A, T** B, unsigned int size) { - T** C = new T*[size]; +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]; } @@ -72,13 +72,12 @@ T** matMul(T** A, T** B, unsigned int size) { } } } - showMatrix(C, size, "C"); return C; } template -T1** outProd(T1* vec1, T2* vec2, unsigned int size) { - T1** C = new T1*[size]; +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]; } @@ -91,9 +90,9 @@ T1** outProd(T1* vec1, T2* vec2, unsigned int size) { } template -T matDet(T** A, unsigned int M) { +T matDet(T **A, unsigned int M) { int det = 0, p, h, k, i, j; - T** temp = new T*[M]; + T **temp = new T*[M]; for (int i = 0; i < M; i++) temp[i] = new T[M]; if(M == 1) { return A[0][0]; diff --git a/SM-MaponiA3.cpp b/SM-MaponiA3.cpp index 207682a..df94ef2 100644 --- a/SM-MaponiA3.cpp +++ b/SM-MaponiA3.cpp @@ -7,20 +7,15 @@ 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; unsigned int *p = new unsigned int[M+1]; - double *breakdown = new double[M+1]; + unsigned int **Id = new unsigned int*[M]; double alpha, beta; - - for (i = 0; i < M+1; i++) { - p[i] = i; - } - - int **Id = new int*[M]; - for (i = 0; i < M; i++) Id[i] = new int[M]; + double **U, *breakdown = new double[M+1]; + double **Al = new double*[M]; + p[0] = 0; for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - if (i != j) Id[i][j] = 0; - else Id[i][j] = 1; - } + p[i+1] = i + 1; + Id[i] = new unsigned int[M]; + Al[i] = new double[M]; } // Declare auxiliary solution matrix ylk @@ -31,6 +26,15 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns ylk[l][k] = new double[M+1]; } } + + // Initialize identity matrix + for (i = 0; i < M; i++) { + for (j = 0; j < M; j++) { + if (i != j) Id[i][j] = 0; + else Id[i][j] = 1; + } + } + // Initialize ylk with zeros for (l = 0; l < M; l++) { for (k = 0; k < M+1; k++) { @@ -73,10 +77,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns } // Construct A-inverse from A0-inverse and the ylk - double **U; - double **Al = new double*[M]; - for (i = 0; i < M; i++) Al[i] = new double[M]; - + // 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. + double **copy = Slater_inv; + for (l = 0; l < M; l++) { k = l+1; U = outProd(ylk[l][p[k]], Id[p[k]-1], M); @@ -86,23 +91,21 @@ 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; - + // Assign the new values of 'Slater_inv' to the old values in 'copy[][]' for (i = 0; i < M; i++) { - delete [] Id[i]; - delete [] U[i]; - delete [] Al[i]; + for (j = 0; j < M; j++) { + copy[i][j] = Slater_inv[i][j]; + } } for (l = 0; l < M; l++) { for (k = 0; k < M+1; k++) { delete [] ylk[l][k]; } - delete [] ylk[l]; + delete [] ylk[l], Id[l], U[l], Al[l], Slater_inv[l]; } + delete [] p, breakdown; } \ No newline at end of file diff --git a/main.cpp b/main.cpp index 7853531..55af456 100644 --- a/main.cpp +++ b/main.cpp @@ -66,7 +66,6 @@ int main() { 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); - showMatrix(A0_inv, M, "A0_inv"); // Deallocate all vectors and matrices