diff --git a/Makefile b/Makefile index ccb052d..24460e3 100644 --- a/Makefile +++ b/Makefile @@ -5,12 +5,12 @@ CXXFLAGS=-O0 -debug full -traceback # ARCH=-xCORE-AVX2 DEPS = SM-MaponiA3.cpp -OBJ = SM-MaponiA3.o +OBJ = SM-MaponiA3.o main.o -%.o: %.c $(DEPS) +%.o: %.cpp $(DEPS) $(CXX) $(ARCH) -c -o $@ $< $(CFLAGS) -SM-MaponiA3: $(OBJ) +Sherman-Morrison: $(OBJ) $(CXX) $(ARCH) -o $@ $^ $(CFLAGS) clean: diff --git a/SM-MaponiA3.cpp b/SM-MaponiA3.cpp index 85fec09..f4209b9 100644 --- a/SM-MaponiA3.cpp +++ b/SM-MaponiA3.cpp @@ -1,107 +1,5 @@ -// Algorithm 3 from P. Maponi, -// p. 283, doi:10.1016/j.laa.2006.07.007 - -#include -#include -#include -#include -#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); - -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 - - // Declare and allocate all vectors and matrices - uint *indices_of_updates = new uint[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 - double **A0inv = new double*[M]; // Inverse of A0 - for (i = 0; i < M; i++) { - A[i] = new int[M]; - A0[i] = new int[M]; - Ar[i] = new int[M]; - A0inv[i] = new double[M]; - } - - // Initialize all matrices with zeros - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - A0[i][j] = 0; - Ar[i][j] = 0; - A0inv[i][j] = 0; - } - } - - // Initialize A with M=3 and 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 { - // for (i = 0; i < M; i++) { - // for (j = 0; j < M; j++) { - // A[i][j] = rand()%(2*randRange+1)-randRange; - // } - // } - // for (i = 0; i < M; i++) { - // A0[i][i] = A[i][i]; - // } - // } while (matDet(A, M) == 0 || matDet(A0, M) == 0); - - showMatrix(A, M, "A"); - - // Init Ar: the update matrix - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - Ar[i][j] = A[i][j] - A0[i][j]; - } - } - - // Define A0inv - for (i = 0; i < M; i++) { - A0inv[i][i] = 1.0/A[i][i]; - } - showMatrix(A0inv, M, "A0inv"); - - // Populate indices_of_updates - for (i = 0; i < M; i++) { - indices_of_updates[i] = i; - } - - uint *dim = new uint(M); - Sherman_Morrison(A0, A0inv, dim, dim, Ar, indices_of_updates); - - showMatrixT(A0inv, M, "A0inv"); - - // Deallocate all vectors and matrices - for (i = 0; i < M; i++) { - delete [] A[i]; - delete [] A0[i]; - delete [] A0inv[i]; - delete [] Ar[i]; - } - - delete [] A, A0, A0inv, Ar; - - return 0; -} +// SM-MaponiA3.cpp +#include "SM-MaponiA3.hpp" uint getMaxIndex(double *arr, uint size) { uint i; @@ -223,8 +121,6 @@ T matDet(T** A, int M) { } void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index) { - cout << "We just entered Sherman-Morrison" << endl; - showMatrix(Slater_inv, *Dim, "Slater_inv"); uint k, l, lbar, i, j, tmp, M = *Dim; uint *p = new uint[M+1]; double *breakdown = new double[M+1]; @@ -292,13 +188,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd } } - cout << "Just before the final construction of the inverse..." << endl; - showMatrix(Slater_inv, M, "Slater_inv"); - // Construct A-inverse from A0-inverse and the ylk - double** U; - double** Al = new double*[M]; + double **U; + double **Al = new double*[M]; for (i = 0; i < M; i++) Al[i] = new double[M]; + for (l = 0; l < M; l++) { k = l+1; U = outProd(ylk[l][p[k]], Id[p[k]-1], M); @@ -310,10 +204,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd } Slater_inv = matMul(Al, Slater_inv, M); } - cout << "Just after the final construction of the inverse..." << endl; - showMatrix(Slater_inv, M, "Slater_inv"); + + delete [] p, breakdown; for (i = 0; i < M; i++) { + delete [] Id[i]; delete [] U[i]; delete [] Al[i]; } @@ -324,4 +219,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 new file mode 100644 index 0000000..e5b742a --- /dev/null +++ b/SM-MaponiA3.hpp @@ -0,0 +1,15 @@ +// SM-MaponiA3.hpp +#include +#include +#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); diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..3889f9e --- /dev/null +++ b/main.cpp @@ -0,0 +1,84 @@ +// Algorithm 3 from P. Maponi, +// p. 283, doi:10.1016/j.laa.2006.07.007 + +#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 + + // Declare and allocate all vectors and matrices + uint *Ar_index = new uint[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 + double **A0_inv = new double*[M]; // Inverse of A0 + for (i = 0; i < M; i++) { + A[i] = new int[M]; + A0[i] = new int[M]; + Ar[i] = new int[M]; + A0_inv[i] = new double[M]; + } + + // Initialize all matrices with zeros + for (i = 0; i < M; i++) { + for (j = 0; j < M; j++) { + A0[i][j] = 0; + Ar[i][j] = 0; + A0_inv[i][j] = 0; + } + } + + // Initialize A with M=3 and 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 { + // for (i = 0; i < M; i++) { + // for (j = 0; j < M; j++) { + // A[i][j] = rand()%(2*randRange+1)-randRange; + // } + // } + // for (i = 0; i < M; i++) { + // 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 + for (i = 0; i < M; i++) { + A0_inv[i][i] = 1.0/A[i][i]; + Ar_index[i] = i; + for (j = 0; j < M; j++) { + Ar[i][j] = A[i][j] - A0[i][j]; + } + } + + uint *dim = new uint(M); + uint *n_updates = new uint(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, A0, A0_inv, Ar, Ar_index; + delete dim, n_updates; + + return 0; +}