Merge pull request #2 from fmgjcoppens/fix/upd-inverse

Fix/upd inverse

Helpers.hpp

@@ -0,0 +1,126 @@
+// Helpers.hpp
+// Some usefull helper functions to support the Maponi algorithm.
+#include <iostream>
+#include <cmath>
+#include <string>
+using namespace std;
+
+template<typename T>
+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<typename T>
+void showScalar(T scalar, string name) {
+    cout << name << " = " << scalar << endl << endl;
+}
+
+template<typename T>
+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<typename T>
+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<typename T>
+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<typename T>
+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<typename T1, typename T2>
+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<typename T>
+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;
+}

SM-MaponiA3.cpp

@@ -1,23 +1,21 @@
 // SM-MaponiA3.cpp
+// Algorithm 3 from P. Maponi,
+// p. 283, doi:10.1016/j.laa.2006.07.007
 #include "SM-MaponiA3.hpp"
+#include "Helpers.hpp"
 
 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];
+    unsigned int **Id = new unsigned int*[M];
     double alpha, beta;
-
+    double **U, *breakdown = new double[M+1];
-    for (i = 0; i < M+1; i++) {
+    double **Al = new double*[M];
-        p[i] = i;
+    p[0] = 0;
-    }
-
-    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++) {
+        p[i+1] = i + 1;
-            if (i != j) Id[i][j] = 0;
+        Id[i] = new unsigned int[M];
-            else Id[i][j] = 1;
+        Al[i] = new double[M];
-        }
     }
 
     // Declare auxiliary solution matrix ylk
@@ -28,6 +26,15 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
             ylk[l][k] = new double[M+1];
         }
     }
+
+    // Initialize identity matrix
+    for (i = 0; i < M; i++) {
+        for (j = 0; j < M; j++) {
+            if (i != j) Id[i][j] = 0;
+            else Id[i][j] = 1;
+        }
+    }
+    
     // Initialize ylk with zeros
     for (l = 0; l < M; l++) {
         for (k = 0; k < M+1; k++) {
@@ -70,10 +77,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
     }
 
     // Construct A-inverse from A0-inverse and the ylk
-    double **U;
+    // Keep the memory location of the passed array 'Slater_inv' before 'Slater_inv'
-    double **Al = new double*[M];
+    // gets reassigned by 'matMul(...)' in the next line, by creating a new
-    for (i = 0; i < M; i++) Al[i] = new double[M];
+    // pointer 'copy' that points to whereever 'Slater_inv' points to now.
-  
+    double **copy = Slater_inv;
+
     for (l = 0; l < M; l++) {
         k = l+1;
         U = outProd(ylk[l][p[k]], Id[p[k]-1], M);
@@ -86,18 +94,18 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
         Slater_inv = matMul(Al, Slater_inv, M);
     }
 
-    delete [] p, breakdown;
+    // Assign the new values of 'Slater_inv' to the old values in 'copy[][]'
-
     for (i = 0; i < M; i++) {
-        delete [] Id[i];
+        for (j = 0; j < M; j++) {
-        delete [] U[i];
+            copy[i][j] = Slater_inv[i][j];
-        delete [] Al[i];
+        }
     }
 
     for (l = 0; l < M; l++) {
         for (k = 0; k < M+1; k++) {
             delete [] ylk[l][k];
         }
-        delete [] ylk[l];
+        delete [] ylk[l], Id[l], U[l], Al[l], Slater_inv[l];
     }
+    delete [] p, breakdown;
 }

SM-MaponiA3.hpp

@@ -1,127 +1,3 @@
 // SM-MaponiA3.hpp
-#include <iostream>
-#include <cmath>
-#include <string>
-using namespace std;
 
-void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index);
+void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index);
-
-template<typename T>
-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<typename T>
-void showScalar(T scalar, string name) {
-    cout << name << " = " << scalar << endl << endl;
-}
-
-template<typename T>
-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<typename T>
-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<typename T>
-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<typename T>
-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<typename T1, typename T2>
-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<typename T>
-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;
-}

main.cpp

@@ -1,7 +1,6 @@
-// Algorithm 3 from P. Maponi,
+// main.cpp
-// p. 283, doi:10.1016/j.laa.2006.07.007
-
 #include "SM-MaponiA3.hpp"
+#include "Helpers.hpp"
 #include <cstdlib>
 #include <ctime>
 
@@ -12,7 +11,7 @@ int main() {
     unsigned int M = 3; // Dimension of the Slater-matrix
     unsigned int i, j; // Indices for iterators
 
-    // Declare and allocate all vectors and matrices
+    // Declare, allocate all vectors and matrices and fill them with zeros
     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
@@ -24,8 +23,7 @@ int main() {
         Ar[i] = new int[M];
         A0_inv[i] = new double[M];
     }
-
+    // Fill with zeros
-    // Initialize all matrices with zeros
     for (i = 0; i < M; i++) {
         for (j = 0; j < M; j++) {
             A0[i][j] = 0;
@@ -34,11 +32,10 @@ int main() {
         }
     }
 
-    // Initialize A with M=3 and Eq. (17) from paper
+    // Initialize A with M=3 and fill acc. to 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 {
@@ -51,11 +48,13 @@ int main() {
     //         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
+    // Initialize the diagonal matrix A0,
+    // the inverse of A0_inv of diagonal matrix A0_inv
+    // and the update matrix Ar
     for (i = 0; i < M; i++) {
+        A0[i][i] = A[i][i];
         A0_inv[i][i] = 1.0/A[i][i];
         Ar_index[i] = i;
         for (j = 0; j < M; j++) {
@@ -63,21 +62,18 @@ int main() {
         }
     }
 
+    // Define pointers dim and n_updates to use in Sherman-Morrison(...) function call
     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");
     
     // Deallocate all vectors and matrices 
     for (i = 0; i < M; i++) {
-        delete [] A[i];
+        delete [] A[i], A0[i], A0_inv[i], Ar[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;
-}
+}