Put helper functions in separate header file from the Sherman-Morrison algo.

Helpers.hpp

@@ -0,0 +1,127 @@
+// 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];
+            }
+        }
+    }
+    showMatrix(C, size, "C");
+    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,5 +1,8 @@
 // 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;
@@ -83,7 +86,9 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
                 Al[i][j] = Id[i][j] - U[i][j] / beta;
             }
         }
+        showMatrix(Slater_inv, M, "Slater_inv");
         Slater_inv = matMul(Al, Slater_inv, M);
+        showMatrix(Slater_inv, M, "Slater_inv");
     }
 
     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);
-
-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;
-}
+void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index);

main.cpp

@@ -1,7 +1,6 @@
-// Algorithm 3 from P. Maponi,
-// p. 283, doi:10.1016/j.laa.2006.07.007
-
+// main.cpp
 #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];
     }
-
-    // Initialize all matrices with zeros
+    // Fill 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,6 +62,7 @@ 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);
@@ -71,13 +71,10 @@ int main() {
     
     // 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[i], A0[i], A0_inv[i], Ar[i];
     }
     delete [] A, A0, A0_inv, Ar, Ar_index;
     delete dim, n_updates;
 
     return 0;
-}
+}