C++ redesign of data-structures

- Use flat arrays
- Use real type for all matrices
- Merge _f MaponiA3 files

Helpers.hpp

@@ -7,8 +7,8 @@ using namespace std;
 
 template<typename T>
 unsigned int getMaxIndex(T *vector, unsigned int size) {
-    unsigned int i;     
+    unsigned int i;
-    unsigned int max = vector[0]; 
+    unsigned int max = vector[0];
     unsigned int maxi = 0;
     for (i = 1; i < size; i++) {
         if (vector[i] > max) {
@@ -16,7 +16,7 @@ unsigned int getMaxIndex(T *vector, unsigned int size) {
             maxi = i;
         }
     }
-    return maxi; 
+    return maxi;
 }
 
 template<typename T>
@@ -34,12 +34,12 @@ void showVector(T *vector, unsigned int size, string name) {
 }
 
 template<typename T>
-void showMatrix(T **matrix, unsigned int size, string name) {
+void showMatrix(T *matrix, unsigned int M, string name) {
     cout << name << " = " << endl;
-    for (unsigned int i = 0; i < size; i++) {
+    for (unsigned int i = 0; i < M; i++) {
         cout << "[ ";
-        for (unsigned int j = 0; j < size; j++) {
+        for (unsigned int j = 0; j < M; j++) {
-            cout << matrix[i][j] << " ";
+            cout << matrix[i*M+j] << " ";
         }
         cout << " ]" << endl;
     }
@@ -60,31 +60,12 @@ void showMatrixT(T **matrix, unsigned int size, string name) {
 }
 
 template<typename T>
-T **matMul(T **A, T **B, unsigned int size) {
+T *matMul(T *A, T *B, unsigned int M) {
-    T **C = new T*[size];
+    T *C = new T[M*M];
-    for (unsigned int i = 0; i < size; i++) {
+    for (unsigned int i = 0; i < M; i++) {
-        C[i] = new T[size];
+        for (unsigned int j = 0; j < M; j++) {
-    }
+            for (unsigned int k = 0; k < M; k++) {
-    for (unsigned int i = 0; i < size; i++) {
+                C[i*M+j] += A[i*M+k] * B[k*M+j];
-        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 T>
-T **matMul2(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];
             }
         }
     }
@@ -93,14 +74,11 @@ T **matMul2(T **A, T (*B)[], unsigned int size) {
 
 
 template<typename T1, typename T2>
-T1 **outProd(T1 *vec1, T2 *vec2, unsigned int size) {
+T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
-    T1 **C = new T1*[size];
+    T1 *C = new T1[M*M];
-    for (unsigned int i = 0; i < size; i++) {
+    for (unsigned int i = 0; i < M; i++) {
-        C[i] = new T1[size];
+        for (unsigned int j = 0; j < M; j++) {
-    }
+            C[i*M+j] = vec1[i+1] * vec2[j];
-    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;
@@ -140,4 +118,4 @@ T matDet(T **A, unsigned int M) {
         return det;
     }
     delete [] temp;
-}
+}

Makefile

@@ -3,9 +3,9 @@ CXX = icpc
 FC = ifort
 
 ## Compiler flags
-CXXFLAGS = -O0 -debug full -traceback
+CXXFLAGS = -O0 #-debug full -traceback
-FFLAGS = -O0 -debug full -traceback
+FFLAGS = -O0 #-debug full -traceback
-# ARCH = -xCORE-AVX2 
+# ARCH = -xCORE-AVX2
 
 ## Deps & objs for the C++ stuff
 cppDEPS = cppmain.cpp SM_MaponiA3.cpp SM_MaponiA3.hpp Helpers.hpp
@@ -13,7 +13,7 @@ cppOBJ = cppmain.o SM_MaponiA3.o
 
 ## Deps & objs for the Fortran stuff
 fDEPS = fmain.f90 SM_MaponiA3_mod.f90
-fOBJ = SM_MaponiA3_f.o SM_MaponiA3_mod.o fmain.o
+fOBJ = SM_MaponiA3.o SM_MaponiA3_mod.o fmain.o
 fLIBS = -lstdc++
 
 ## Compile recipes for C++ stuff
@@ -31,7 +31,7 @@ all: cppSherman-Morrison fSherman-Morrison
 
 clean:
 	@rm -vf *.o *.mod
-	
+
 distclean: clean
 	@rm -vf cppSherman-Morrison fSherman-Morrison
 

SM_MaponiA3.cpp

@@ -1,21 +1,20 @@
-// SM-MaponiA3.cpp
+// SM-MaponiA3_f.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) {
+void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
-    unsigned int k, l, lbar, i, j, tmp, M = *Dim;
+
+    unsigned int k, l, lbar, i, j, tmp = M;
     unsigned int *p = new unsigned int[M+1];
-    unsigned int **Id = new unsigned int*[M];
+    double *Id = new double[M*M];
     double alpha, beta;
-    double **U, *breakdown = new double[M+1];
+    double *breakdown = new double[M+1];
-    double **Al = new double*[M];
+    double *Al = new double[M*M];
     p[0] = 0;
     for (i = 0; i < M; i++) {
         p[i+1] = i + 1;
-        Id[i] = new unsigned int[M];
-        Al[i] = new double[M];
     }
 
     // Declare auxiliary solution matrix ylk
@@ -30,11 +29,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
     // Initialize identity matrix
     for (i = 0; i < M; i++) {
         for (j = 0; j < M; j++) {
-            if (i != j) Id[i][j] = 0;
+            if (i != j) Id[i*M+j] = 0;
-            else Id[i][j] = 1;
+            else Id[i*M+j] = 1;
         }
     }
-    
+
     // Initialize ylk with zeros
     for (l = 0; l < M; l++) {
         for (k = 0; k < M+1; k++) {
@@ -47,7 +46,7 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
     // 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];
+            ylk[0][k][i] = Slater_inv[(i-1)*M+(i-1)] * Updates[(i-1)*M+(k-1)];
         }
     }
 
@@ -76,19 +75,19 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
         }
     }
 
-    // Construct A-inverse from A0-inverse and the ylk
     // Keep the memory location of the passed array 'Slater_inv' before 'Slater_inv'
     // gets reassigned by 'matMul(...)' in the next line, by creating a new
     // pointer 'copy' that points to whereever 'Slater_inv' points to now.
-    double **copy = Slater_inv;
+    double *copy  = Slater_inv;
 
+    // Construct A-inverse from A0-inverse and the ylk
     for (l = 0; l < M; l++) {
         k = l+1;
-        U = outProd(ylk[l][p[k]], Id[p[k]-1], M);
+        double * U = outProd(ylk[l][p[k]], (Id + (p[k]-1)*M), 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;
+                Al[i*M+j] = Id[i*M+j] - U[i*M+j] / beta;
             }
         }
         Slater_inv = matMul(Al, Slater_inv, M);
@@ -96,16 +95,24 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, uns
 
     // Assign the new values of 'Slater_inv' to the old values in 'copy[][]'
     for (i = 0; i < M; i++) {
-        for (j = 0; j < M; j++) {
+      for (j = 0; j < M; j++) {
-            copy[i][j] = Slater_inv[i][j];
+        copy[i*M+j] = Slater_inv[i*M+j];
-        }
+      }
     }
 
     for (l = 0; l < M; l++) {
         for (k = 0; k < M+1; k++) {
             delete [] ylk[l][k];
         }
-        delete [] ylk[l], Id[l], U[l], Al[l], Slater_inv[l];
+        delete [] ylk[l];
     }
+    delete [] Id, Al;
     delete [] p, breakdown;
-}
+}
+
+extern "C" {
+  void MaponiA3_f(double **linSlater0, double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, double **linUpdates, unsigned int **Updates_index)
+  {
+    MaponiA3(*linSlater0, *linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
+}

SM_MaponiA3.hpp

@@ -1,2 +1 @@
-// SM-MaponiA3.hpp
+void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int N_updates, double *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);

SM_MaponiA3_f.cpp

@@ -1,144 +0,0 @@
-// SM-MaponiA3_f.cpp
-// Algorithm 3 from P. Maponi,
-// p. 283, doi:10.1016/j.laa.2006.07.007
-#include "SM_MaponiA3_f.hpp"
-#include "Helpers.hpp"
-
-void MaponiA3(int **linSlater0, double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, int **linUpdates, unsigned int *Updates_index) {
-
-    // Define new 2D arrays and copy the elements of the 
-    // linear passed Fortran arrays. This block needs to
-    // be replaced with a suitable casting mechanism to
-    // avoid copying of memory.
-    int **Slater0 = new int*[*Dim];
-    int **Updates = new int*[*Dim];
-    double **Slater_inv = new double*[*Dim];
-    for (int i = 0; i < *Dim; i++) {
-        Slater0[i] = new int[*Dim];
-        Updates[i] = new int[*Dim];
-        Slater_inv[i] = new double[*Dim];
-    }
-    for (unsigned int i = 0; i < *Dim; i++) {
-        for (unsigned int j = 0; j < *Dim; j++) {
-            Slater0[i][j] = linSlater0[0][i+*Dim*j];
-            Slater_inv[i][j] = linSlater_inv[0][i+*Dim*j];
-            Updates[i][j] = linUpdates[0][i+*Dim*j];
-        }
-    }
-    // Possible casting candidates
-    // int (*Slater0)[*Dim] = (int(*)[*Dim])linSlater0[0];
-    // double (*Slater_inv)[*Dim] = (double(*)[*Dim])linSlater_inv[0];
-    // int (*Updates)[*Dim] = (int(*)[*Dim])linUpdates[0];
-    ////////////////////////////////////////////////////////////////////////
-
-    unsigned int k, l, lbar, i, j, tmp, M = *Dim;
-    unsigned int *p = new unsigned int[M+1];
-    unsigned int **Id = new unsigned int*[M];
-    double alpha, beta;
-    double **U, *breakdown = new double[M+1];
-    double **Al = new double*[M];
-    p[0] = 0;
-    for (i = 0; i < M; i++) {
-        p[i+1] = i + 1;
-        Id[i] = new unsigned int[M];
-        Al[i] = new double[M];
-    }
-
-    // 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 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++) {
-            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];
-            }
-        }
-    }
-
-    // Keep the memory location of the passed array 'Slater_inv' before 'Slater_inv'
-    // gets reassigned by 'matMul(...)' in the next line, by creating a new
-    // pointer 'copy' that points to whereever 'Slater_inv' points to now.
-    // double **copy  = Slater_inv;
-
-    // Construct A-inverse from A0-inverse and the ylk
-    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);
-    }
-
-    // Overwrite the old values in 'copy' with the new ones in Slater_inv
-    // for (i = 0; i < M; i++) {
-    //     for (j = 0; j < M; j++) {
-    //         copy[i][j] = Slater_inv[i][j];
-    //     }
-    // }
-
-    // Overwrite the old values in 'linSlater_inv' with the new values in Slater_inv
-    for (i = 0; i < M; i++) {
-        for (j = 0; j < M; j++) {
-            linSlater_inv[0][i+*Dim*j] = Slater_inv[i][j];
-        }
-    }
-
-    for (l = 0; l < M; l++) {
-        for (k = 0; k < M+1; k++) {
-            delete [] ylk[l][k];
-        }
-        delete [] ylk[l], Id[l], U[l], Al[l], Slater_inv[l];
-    }
-    delete [] p, breakdown;
-}

SM_MaponiA3_f.hpp

@@ -1,4 +0,0 @@
-// SM-MaponiA3_f.hpp
-extern "C" {
-    void MaponiA3(int **linSlater0, double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, int **linUpdates, unsigned int *Updates_index);    
-}

SM_MaponiA3_mod.f90

@@ -1,11 +1,11 @@
-module Sherman_Morrison 
+module Sherman_Morrison
     interface
-        subroutine MaponiA3(Slater0, Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3")
+        subroutine MaponiA3(Slater0, Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_f")
             use, intrinsic :: iso_c_binding, only : c_int, c_double
             integer(c_int), intent(in) :: dim, n_updates
             integer(c_int), dimension(:), allocatable, intent(in) :: Updates_index
-            integer(c_int), dimension(:,:), allocatable, intent(in) :: Slater0, Updates
+            real(c_double), dimension(:,:), allocatable, intent(in) :: Slater0, Updates
             real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
         end subroutine MaponiA3
     end interface
-end module Sherman_Morrison
+end module Sherman_Morrison

cppmain.cpp

@@ -7,73 +7,52 @@
 int main() {
 
     srand((unsigned) time(0));
-    unsigned int randRange = 1; // to get random integers in range [-randRange, randRange]
     unsigned int M = 3; // Dimension of the Slater-matrix
     unsigned int i, j; // Indices for iterators
 
     // 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
+    double *A = new double[M*M]; // The matrix to be inverted
-    int **A0 = new int*[M]; // A diagonal matrix with the digonal elements of A
+    double *A0 = new double[M*M]; // A diagonal matrix with the digonal elements of A
-    int **Ar = new int*[M]; // The update matrix
+    double *Ar = new double[M*M]; // The update matrix
-    double **A0_inv = new double*[M]; // Inverse of A0
+    double *A0_inv = new double[M*M]; // The inverse
-    for (i = 0; i < M; i++) { 
+
-        A[i] = new int[M];
-        A0[i] = new int[M];
-        Ar[i] = new int[M];
-        A0_inv[i] = new double[M];
-    }
     // Fill with zeros
     for (i = 0; i < M; i++) {
         for (j = 0; j < M; j++) {
-            A0[i][j] = 0;
+            A0[i*M+j] = 0;
-            Ar[i][j] = 0;
+            Ar[i*M+j] = 0;
-            A0_inv[i][j] = 0;
+            A0_inv[i*M+j] = 0;
         }
     }
 
     // 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[0] = 1;    A[3] = 1;    A[6] = -1;
-    A[1][0] = 1;    A[1][1] = 1;    A[1][2] = 0;
+    A[1] = 1;    A[4] = 1;    A[7] = 0;
-    A[2][0] = -1;   A[2][1] = 0;    A[2][2] = -1;
+    A[2] = -1;   A[5] = 0;    A[8] = -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");
 
     // 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[i*M+i] = A[i*M+i];
-        A0_inv[i][i] = 1.0/A[i][i];
+        A0_inv[i*M+i] = 1.0/A[i*M+i];
         Ar_index[i] = i;
         for (j = 0; j < M; j++) {
-            Ar[i][j] = A[i][j] - A0[i][j];
+            Ar[i*M+j] = A[i*M+j] - A0[i*M+j];
         }
     }
 
     // Define pointers dim and n_updates to use in Sherman-Morrison(...) function call
-    unsigned int *dim = new unsigned int(M);
+    unsigned int dim = M;
-    unsigned int *n_updates = new unsigned int(M);
+    unsigned int n_updates = M;
-    Sherman_Morrison(A0, A0_inv, dim, n_updates, Ar, Ar_index);
+    MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index);
     showMatrix(A0_inv, M, "A0_inv");
-    
+
-    // Deallocate all vectors and matrices 
+    // Deallocate all vectors and matrices
-    for (i = 0; i < M; 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;
-}
+}

fmain.f90

@@ -6,7 +6,7 @@ program Interface_test
     integer i, j !! Iterators
     integer(c_int) :: dim, N_updates
     integer(c_int), dimension(:), allocatable :: Ar_index
-    integer(c_int), dimension(:,:), allocatable :: A, A0, Ar
+    real(c_double), dimension(:,:), allocatable :: A, A0, Ar
     real(c_double), dimension(:,:), allocatable :: A0_inv
 
     dim = 3
@@ -14,15 +14,15 @@ program Interface_test
     allocate(Ar_index(dim), A(dim,dim), A0(dim,dim), Ar(dim,dim), A0_inv(dim,dim))
 
     !! Initialize A with M=3 and fill acc. to Eq. (17) from paper
-    A(1,1) = 1
+    A(1,1) = 1.0d0
-    A(1,2) = 1
+    A(1,2) = 1.0d0
-    A(1,3) = -1
+    A(1,3) = -1.0d0
-    A(2,1) = 1
+    A(2,1) = 1.0d0
-    A(2,2) = 1
+    A(2,2) = 1.0d0
-    A(2,3) = 0
+    A(2,3) = 0.0d0
-    A(3,1) = -1
+    A(3,1) = -1.0d0
-    A(3,2) = 0
+    A(3,2) = 0.0d0
-    A(3,3) = -1
+    A(3,3) = -1.0d0
 
     !! Prepare the diagonal matrix A0 and the update matrix Ar
     do i=1,dim
@@ -32,13 +32,13 @@ program Interface_test
                 A0(i,j) = A(i,j)
                 A0_inv(i,j) = 1.0d0 / A0(i,j)
             else
-                A0(i,j) = 0
+                A0(i,j) = 0.0d0
                 A0_inv(i,j) = 0.0d0
             end if
             Ar(i,j) = A(i,j) - A0(i,j)
         end do
     end do
-    
+
     call MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index)
 
     do i=1,dim