diff --git a/.gitignore b/.gitignore index f13e716..b0ac550 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ *.o SM-MaponiA3 +Sherman-Morrison .vscode \ No newline at end of file diff --git a/SM-MaponiA3.cpp b/SM-MaponiA3.cpp index f4209b9..abae3a7 100644 --- a/SM-MaponiA3.cpp +++ b/SM-MaponiA3.cpp @@ -1,128 +1,9 @@ // SM-MaponiA3.cpp #include "SM-MaponiA3.hpp" -uint getMaxIndex(double *arr, uint size) { - uint i; - uint max = arr[0]; - uint maxi = 0; - for (i = 1; i < size; i++) { - if (arr[i] > max) { - max = arr[i]; - maxi = i; - } - } - return maxi; -} - -template -void showScalar(T scalar, string name) { - cout << name << " = " << scalar << endl << endl; -} - -template -void showVector(T* vector, uint size, string name) { - cout << name << " = " << endl; - for (uint i = 0; i < size; i++) { - cout << "[ " << vector[i] << " ]" << endl; - } - cout << endl; -} - -template -void showMatrix(T** matrix, uint size, string name) { - cout << name << " = " << endl; - for (uint i = 0; i < size; i++) { - cout << "[ "; - for (uint j = 0; j < size; j++) { - cout << matrix[i][j] << " "; - } - cout << " ]" << endl; - } - cout << endl; -} - -template -void showMatrixT(T** matrix, uint size, string name) { - cout << name << " = " << endl; - for (uint i = 0; i < size; i++) { - cout << "[ "; - for (uint j = 0; j < size; j++) { - cout << matrix[j][i] << " "; - } - cout << " ]" << endl; - } - cout << endl; -} - -template -T** matMul(T** A, T** B, uint size) { - T** C = new T*[size]; - for (uint i = 0; i < size; i++) { - C[i] = new T[size]; - } - for (uint i = 0; i < size; i++) { - for (uint j = 0; j < size; j++) { - for (uint k = 0; k < size; k++) { - C[i][j] += A[i][k] * B[k][j]; - } - } - } - return C; -} - -template -T1** outProd(T1* vec1, T2* vec2, uint size) { - T1** C = new T1*[size]; - for (uint i = 0; i < size; i++) { - C[i] = new T1[size]; - } - for (uint i = 0; i < size; i++) { - for (uint j = 0; j < size; j++) { - C[i][j] = vec1[i+1] * vec2[j]; - } - } - return C; -} - -template -T matDet(T** A, 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; -} - -void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index) { - uint k, l, lbar, i, j, tmp, M = *Dim; - uint *p = new uint[M+1]; +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]; double alpha, beta; @@ -219,4 +100,4 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd } delete [] ylk[l]; } -} +} \ No newline at end of file diff --git a/SM-MaponiA3.hpp b/SM-MaponiA3.hpp index e5b742a..7a501ec 100644 --- a/SM-MaponiA3.hpp +++ b/SM-MaponiA3.hpp @@ -4,12 +4,124 @@ #include using namespace std; -uint getMaxIndex(double *arr, uint size); -templatevoid showScalar(T scalar, string name); -templatevoid showVector(T *vector, uint size, string name); -templatevoid showMatrix(T **matrix, uint size, string name); -templatevoid showMatrixT(T **matrix, uint size, string name); -templateT **matMul(T **A, T **B, uint size); -templateT1 **outProd(T1 *vec1, T2 *vec2, uint size); -templateT matDet(T **A, int M); -void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index); +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 diff --git a/main.cpp b/main.cpp index 3889f9e..eec6795 100644 --- a/main.cpp +++ b/main.cpp @@ -1,20 +1,19 @@ // Algorithm 3 from P. Maponi, // p. 283, doi:10.1016/j.laa.2006.07.007 +#include "SM-MaponiA3.hpp" #include #include -#include "SM-MaponiA3.hpp" -using namespace std; int main() { srand((unsigned) time(0)); - uint randRange = 1; // to get random integers in range [-randRange, randRange] - uint M = 3; // Dimension of the Slater-matrix - uint i, j; // Indices for iterators + unsigned int randRange = 1; // to get random integers in range [-randRange, randRange] + unsigned int M = 3; // Dimension of the Slater-matrix + unsigned int i, j; // Indices for iterators // Declare and allocate all vectors and matrices - uint *Ar_index = new uint[M]; + 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 int **Ar = new int*[M]; // The update matrix @@ -64,8 +63,8 @@ int main() { } } - uint *dim = new uint(M); - uint *n_updates = new uint(M); + 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");