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