2
$\begingroup$

I am writing a Fortran function that needs to receive two vectors of reals, and as an output, returns a vector. The function could be for example the sum of the two vectors:

function sumtwovectors( n, v1, v2 )

implicit none
real(kind=8), dimension(100000)             :: sumtwovectors
integer, intent(in)                         :: n
real(kind=8), dimension(100000), intent(in) :: v1, v2

integer :: i

if (n .gt. 100000) then
    sumtwovectors(:) = -99
    return
endif

do i = 1, n
    sumtwovectors(i) = v1(i) + v2(i)
enddo

return

end function sumtwovectors

which is then compiled without error into a dll with

gfortran -c sumtwovectors.f90 -o sumtwovectors.o
gfortran -shared -mrtd -o test.dll sumtwovectors.o

This dll should be NET-compatible, so I load it into Mathematica ( path contains the path where the dll is stored):

InstallNET[];
sumtwovectors = DefineDLLFunction["sumtovectors_", path, "double[]", {"short*", "double[]", "double[]"}]

but when I run this, I have an error:

sumtwovectors[3, {1, 3, 2}, {3, 4, 5}]

NET::netexcptn: A .NET exception occurred: System.Runtime.InteropServices.MarshalDirectiveException: Impossible de marshaler 'return value' : Combinaison de types managés/non managés non valide.
à Wolfram.NETLink.DynamicDLLNamespace.DLLWrapper1.sumtovectors_(Int16& , Double[] , Double[] ).
  1. How can I make this work in Mathematica? I tried enclosing the two list into MakeNETObject[] but to no avail.
  2. Instead of having a maximum size (1000000), could it work with no defined maximum size?
$\endgroup$
3
$\begingroup$

I'm not familiar with .NET interface for WM, but your problem can be solved with LibraryLink. With LibraryLink it is possible to create a shared lib which can be loaded with LibraryFunctionLoad[] or, perhaps, this lib can be accessed with .NET interface.

Fortran has interoperability with C. While it is possible to call Fortran functions with <name>_, the modern way is to use the intrinsic iso_c_binding module.

You want to interface a function that returns an array. It's OK in Fortran, but C functions do not return arrays. A solution is to use a subroutine in Fortran and a void function in C.

Create funs.f90, test.c, test.f90 and interface.cpp in your notebook directory (see files contents below or download all files forum.tar.gz)

(* set directory *)
SetDirectory[NotebookDirectory[]] ;

You can define your Fortran procedures in a module.

(* fortran module *)
FilePrint["funs.f90"]
module funs
  use, intrinsic :: iso_c_binding, only : ik => c_int, rk => c_double
  implicit none
  private
  public :: ik
  public :: rk
  public :: vsum
  contains
  ! void vsum(int, double*, double*, double*) ;
  pure subroutine vsum(n, v1, v2, v3) bind(c, name = "vsum")
    integer(ik), value, intent(in) :: n
    real(rk), dimension(n), intent(in) :: v1
    real(rk), dimension(n), intent(in) :: v2
    real(rk), dimension(n), intent(out) :: v3
    v3 = v1 + v2
  end subroutine vsum
end module funs

Note, array dims are set by the input arg n. Array length n also has value keyword which allows to pass it by value in C.

(* compile fortran module *)
Run["gfortran -c -std=f2008 -O3 funs.f90"]
(* 0 *)

Here is a test Fortran program that uses this module.

(* fortran test program *)
FilePrint["test.f90"]
program test
  use funs
  implicit none
  integer(ik), parameter :: n = 16_ik
  real(rk), dimension(n) :: v1, v2, v3
  integer(ik) :: i
  v1 = real([(i, i = 1_ik, n, 1_ik)], rk)
  v2 = v1
  call vsum(n, v1, v2, v3)
  do i = 1_ik, n, 1_ik
    write(*,*) v3(i)
  end do
end program test
(* compile and run fortran test program *)
Run["gfortran -c -std=f2008 -O3 test.f90"]
Run["gfortran -o test -std=f2008 -O3 funs.o test.o"]
ImportString[RunProcess["./test", "StandardOutput"],"List"]
(*0*)
(*0*)
(*{2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,26.,28.,30.,32.}*)

Here is a test C program that uses this module.

(* c test program *)
FilePrint["test.c"]
#include <stdio.h>
void vsum(int, double*, double*, double*) ;
int main(){
  int n = 16 ;
  double x[n], y[n], z[n] ;
  for(int i = 0 ; i < n ; i++){
    x[i] = (double) i + 1 ;
    y[i] = (double) i + 1;
    z[i] = 0.0 ;
  }
  vsum(n, x, y, z) ;
  for(int i = 0 ; i < n ; i++){
    printf("%lf\n",z[i]) ;
  }
  return 0 ;
}
(* compile and run c test program *)
Run["gcc -c -std=c99 -O3 test.c"]
Run["gcc -o test -std=c99 -O3 funs.o test.o"]
ImportString[RunProcess["./test", "StandardOutput"],"List"]
(*0*)
(*0*)
(*{2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,26.,28.,30.,32.}*)

For LibraryLink you need to define a wrapper.

(* library link *)
FilePrint["interface.cpp"]
#include "WolframLibrary.h"
#include "WolframCompileLibrary.h"
DLLEXPORT mint WolframLibrary_getVersion(){ return WolframLibraryVersion ; }
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){ return 0 ; }
extern "C" { void vsum(int, double*, double*, double*) ;}
EXTERN_C DLLEXPORT int vsum_ll(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
  mint n = MArgument_getInteger(Args[0]) ;
  MTensor v1 = MArgument_getMTensor(Args[1]) ;
  MTensor v2 = MArgument_getMTensor(Args[2]) ;
  MTensor v3 = 0 ;
  libData->MTensor_clone(v1, &v3);
  vsum(
    n,
    MTensor_getRealDataMacro(v1),
    MTensor_getRealDataMacro(v2),
    MTensor_getRealDataMacro(v3)
  ) ;
  MArgument_setMTensor(Res,v3) ;
  return LIBRARY_NO_ERROR ;
}
(* create library *)
Needs["CCompilerDriver`"] ;
CreateLibrary[
    {"interface.cpp", "funs.o"},
    "vsum_ll",
    "TargetDirectory" -> NotebookDirectory[],
    "CleanIntermediate" -> True,
    "Debug" -> False,
    "CompileOptions" -> "-O3"
] ;

This creates vsum_ll lib with vsum_ll function that calls vsum procedure. You can try to load it with .NET or with LibraryFunctionLoad[].

(* load library *)
Quiet[LibraryFunctionUnload[lib]] ;
ClearAll[lib] ;
lib = LibraryFunctionLoad[
    FileNameJoin[{NotebookDirectory[],"vsum_ll"}],
    "vsum_ll",
    {Integer,{Real,1},{Real,1}},
    {Real,1}
] ;

This function can now be used in WM.

(* test *)
n = 16 ;
v1 = v2 = N[Range[n]] ;
lib[n,v1,v2]
(*{2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,26.,28.,30.,32.}*)

As an extra step, you can wrap it into Compile[], this allows to use listability and also provides some safety.

(* wrapper *)
ClearAll[vsum] ;
vsum = Compile[
    {{n,_Integer},{v1,_Real,1},{v2,_Real,1}},
    lib[n,v1,v2],
    CompilationTarget -> "C",
    CompilationOptions -> {"InlineExternalDefinitions" -> True, "InlineCompiledFunctions" -> True},
    RuntimeOptions -> "Speed",
    Parallelization -> True,
    RuntimeAttributes -> {Listable}
] ;

Test wrapped version and instability.

(* test *)
(* sequential *)
vsum[n,v1,v2]
(* listable *)
vsum[n,{v1,v1},{v2,v2}]
(*{2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,26.,28.,30.,32.}*)
(*{{2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,26.,28.,30.,32.},{2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,26.,28.,30.,32.}}*)

Unload lib when done.

(* unload *)
LibraryFunctionUnload[lib]

You can find more examples of Fortran in WM here or here.

$\endgroup$
4
  • $\begingroup$ Thank you very much for you very complete answer. It was very helpful. However, the complexity of the process (using a C wrapper to a fortran function) is precisely something I wanted to avoid. $\endgroup$ – Denis Cousineau Dec 7 '20 at 2:46
  • $\begingroup$ "... but C functions do not return arrays" Ahah, I didn't know about this. Modifying the fortran code to have two input vectors and one ouptut vectors should make things easier? $\endgroup$ – Denis Cousineau Dec 7 '20 at 2:47
  • $\begingroup$ @DenisCousineau, yes, Fortran code is changed (function is replaced by subroutine) to match void C function. The wrapper is for Cpp which is used in LibraryLink, it uses an external C function which in fact is a Fortran sub. Basically, you just need to interface an external C in C#. $\endgroup$ – I.M. Dec 7 '20 at 6:02
  • $\begingroup$ @DenisCousineau, if you insist on using function in Fortran, you can wrap it in subroutine and then use it in C. Perhaps, this can be done with pointers, but I do not have such experience. $\endgroup$ – I.M. Dec 7 '20 at 6:04
2
$\begingroup$

The comment from I.M. "C functions do not return arrays" made me try a different approach with a subroutine rather than a function.

Let the Fortran subroutine be

subroutine vsum( n, v1, v2, v3 )

implicit none
integer, intent(in)                          :: n
real(kind=8), dimension(100000), intent(in)  :: v1, v2
real(kind=8), dimension(100000), intent(out) :: v3

integer :: i

do i = 1, n
    v3(i) = v1(i) + v2(i)
enddo

return

end subroutine vsum

in which v3 is an ouptut argument that will be used to return the answer. This file is compiled to a dynamically-link library (DLL) with

gfortran -c vsum.f90 -o vsum.o
gfortran -shared - mrtd -o vsum.dll vsum.o

In Mathematica, load and install the NETLink library, define the function, and initialize a container NETres that will receive the result:

Needs["NETLink`"]
InstallNET[];
sumtwovectors = 
  DefineDLLFunction["vsum_", path <> "vsum.dll", 
   "void", {"short*", "double[]", "double[]", "double[]"}
];
NETres = MakeNETObject[{0., 0., 0.}];

You are ready to use the function and unpack the NETres container:

sumtwovectors[3, {N[Pi], 10., 1000.}, {100., 100., 100.}, NETres]
NETObjectToExpression[NETres]
$\endgroup$
1
  • $\begingroup$ This approach was also mentioned 8 years ago here, though I am not so sure how well supported and robust NETLink is (on Windows, but also on Linux and MacOSX, as Mono). $\endgroup$ – Rolf Mertig Dec 10 '20 at 0:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.