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

@@ -8,160 +8,87 @@
 #include <ctime>
 using namespace std;
 
-uint getMaxIndex(double* arr, uint size);
+uint getMaxIndex(double *arr, uint size);
 template<typename T>void showScalar(T scalar, string name);
-template<typename T>void showVector(T* vector, uint size, string name);
+template<typename T>void showVector(T *vector, uint size, string name);
-template<typename T>void showMatrix(T** matrix, uint size, string name);
+template<typename T>void showMatrix(T **matrix, uint size, string name);
-template<typename T>void showMatrixT(T** matrix, uint size, string name);
+template<typename T>void showMatrixT(T **matrix, uint size, string name);
-template<typename T>T** matMul(T** A, T** B, uint size);
+template<typename T>T **matMul(T **A, T **B, uint size);
-template<typename T1, typename T2>T1** outProd(T1* vec1, T2* vec2, uint size);
+template<typename T1, typename T2>T1 **outProd(T1 *vec1, T2 *vec2, uint size);
-template<typename T>T matDet(T** A, int M);
+template<typename T>T 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 M = 3; // Dimension of the Slater-matrix
-    uint i, j, k, l, lbar, tmp;
+    uint i, j; // Indices for iterators
-    double alpha, beta;
 
     // Declare and allocate all vectors and matrices
-    uint* p = new uint[M+1];
+    uint *indices_of_updates = new uint[M];
-    double* breakdown = new double[M+1];
+    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** A = new int*[M];
+    int **Ar = new int*[M]; // The update matrix
-    int** A0 = new int*[M];
+    double **A0inv = new double*[M]; // Inverse of A0
-    int** Ar = new int*[M];
-    int** Id = new int*[M];
-    double** A0inv = new double*[M];
-    double** Ainv; //= new double*[M];
     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];
+        A0inv[i] = new double[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];
-        }
     }
 
     // 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
+    // Initialize A with M=3 and Eq. (17) from paper
-    for (l = 0; l < M; l++) {
+    A[0][0] = 1;    A[0][1] = 1;    A[0][2] = -1;
-        for (k = 0; k < M+1; k++) {
+    A[1][0] = 1;    A[1][1] = 1;    A[1][2] = 0;
-            for (i = 0; i < M+1; i++) {
+    A[2][0] = -1;   A[2][1] = 0;    A[2][2] = -1;
-                ylk[l][k][i] = 0;
-            }
-        }
-    }
 
-    // // Initialize A with M=3 and Eq. (17) from paper
+    // // Fill A with random numbers from [-randRange,randRange] 
-    // A[0][0] = 1;    A[0][1] = 1;    A[0][2] = -1;
+    // // and check if A and A0 are invertable
-    // A[1][0] = 1;    A[1][1] = 1;    A[1][2] = 0;
+    // do {
-    // A[2][0] = -1;   A[2][1] = 0;    A[2][2] = -1;
+    //     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);
  
-    // Fill A with random numbers from [-1,1]
+    showMatrix(A, M, "A");
-    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
+    // Init Ar: the update matrix
-    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
     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");
+    // Define A0inv
-    // showMatrix(A0, M, M, "A0");
+    for (i = 0; i < M; i++) {
-    // showMatrix(Ar, M, M, "Ar");
+        A0inv[i][i] = 1.0/A[i][i];
-    // 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");
     }
-    // showMatrixT(ylk[0], M+1, "y0k");
+    showMatrix(A0inv, M, "A0inv");
 
-    // Calculate all the ylk from the y0k
+    // Populate indices_of_updates
-    // showVector(p, M+1, "p");
+    for (i = 0; i < M; i++) {
-    for (l = 1; l < M; l++) {
+        indices_of_updates[i] = i;
-        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");
-        }
     }
-    // showMatrixT(ylk[1], M+1, "y1k");
-    // showMatrixT(ylk[2], M+1, "y2k");
 
-    // Construct A-inverse from A0-inverse and the ylk
+    uint *dim = new uint(M);
-    double** U;
+    Sherman_Morrison(A0, A0inv, dim, dim, Ar, indices_of_updates);
-    double** Ap = new double*[M];
+    
-    for (i = 0; i < M; i++) Ap[i] = new double[M];
+    showMatrixT(A0inv, M, "A0inv");
-    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");
     
     // 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++) {
+    delete [] A, A0, A0inv, Ar;
-        for (k = 0; k < M+1; k++) {
-            delete [] ylk[l][k];
-        }
-        delete [] ylk[l];
-    }
-
-    delete [] p, breakdown, A, Ainv, A0, A0inv, Ar, Id, ylk;
 
     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) {
@@ -304,3 +221,107 @@ T matDet(T** A, int M) {
     }
     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];
+    }
+}