Merge pull request #4 from fmgjcoppens/feature/c-from-fort

Calling C++ function MaponiA3 from Fortran

.gitignore

@@ -1,4 +1,5 @@
 *.o
-SM-MaponiA3
-Sherman-Morrison
+*.mod
+cppSherman-Morrison
+fSherman-Morrison
 .vscode

Helpers.hpp

@@ -75,6 +75,23 @@ T **matMul(T **A, T **B, unsigned int size) {
     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];
+            }
+        }
+    }
+    return C;
+}
+
+
 template<typename T1, typename T2>
 T1 **outProd(T1 *vec1, T2 *vec2, unsigned int size) {
     T1 **C = new T1*[size];

Makefile

@@ -1,17 +1,44 @@
-CC=icc
-CXX=icpc
-CFLAGS=-O0 -debug full
-CXXFLAGS=-O0 -debug full -traceback
-# ARCH=-xCORE-AVX2 
+## Used compilers
+CXX = icpc
+FC = ifort
 
-DEPS = SM-MaponiA3.cpp
-OBJ = SM-MaponiA3.o main.o
+## Compiler flags
+CXXFLAGS = -O0 -debug full -traceback
+FFLAGS = -O0 -debug full -traceback
+# ARCH = -xCORE-AVX2 
 
-%.o: %.cpp $(DEPS)
-	$(CXX) $(ARCH) -c -o $@ $< $(CFLAGS)
+## Deps & objs for the C++ stuff
+cppDEPS = cppmain.cpp SM_MaponiA3.cpp SM_MaponiA3.hpp Helpers.hpp
+cppOBJ = cppmain.o SM_MaponiA3.o
 
-Sherman-Morrison: $(OBJ)
-	$(CXX) $(ARCH) -o $@ $^ $(CFLAGS)
+## Deps & objs for the Fortran stuff
+fDEPS = fmain.f90 SM_MaponiA3_mod.f90
+fOBJ = SM_MaponiA3_f.o SM_MaponiA3_mod.o fmain.o
+fLIBS = -lstdc++
+
+## Compile recipes for C++ stuff
+%.o: %.cpp $(cppDEPS)
+	$(CXX) $(ARCH) $(CXXFLAGS) -c -o $@ $<
+
+## Compile recepies for Fortran stuff
+%.o: %.f90 $(fDEPS)
+	$(FC) $(ARCH) $(FFLAGS) -c -o $@ $<
+
+## Build tagets
+.PHONY: all clean distclean
+
+all: cppSherman-Morrison fSherman-Morrison
 
 clean:
-	@rm -vf *.o
+	@rm -vf *.o *.mod
+	
+distclean: clean
+	@rm -vf cppSherman-Morrison fSherman-Morrison
+
+## Linking the C++ example program
+cppSherman-Morrison: $(cppOBJ)
+	$(CXX) $(ARCH) $(CXXFLAGS) -o $@ $^
+
+## Linking Fortran example program calling the C++ function 'Sherman_Morrison()'
+fSherman-Morrison: $(fOBJ)
+	$(FC) $(ARCH) $(FFLAGS) $(fLIBS) -o $@ $^

SM_MaponiA3.cpp

@@ -1,7 +1,7 @@
 // SM-MaponiA3.cpp
 // Algorithm 3 from P. Maponi,
 // p. 283, doi:10.1016/j.laa.2006.07.007
-#include "SM-MaponiA3.hpp"
+#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) {

SM_MaponiA3.hpp

@@ -1,3 +1,2 @@
 // SM-MaponiA3.hpp
-
-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);

SM_MaponiA3_f.cpp

@@ -0,0 +1,144 @@
+// 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

@@ -0,0 +1,4 @@
+// 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

@@ -0,0 +1,11 @@
+module Sherman_Morrison 
+    interface
+        subroutine MaponiA3(Slater0, Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3")
+            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 out) :: Slater_inv
+        end subroutine MaponiA3
+    end interface
+end module Sherman_Morrison

cppmain.cpp

@@ -1,5 +1,5 @@
 // main.cpp
-#include "SM-MaponiA3.hpp"
+#include "SM_MaponiA3.hpp"
 #include "Helpers.hpp"
 #include <cstdlib>
 #include <ctime>

fmain.f90

@@ -0,0 +1,52 @@
+program Interface_test
+    use Sherman_Morrison, only : MaponiA3
+    use, intrinsic :: iso_c_binding, only : c_int, c_double
+    implicit none
+
+    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 :: A0_inv
+
+    dim = 3
+    N_updates = dim
+    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,2) = 1
+    A(1,3) = -1
+    A(2,1) = 1
+    A(2,2) = 1
+    A(2,3) = 0
+    A(3,1) = -1
+    A(3,2) = 0
+    A(3,3) = -1
+
+    !! Prepare the diagonal matrix A0 and the update matrix Ar
+    do i=1,dim
+        Ar_index(i) = i
+        do j=1,dim
+            if (i == j) then
+                A0(i,j) = A(i,j)
+                A0_inv(i,j) = 1.0d0 / A0(i,j)
+            else
+                A0(i,j) = 0
+                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
+        do j=1,dim
+            write(*,"(F3.0,3X)", advance="no") A0_inv(i,j)
+        end do
+        write(*,*)
+    end do
+
+    deallocate(Ar_index, A, A0, Ar, A0_inv)
+end program

mwe/compile.sh

@@ -0,0 +1 @@
+icpc -c worker.cpp && ifort -c main.f90 && ifort -lstdc++ worker.o main.o -o test

mwe/main.cpp

@@ -0,0 +1,23 @@
+#include <iostream>
+
+
+int main()
+{
+
+    typedef int (*to2D)[3]; //pint2 is a pointer to an array of 2 ints
+    int linArray[9] = {0,1,2,3,4,5,6,7,8};
+    
+    to2D dArray = (to2D)linArray;
+    
+    std::cout << dArray[0][0] << std::endl;
+    std::cout << dArray[0][1] << std::endl;
+    std::cout << dArray[0][2] << std::endl;
+    std::cout << dArray[1][0] << std::endl;
+    std::cout << dArray[1][1] << std::endl;
+    std::cout << dArray[1][2] << std::endl;
+    std::cout << dArray[2][0] << std::endl;
+    std::cout << dArray[2][1] << std::endl;
+    std::cout << dArray[2][2] << std::endl;
+
+    return 0;
+}

mwe/main.f90

@@ -0,0 +1,12 @@
+program test
+  use iso_c_binding
+  implicit none
+
+  interface
+    subroutine hello() bind(C, name="worker")
+    end subroutine
+  end interface
+
+  call hello()
+
+end program test

mwe/worker.cpp

@@ -0,0 +1,7 @@
+#include "worker.h"
+
+void worker()
+{
+  std::cout << "Hello, World!" << std::endl;
+}
+

mwe/worker.h

@@ -0,0 +1,6 @@
+#include <iostream>
+
+extern "C"
+{
+  void worker();
+}