The Fortran interface to C++ fuction MaponiA3() works but the mechanism of passing the 2D arrays from Fortran to C++ must be improved. Now the passed 'linear' 2D arrays are copied and reshaped into usable 2D arrays. This needs to be replaced with some suitable casting mechanism.

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

@@ -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.o SM_MaponiA3_mod.o fmain.o
+fOBJ = SM_MaponiA3_f.o SM_MaponiA3_mod.o fmain.o
 fLIBS = -lstdc++
 
 ## Compile recipes for C++ stuff

SM_MaponiA3.hpp

@@ -1,5 +1,2 @@
 // SM-MaponiA3.hpp
-extern "C" {
-    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,145 @@
+// 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 some casting mechanism to avoid
+    // copying of arrays.
+    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];
+            }
+        }
+    }
+
+    // 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;
+
+    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);
+    }
+
+    // // 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++) {
+    //         copy[i][j] = Slater_inv[i][j];
+    //     }
+    // }
+
+    // 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++) {
+            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

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

fmain.f90

@@ -1,6 +1,6 @@
 program Interface_test
+    use Sherman_Morrison, only : MaponiA3
     use, intrinsic :: iso_c_binding, only : c_int, c_double
-    use MYMODULE, only : MYSUBROUTINE
     implicit none
 
     integer i, j !! Iterators
@@ -39,14 +39,14 @@ program Interface_test
         end do
     end do
     
-    ! do i=1,dim
-    !     do j=1,dim
-    !         write(*,"(I)", advance="no") Ar_index(i)
-    !     end do
-    !     write(*,*)
-    ! end do
+    call MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index)
 
-    call MYSUBROUTINE(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/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;
+}