Call an Intel Fortran DLL from MATLAB

[eos_pengwern]

I need to call aDLL written in IVF V10.1 from MATLAB. In principle everything works well, except I can't seem to pass real arrays back and forth.

MATLAB's arrays are pretty much like Fortran allocatable arrays; indices start at 1, the arrays can be multi-dimensional, and are stored in column-major order. However, MATLAB expects all external DLLs to be written in C, so converts all data to C-format on the way out and on the way back. My Fortran DLL, in turn, uses ISO_C_BINDING and declares all data types as C-compatible.

C of course doesn't support multidimensional allocatable arrays, but this shouldn't in theory be a problem. All that should be passing back and forth between MATLAB and Fortran are pointers to the first element of each array, and the arrays themselves should never actually need to be converted. Naturally the parameters that I pass from MATLAB to Fortran include the array dimensions, so that Fortran can allocate the arrays correctly.

In any case, what actually happens is that integer values are passed correctly, so I can allocate arrays with the correct dimensions. However, the real values within the matrices themselved get completely corrupted. The values that Fortran receives from MATLAB are all wrong (typically by 199 orders of magnitude), and the values that MATLAB received back from Fortran are similarly wrong (though not quite so dramatically).

As well as using ISO_C_BINDING, I have tried out the compiler attributes C and STDCALL, but neither seem to affect the behavior in any way. So, what am I missing?

[Steven_L_Intel1]

I have never used MATLAB, but some of the comments you make stand out to me. Are you sure that MATLAB addresses arrays in column-major order? That is, as you say, natural for Fortran but not for C. The other thing I note is that in C, a multidimensional array is actually an array of array pointers.

What I do when faced with a problem such as this is, on the Fortrtan side, write out the hex value of LOC(arg) and then write hex values for elements of the array. This usually gives me a clue as to what I'm dealing with - whether it be elements in the wrong order, or pointers instead of data, etc. You may also have a datatype mismatch, for example, single precision vs. double. Your comments about the results being off strongly suggest that.

[eos_pengwern]

Thank you for your response.

Please forgive my ignorance, but what is the best way to output the hex values? So far I have been troubleshooting this by popping up a console window from within my DLL to print out the array values. If I try doing this in hex (using the 'Z' format specifier), I get some strange repeated patterns (e.g. trying to print out an array element in 'Z30' format gives me 12CE017812CE0178). Also, I can't find a format which willhandle the pointer value returned by c_log(array_name) without giving me a compile-time error.

[Steven_L_Intel1]

You would use a Z format with a width of 2X the number of bytes in the data type. So that would be Z8.8 for REAL(4) and Z16.16 for REAL(8). That's for one element. The 8.8 provides for leading zeroes.

The value you show looks a lot like an address to me.

You can't use C_LOC for this purpose - use LOC instead.

[eos_pengwern]

Again, thank you.

OK, the location of a one-dimensionalarray passed from MATLAB was 12EB27D0, or (reported as an integer) 317401040; the hex values of the first two elements were:

12EB017812EB0178, 69706D662062696C

For a two-dimension array passed at the same time, the location was 12EB27E8, or (reported as an integer) 317401064, and the hex values of the first two elements were:

12EB017812EB0178 (yes, the same as before), 6C2E6C736D693D53

Just for fun, I repeated all this after forcing MATLAB to pass the values in single-precision rather thandouble. Thelocations were then 12CD27D0 (315434960) and 12CD27E8 (315434984), the first elements of each array became 12CD017812CD0178, and the second elements were exactly the same as before...

All of which I'm sure tells us something, but I'm not sure what. Any further advice would be appreciated.

[Steven_L_Intel1]

Here's what I interpret from what you have shown.

12EB27D0 is the address passed to you by MATLAB. This points to an array of two pointers with values of 12EB0178 that presumably point to data. The next two "elements" are some random data that may or may not be part of what you want to see.

It is not clear to me why you get passed two pointers, but as I said, I don't know MATLAB. In any event, what MATLAB is passing to you is not an array of reals in the Fortran sense but rather a pointer to an array of reals (presumably).

Since you're using ISO_C_BINDING, try declaring the argument received as type C_PTR and then use C_F_POINTER to "convert" that to a Fortran pointer to an array of the appropriate datatype and shape.

[eos_pengwern]

Thank you once again; so, if I understand correctly, my routine should start to look like...

subroutine DLLTest1(output_ptr, ncols, nrows, &
                    one_dimension_input_ptr, &
                    two_dimension_input_ptr) &
                    bind(c, name='dlltest1')    !DEC$ ATTRIBUTES DLLEXPORT :: DLLTest1    use iso_c_binding    implicit none    integer(c_int), value :: ncols, nrows    type(c_ptr) :: one_dimension_input_ptr, two_dimension_input_ptr     type(c_ptr), intent(out) :: output_ptr    real(c_double), dimension(:), pointer :: one_dimension_array    real(c_double), dimension(:,:), pointer :: two_dimension_array    real(c_double), allocatable, dimension(:,:) :: outputarray    integer i,j    call c_f_pointer(one_dimension_input_ptr, one_dimension_array, (/ nrows /))    call c_f_pointer(two_dimension_input_ptr, two_dimension_array, (/ nrows, ncols /))    allocate (outputarray(nrows,ncols))    ! For now we just set the matrix elements to arbitrary values    forall (i=1:nrows, j=1:ncols) outputarray(i,j)=real(i+j, c_double)
    output_ptr=c_loc(outputarray)end subroutine DllTest1
...which compiles and runs fine so long as I don't try to access the values of the input arrays within the code.
However, MATLAB still receives the wrong answers back (which may be something I need to take up with them), and I can't tell whether Fortran is getting the right matrices because any attempt to read an element leads to an exception.
Am I at least on the right track?

[Steven_L_Intel1]

If you're getting access violations, then we still don't understand what's being passed.

I'll comment that your output array gets deallocated on exit from the routine. You want that to be POINTER instead. I assume MATLAB is expecting the output array to be allocated by the routine?

You need to do some more analysis of what MATLAB is passing. Let's start with this:

subroutine DLLTest1(output_ptr, ncols, nrows, &
one_dimension_input_ptr, &
two_dimension_input_ptr) &
bind(c, name='dlltest1')

!DEC$ ATTRIBUTES DLLEXPORT :: DLLTest1

integer output_ptr, ncols, nrows, one_dimension_input_ptr, two_dimension_input_ptr)

! First, let's look at the argument list

print 101, 'Arguments', loc(output_ptr), loc(ncols), loc(nrows), &
loc(one_dimension_input_ptr), loc(two_dimension_input_ptr)
101 format (A,5Z9.8)

! Assume that we're being passed good addresses - let's see what's there

print 101, 'Values', output_ptr, one_dimension_input_ptr, two_dimension_input_ptr
end subroutine DLLtest1

[eos_pengwern]

Thank you; the output was:

Arguments: 127F7420 00000002 000000040300DA20 127FAC20

Values 00000000 00000000 00000000

N.B. '2 and '4' were the integer values of ncols and nrows that I passed.

[Steven_L_Intel1]

That's interesting - and doesn't match what you provided earlier. Perhaps I don't understand how you wrote the code to display the array values.

Now let's try it this way.

subroutine DLLTest1(output_ptr, ncols, nrows, &
one_dimension_input_ptr, &
two_dimension_input_ptr) &
bind(c, name='dlltest1')

!DEC$ ATTRIBUTES DLLEXPORT :: DLLTest1

integer(8), dimension(8) output_ptr, one_dimension_input_ptr, two_dimension_input_ptr
integer ncols,nrows

101 format (A,8Z16.17)

print 101, '1', one_dimension_input_ptr
print 101, '2', two_dimension_input_ptr
print 101, '3', output_ptr
end subroutine DLLtest1

[jimdempseyatthecove]

A couple of possibilities come to mind.

You have

real(c_double), allocatable, dimension(:,:) :: outputarray

The specification says:

If an allocatable array does not have the SAVE attribute, it has the status of "not currently allocated" at the beginning of each invocation of the procedure. If the array's status changes to "currently allocated", it is deallocated if the procedure is terminated by execution of a RETURN or END statement.

That is to sayyour array is deallocated at the end of the subroutine (not what you wanted)

To make it persistent then either add SAVE to the attribues

real(c_double), allocatable, save, dimension(:,:) :: outputarray

Or make it a pointer

real(c_double), pointer, dimension(:,:) :: outputarray

You will most likely want it as a pointer if you are allocating multiple persistent arrays.

Eventually, at some point in time you will need to return the allocated memory, the C routine is not maintaining the Fortran descriptor for the array. therefore the deallocation becomes problematic. The suggestion I have for you is for each collection of allocations that you create a Fortran derived type containing the pointers to the array(s) plus a link pointer to a self type.

type AllocationNode
type(AllocationNode), pointer :: next
real(c_double), pointer, dimension(:,:) :: outputarray
end type AllocationNode

Then in a module data area create a head of list pointer

type(AllocationNode), pointer :: AllocationNodeHead => NULL()

Then maintain the list of allocations. i.e. first allocate anAllocationNode and link to list as head. Then allocate using descriptor held by pointer.When time comes to deallocate, the C routine will pass back the location of the 1st element of outputarray. Use the address of the first element provide by C as a search key in the list of nodes testing the address of the 1st element of the arrays held by the nodes in the list. Perform the deallocation of the memory using the descriptor held in the node Then unlink and return the node.

Note, earlier verstions of IVF had problems using allocatable within the AllocationNode, that is why pointer is used. Your preference may be to use allocatable in place of pointer.

And just for kicks try adding a diagnostic test

output_ptr=c_loc(outputarray)
if(output_ptr .ne. loc%(outputarray(1,1)) stop

Jim Dempsey

[eos_pengwern]

If the results were different before, maybe it was because I'm no longer using the C-format variable types and iso_c_binding.

Anyway, the output this time was:

1 3FF0000000000000 3FF00000000000000 EADB2B2500000000 FFFFFFE41270AD9F F14C8F2800000000 000000111270315F 0000000000000000 FFFFFFFF00000000

2 4000000000000000 40000000000000000 etc. (i.e. all eight numbers the same)

3 4008000000000000 40080000000000000 etc. (i.e. all eight numbers the same)

[eos_pengwern]

In reply to Jim's mail, in this case I actually don't need the array values to be persistent and it suits me very well for the array to be automatically deallocated when the Fortran routine returns; each call from MATLAB to Fortran will be entirely self-contained.

I can see, however, that there will be a problem if the Fortran routine is just passing back to Matlab the C address of a Fortran array which, by the time MATLAB tries to read it, will have ceased to exist.

This might explain why I'm having difficulty sending results back to MATLAB, but not why I'm having difficulty receiving the array data in the first place.

Thanks, in any case, and I'll try making the output array a SAVEd array when I've solved the other issue.

[Steven_L_Intel1]

Ok. Now it looks as if MATLAB is just passing you simple arrays of double-precision numbers and that the output array is already allocated, so it's not correct for you to try reallocating it.

Try the same program except change the INTEGER(8) to REAL(8) and replace the Z17.16 with E15.6E3 . You should see some ordinary-looking floating point numbers that will tell you in what order things are passed.

You should be able to use REAL(C_DOUBLE) for these with ISO_C_BINDING.

[eos_pengwern]

OK, so that worked perfectly, both with and without the ISO_C_BINDING, although I haven't tried passing any values back to MATLAB yet.

For the actual application, however, the incoming arrays may be any size at all so the arrays will need to be allocatable, not fixed-size as in this example. How do I get back to that? What exactly was wrong with the way I started out...

[Steven_L_Intel1]

You don't want them allocatable - you want adjustable sizes. Since the ncols and nrows are being passed in, use them, for example:

REAL(C_DOUBLE) ARRAY1(NROWS,NCOLS)

An alternative is to use (*) for the bounds, but that works for a single dimension only.

For your output array you need to determine if MATLAB has preallocated it or if you need to do so. If you need to allocate it, then I suggest creating a local POINTER array, allocating it, then assigning C_LOC of the array to the C_PTR argument. Given the test case you have so far, I'd say the output array is already allocated and all you have to do is assign into it.

[eos_pengwern]

Yes, that works a treat! The arrays are correctly passed both in and out of the Fortran subroutine.

So it seems that at the end of the day, apart from the need to use the C-consistent data types and the C binding for the subroutine name (and of course to supply MATLAB with a C-style .h file), writing a routine like this to be called from MATLAB is exactly the same as writing a routine to be called from Fortran... I'm afraid I was sucked too far into the topic of converting C pointer arrays (on the assumption that the array structure information held by MATLAB would be lost in going through the C-style interface), and that's obviously where I went wrong!

Thank you for your patience in resolving this.

[kevin_k]

Hi Steve,

Would you mind sharing a working DLL code that's callable from Matlab?

Thanks a bunch!

Kevin

[Steven_L_Intel1]

If you're asking me, I don't have MATLAB so can't help with this.

[eos_pengwern]

Kevin,

Here's an example of an 'interface module' I wrote to transfer data to and fro between MATLAB and a completely ordinary single-precision Fortran subroutine (which isn't included in the sample):

module c_interfaces
use iso_c_binding
implicit none
interface
subroutine S_pol_single_prec(nthetas, nlambdas, reflectance, &
  thetas, wavelengths, nparms, parmvector) 
integer, value :: nthetas, nlambdas, nparms
real(kind(1e0)), intent(in), dimension(nthetas) :: thetas 
real(kind(1e0)), intent(in), dimension(nlambdas) :: wavelengths
real(kind(1e0)), intent(in), dimension(nparms) :: parmvector
real(kind(1e0)), intent(inout), dimension(nthetas, nlambdas) :: reflectance
end subroutine S_pol_single_prec 
end interface
contains
subroutine calc_s_reflectance_single(reflectance, nthetas, nlambdas, &
 thetas, wavelengths, intensities, &
 nparms, parmvector) &
 bind(c, name='calc_s_reflectance_single')
!DEC$ ATTRIBUTES DLLEXPORT :: calc_s_reflectance_single
! The prototype to go into the C header file will be:
!
! void calc_s_reflectance_single(double *, int, int, double *, double *, 
!  double *, int, double *);
!
 
! MATLAB will parse the arguments list and treat any arrays passed by reference 
! as function return values. Hence the MATLAB version of the routine (as seen
! by using 'libfunctions' with the '-full' option after the library name) 
! will be
!
! [doublePtr, doublePtr, doublePtr, doublePtr] = calc_s_reflectance_single ...
! (doublePtr, int32, int32, doublePtr, doublePtr, int32, doublePtr)
!
! ...and of course as usual in MATLAB the first returned value can be taken and
! the rest ignored.
! These are the C-compatible input and output variables
integer(c_int), value :: nthetas, nlambdas, nparms
real(c_double), intent(in), dimension(nthetas) :: thetas 
real(c_double), intent(in), dimension(nlambdas) :: wavelengths
real(c_double), intent(in), dimension(nthetas,nlambdas) :: intensities
real(c_double), intent(in), dimension(nparms) :: parmvector
real(c_double), intent(inout), dimension(nthetas, nlambdas) :: reflectance
! These are their native Fortran counterparts
integer :: nthetas_F, nlambdas_F, nparms_F
real(kind(1e0)), dimension(nthetas) :: thetas_F 
real(kind(1e0)), dimension(nlambdas) :: wavelengths_F
real(kind(1e0)), dimension(nparms) :: parmvector_F
real(kind(1e0)), dimension(nthetas, nlambdas) :: reflectance_F 
! Convert the input variables to native Fortran
nthetas_F=int(nthetas)
nlambdas_F=int(nlambdas)
nparms_F=int(nparms)
thetas_F=real(thetas, kind(1e0))
wavelengths_F=real(wavelengths, kind(1e0))
parmvector_F=real(parmvector, kind(1e0))
reflectance_F=real(reflectance, kind(1e0))
! Call the Fortran subroutine
call S_pol_single_prec(nthetas_F, nlambdas_F, reflectance_F, &
thetas_F, wavelengths_F, nparms_F, parmvector_F)
! Convert the output var
iable to C format
reflectance=real(reflectance_F, c_double)
end subroutine calc_s_reflectance_single
end module c_interfaces
This should be self-explanatory, and I find it works well.
Stephen.

[kevin_k]

Stephen,

Thanks a lot for sharing your code. It will save me a lot of time.

Best,

Kevin