From 9d47bf9bb30b7cb10806ed65082a2ca061d292e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Tue, 2 Feb 2021 15:05:47 +0100 Subject: [PATCH] Extracted the SM algorithm and wrote it as a void function that updates the passed old inverse. Does not work yet. For some reason the changes in the inverse do not persist when the function is finished. Also the inverse is not correct. --- SM-MaponiA3.cpp | 287 ++++++++++++++++++++++++++---------------------- 1 file changed, 154 insertions(+), 133 deletions(-) diff --git a/SM-MaponiA3.cpp b/SM-MaponiA3.cpp index d6c6346..85fec09 100644 --- a/SM-MaponiA3.cpp +++ b/SM-MaponiA3.cpp @@ -8,160 +8,87 @@ #include using namespace std; -uint getMaxIndex(double* arr, uint size); +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); +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; - uint i, j, k, l, lbar, tmp; - double alpha, beta; + uint M = 3; // Dimension of the Slater-matrix + uint i, j; // Indices for iterators // Declare and allocate all vectors and matrices - uint* p = new uint[M+1]; - double* breakdown = new double[M+1]; - - int** A = new int*[M]; - int** A0 = new int*[M]; - int** Ar = new int*[M]; - int** Id = new int*[M]; - double** A0inv = new double*[M]; - double** Ainv; //= new double*[M]; + 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]; - //Ainv[i] = new double[M]; A0[i] = new int[M]; - A0inv[i] = new double[M]; Ar[i] = new int[M]; - Id[i] = new int[M]; - } - - double*** ylk = new double**[M]; - for (l = 0; l < M; l++) { - ylk[l] = new double*[M+1]; - for (k = 0; k < M+1; k++) { - ylk[l][k] = new double[M+1]; - } + 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; - A0inv[i][j] = 0; Ar[i][j] = 0; - Id[i][j] = 0; + A0inv[i][j] = 0; } } - // Initialize ylk with zeros - for (l = 0; l < M; l++) { - for (k = 0; k < M+1; k++) { - for (i = 0; i < M+1; i++) { - ylk[l][k][i] = 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; - // // 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"); - // Fill A with random numbers from [-1,1] - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - A[i][j] = rand()%(2*randRange+1)-randRange; - } - } - - // Define identity matrix, A0, A0inv and p - p[0] = 0; - for (i = 0; i < M; i++) { - Id[i][i] = 1; - A0[i][i] = A[i][i]; - A0inv[i][i] = 1.0/A[i][i]; - p[i+1] = i+1; - } - - // Init Ar + // 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]; } } - showMatrix(A, M, "A"); - // showMatrix(A0, M, M, "A0"); - // showMatrix(Ar, M, M, "Ar"); - // showMatrix(A0inv, M, M, "A0inv"); - - // Calculate all the y0k in M^2 multiplications instead of M^3 - for (k = 1; k < M+1; k++) { - for (i = 1; i < M+1; i++) { - ylk[0][k][i] = A0inv[i-1][i-1] * Ar[i-1][k-1]; - } - // showVector(ylk[0][k], M+1, "y0k"); + // Define A0inv + for (i = 0; i < M; i++) { + A0inv[i][i] = 1.0/A[i][i]; } - // showMatrixT(ylk[0], M+1, "y0k"); + showMatrix(A0inv, M, "A0inv"); - // Calculate all the ylk from the y0k - // showVector(p, M+1, "p"); - for (l = 1; l < M; l++) { - for (j = l; j < M+1; j++) { - breakdown[j] = abs( 1 + ylk[l-1][p[j]][p[j]] ); - } - // showVector(breakdown, M+1, "break-down vector"); - lbar = getMaxIndex(breakdown, M+1); - // showScalar(lbar, "lbar"); - for (i = 0; i < M; i++) { - breakdown[i] = 0; - } - tmp = p[l]; - p[l] = p[lbar]; - p[lbar] = tmp; - // showVector(p, M+1, "p"); - for (k = l+1; k < M+1; k++) { - beta = 1 + ylk[l-1][p[l]][p[l]]; - if (beta == 0) { - cout << "Break-down condition occured. Exiting..." << endl; - exit; - } - for (i = 1; i < M+1; i++) { - alpha = ylk[l-1][p[k]][p[l]] / beta; - ylk[l][p[k]][i] = ylk[l-1][p[k]][i] - alpha * ylk[l-1][p[l]][i]; - } - // showVector(ylk[l][p[k]], M+1, "ylk"); - } + // Populate indices_of_updates + for (i = 0; i < M; i++) { + indices_of_updates[i] = i; } - // showMatrixT(ylk[1], M+1, "y1k"); - // showMatrixT(ylk[2], M+1, "y2k"); - // Construct A-inverse from A0-inverse and the ylk - double** U; - double** Ap = new double*[M]; - for (i = 0; i < M; i++) Ap[i] = new double[M]; - Ainv = A0inv; - for (l = 0; l < M; l++) { - k = l+1; - U = outProd(ylk[l][p[k]], Id[p[k]-1], M); - beta = 1 + ylk[l][p[k]][p[k]]; - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - Ap[i][j] = Id[i][j] - U[i][j] / beta; - } - } - Ainv = matMul(Ap, Ainv, M); - } - showMatrixT(Ainv, M, "Ainv"); + 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++) { @@ -169,24 +96,14 @@ int main() { delete [] A0[i]; delete [] A0inv[i]; delete [] Ar[i]; - delete [] Id[i]; - delete [] U[i]; - delete [] Ap[i]; } - for (l = 0; l < M; l++) { - for (k = 0; k < M+1; k++) { - delete [] ylk[l][k]; - } - delete [] ylk[l]; - } - - delete [] p, breakdown, A, Ainv, A0, A0inv, Ar, Id, ylk; + delete [] A, A0, A0inv, Ar; return 0; } -uint getMaxIndex(double* arr, uint size) { +uint getMaxIndex(double *arr, uint size) { uint i; uint max = arr[0]; uint maxi = 0; @@ -274,7 +191,7 @@ 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) { + if(M == 1) { return A[0][0]; } else if(M == 2) { @@ -303,4 +220,108 @@ T matDet(T** A, int M) { return det; } delete [] temp; +} + +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]; + 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]; + for (i = 0; i < M; i++) { + for (j = 0; j < M; j++) { + if (i != j) Id[i][j] = 0; + else Id[i][j] = 1; + } + } + + // Declare auxiliary solution matrix ylk + double ***ylk = new double**[M]; + for (l = 0; l < M; l++) { + ylk[l] = new double*[M+1]; + for (k = 0; k < M+1; k++) { + ylk[l][k] = new double[M+1]; + } + } + // Initialize ylk with zeros + for (l = 0; l < M; l++) { + for (k = 0; k < M+1; k++) { + for (i = 0; i < M+1; i++) { + ylk[l][k][i] = 0; + } + } + } + + // Calculate all the y0k in M^2 multiplications instead of M^3 + for (k = 1; k < M+1; k++) { + for (i = 1; i < M+1; i++) { + ylk[0][k][i] = Slater_inv[i-1][i-1] * Updates[i-1][k-1]; + } + } + + // Calculate all the ylk from the y0k + for (l = 1; l < M; l++) { + for (j = l; j < M+1; j++) { + breakdown[j] = abs( 1 + ylk[l-1][p[j]][p[j]] ); + } + lbar = getMaxIndex(breakdown, M+1); + for (i = 0; i < M; i++) { + breakdown[i] = 0; + } + tmp = p[l]; + p[l] = p[lbar]; + p[lbar] = tmp; + for (k = l+1; k < M+1; k++) { + beta = 1 + ylk[l-1][p[l]][p[l]]; + if (beta == 0) { + cout << "Break-down condition occured. Exiting..." << endl; + exit; + } + for (i = 1; i < M+1; i++) { + alpha = ylk[l-1][p[k]][p[l]] / beta; + ylk[l][p[k]][i] = ylk[l-1][p[k]][i] - alpha * ylk[l-1][p[l]][i]; + } + } + } + + 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]; + 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); + beta = 1 + ylk[l][p[k]][p[k]]; + for (i = 0; i < M; i++) { + for (j = 0; j < M; j++) { + Al[i][j] = Id[i][j] - U[i][j] / beta; + } + } + Slater_inv = matMul(Al, Slater_inv, M); + } + cout << "Just after the final construction of the inverse..." << endl; + showMatrix(Slater_inv, M, "Slater_inv"); + + for (i = 0; i < M; i++) { + delete [] U[i]; + delete [] Al[i]; + } + + for (l = 0; l < M; l++) { + for (k = 0; k < M+1; k++) { + delete [] ylk[l][k]; + } + delete [] ylk[l]; + } } \ No newline at end of file