diff --git a/Helpers.hpp b/Helpers.hpp new file mode 100644 index 0000000..646ba6e --- /dev/null +++ b/Helpers.hpp @@ -0,0 +1,126 @@ +// 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]; + } + } + } + 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..df94ef2 100644 --- a/SM-MaponiA3.cpp +++ b/SM-MaponiA3.cpp @@ -1,23 +1,21 @@ // 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; 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 @@ -28,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++) { @@ -70,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,18 +94,18 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns Slater_inv = matMul(Al, Slater_inv, M); } - 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/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..55af456 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,21 +62,18 @@ 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); - showMatrix(A0_inv, M, "A0_inv"); // 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