C interoperability: trivial question

[rudi-gaelzer]

This is similar to the topic discussed here:

https://software.intel.com/en-us/forums/topic/387235

This should be a trivial question for you guys, but it's nagging me anyways...

Please take a look at the very simple pair of Fortran/C code below. 

#include <stdio.h>
#include <math.h>

void sub_f1( float f1, float *f2 )
{
   float f3;
   f3= sin(f1);
//    f2= &f3;
   *f2= f3;
}

program tes_sub_f1
use, intrinsic :: iso_c_binding
implicit none
interface
   subroutine sub_f1(f1, f2) bind(c)
   import :: c_float
   real(c_float), value :: f1
   real(c_float), intent(out) :: f2
   end subroutine sub_f1
end interface
real(c_float) :: ff1, ff2
!
do
   write(*, '(/,a)', advance= 'no')'Enter f1: '
   read(*,*)ff1
   print*, 'F: sin(f1)=', sin(ff1)
   call sub_f1(ff1, ff2)
   print*, 'C: sin(f1)=', ff2
end do
end program tes_sub_f1

I'm trying to understand how to pass values to a C subroutine (void function) and receive results in other variables. 

The program works all right.  I enter any real number, fortran passes the value to the sub_f1 C routine, it evaluates the sin of the number and returns the result to Fortran which displays the correct value.

However, if I uncomment the line

f2= &f3;

i. e., the pointer f2 first points to the address of f3, and after to its value, it will always return 0 to fortran.

Is this behavior expected?

Also: if I want to pass results (scalars, arrays) FROM C to Fortran, the said results must always be passed as pointers? I all the tests I did, this was the only possibility that worked.

Thanks.

[mecej4]

Line-4 of your first listing rests on a convention: that the second argument is the address into which the sine value should be placed. That address is provided by the caller, which will retrieve the value from that specific address. If, in the meantime, the C routine changes the address itself, instead of merely changing its contents, the program will probably fail and the results are garbage even if it does not fail. The situation would be somewhat similar to changing labels on a set of mail boxes during just the few minutes when the postal worker is delivering mail.

Further damage is done by the property of the variable f3 that it is a local variable, allocated on the stack. When the function returns, the stack is freed and other function calls and the caller can reuse that stack space, making it hard to debug programs with this type of bug.

You can declare f2 as a "constant address" of a variable value, using the const modifier:

float *const f2

instead of

float *f2;

In the context of C++ rather than C, there are many useful discussions regarding references versus pointers, e.g., http://www.embedded.com/electronics-blogs/programming-pointers/4023307/References-vs-Pointers .

[rudi-gaelzer]

mecej4 wrote:

Line-4 of your first listing rests on a convention: that the second argument is the address into which the sine value should be placed. That address is provided by the caller, which will retrieve the value from that specific address. If, in the meantime, the C routine changes the address itself, instead of merely changing its contents, the program will probably fail and the results are garbage even if it does not fail. The situation would be somewhat similar to changing labels on a set of mail boxes during just the few minutes when the postal worker is delivering mail.

Further damage is done by the property of the variable f3 that it is a local variable, allocated on the stack. When the function returns, the stack is freed and other function calls and the caller can reuse that stack space, making it hard to debug programs with this type of bug.

You can declare f2 as a "constant address" of a variable value, using the const modifier:

float *const f2

instead of

float *f2;

I tried your suggestion.  Leaving the f2= &f3 line commented, the program worked as before.

After uncommenting this line, I got the following error message from gcc (I don't have icc):

gcc -c sub_f2.c
sub_f2.c: In function ‘sub_f2’:
sub_f2.c:8:6: error: assignment of read-only parameter ‘f2’
    f2= &f3;
      ^

Deitel's "C++ How to Program" book says (sec. 8.8) that "(...) a constant pointer cannot be used on the left side of an assignment operator (...)."  So, I guess this is the case?

[mecej4]

Rudi-Gaelzer wrote:
..a constant pointer cannot be used on the left side of an assignment operator. So, I guess this is the case?

Exactly. The addition of the "const" modifier helps the compiler to warn you about the inconsistency in intent regarding the pointer.

[rudi-gaelzer]

mecej4 wrote:

Quote:

Exactly. The addition of the "const" modifier helps the compiler to warn you about the inconsistency in intent regarding the pointer.

Ok.  Thank you for the tip.

Now, since you were kind (and patient...) enough to reply, what I really want to do is to return a variable array from the C void.

If I have the following routine:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

void c_varr( int a1, float a2, float * const varr )
{
   int i;
   float v;
   varr= (float *) malloc(a1*sizeof(float));
   for (i= 0; i < a1; i++){
      v= (float) i + sqrt(a2);
      varr= v;
   }
}

Is this the correct interface I should use?:

interface
   subroutine c_varr(a1, a2, varr) bind(c)
   import :: c_int, c_float, c_ptr
   integer(c_int), value :: a1
   real(c_float), value  :: a2
   type(c_ptr), intent(out) :: varr
   end subroutine c_varr
end interface

Then declare the vars:

integer(c_int) :: n= 5
real(c_float) :: a= 1.5
type(c_ptr) :: farr
real(c_float), dimension(:), pointer :: farr_p

Then call the routines:

call c_varr(n, a, farr)
call c_f_pointer(farr, farr_p, )

Should this work?

Thanks again.

[mecej4 · Accepted solution]

Rudi-Gaelzer wrote:
what I really want to do is to return a variable array from the C void.

OK, that can be done, and the code to do that is given below. However, I urge you not to do allocations on the C side in a Fortran/C mixed program. Use Fortran allocatable arrays instead -- they are intended precisely for this purpose, and they are much safer to use and easier to debug than free-wheeling C allocations and type casting. Once you fall into it, it is hard to get back from the "C void"!

The C-function, to be called from Fortran:

#include <stdlib.h>
#include <math.h>

void c_varr( int a1, float a2, float **varrp )
{
   int i;
   float v,*varr;
   varr= (float *) malloc(a1*sizeof(float)); // allocate array to be handed back
   for (i= 0; i < a1; i++){
      v= (float) i + sqrt(a2);
      varr= v;                            // set array element values
      }
   *varrp=varr;                              // set handle value
}

The Fortran caller:

program calloc
use, intrinsic :: iso_c_binding
implicit none

integer(c_int) :: n= 5
real(c_float) :: a= 1.5
type(c_ptr) :: farr
real(c_float), dimension(:), pointer :: farr_p

interface
   subroutine c_varr(a1, a2, varr) bind(c)
   import :: c_int, c_float, c_ptr
   integer(c_int), value :: a1
   real(c_float), value  :: a2
   type(c_ptr), intent(out) :: varr
   end subroutine c_varr
end interface

call c_varr(n, a, farr)
call c_f_pointer(farr, farr_p, )
write(*,*)farr_p(1:n)
end

Do note that the malloc'ed array was never given back to C's free() function. In a real application, freeing C-allocated arrays in a timely manner can become a tough job.

[rudi-gaelzer]

Indeed it worked all right.

Thanks a lot for your suggestions.

BTW, when you suggested that I shouldn't do allocations on the C side, were you referring to the use of malloc inside the C function? I'm afraid there's no escaping from that in my case.  The codes above are just tests.  In my real application, the size of the array isn't predetermined, but depends on further computation done inside the C routine.

Anyway, thanks again.

[mecej4]

Fortran provides ALLOCATE/DEALLOCATE, MOVE_ALLOC and ALLOCATED() for dynamic memory management. In addition, there is automatic allocation of suitably declared variables when such a variable is set equal to an array expression.

Using the Fortran dynamic memory facilities alongside malloc/free in the same program is similar to having two finance ministers in, say, Greece.

[jimdempseyatthecove]

Array allocation in Fortran can be thought of as analogous to a C/C++ container, where the container has (unseen) information about the array that is allocated as POD outside of the container, this is called an array descriptor. The call to c_f_pointer creates an analog to the container (array descriptor), but it cannot provide all the attribute information. As an example, a c_f_pointer constructed cannot contain an indication that there is no alias to the data. Compiler optimization rules differ when it is known there cannot be an alias to the data.

The allocation on the Fortran side is preferable but not mandatory. The design of your application may not lend itself to allocation on the Fortran side (e.g. Fortran DLL/.so that is processing data that "lives" in the C++ application).

Jim Dempsey

[Steven_L_Intel1]

You might want to look at the "Further C Interoperability" features from Fortran 2015 Draft that are supported in ifort 16. This includes the ability for C code to play with Fortran allocatable arrays through the use of "C descriptors" and dedicated functions. You still can't use malloc/free, but it does give you more flexibility on the C side.

[rudi-gaelzer]

Steve Lionel (Intel) wrote:

You might want to look at the "Further C Interoperability" features from Fortran 2015 Draft that are supported in ifort 16. This includes the ability for C code to play with Fortran allocatable arrays through the use of "C descriptors" and dedicated functions. You still can't use malloc/free, but it does give you more flexibility on the C side.

I downloaded ifort 16 but still did not install it.  Is there some sample code on C descriptors in the docs, or anywhere else?  I have the WD 1539-9 document release in 12/24/2014, but the description given from section 15.4 is not too clear...

I've managed to implement routines that transfer dynamic arrays from C to F along the lines of the sample codes above.  However, I noticed some loss of accuracy during the number conversions, and so now I'm thinking on transfering the results as an array of strings into F and then convert them to floats via internal files.

[Steven_L_Intel1]

Here's an example.

use, intrinsic :: iso_c_binding
interface
  function c_alloc (array) bind(C)
  import
  integer(C_INT) :: c_alloc
  real(C_FLOAT), intent(out), allocatable, dimension(:) :: array
  end function c_alloc
end interface

real(C_FLOAT), allocatable, dimension(:) :: my_array
if (c_alloc(my_array) == 0) then
  print *, lbound(my_array), ubound(my_array); print *, my_array
  end if
end

#include "ISO_Fortran_binding.h"
extern int c_alloc (CFI_cdesc_t * descr) {
	int ret, i;
	float * array;
	CFI_index_t lower = 0, upper = 10;
	ret = CFI_allocate (descr, &lower, &upper, 0); // No elem_len
	if (ret == CFI_SUCCESS) {
		array = descr->base_addr;
		for (i=lower;i<=upper;i++) {
			array = (float) i;
		}
	}
	return ret;
}

[mecej4]

Rudi-Gaelzer wrote:
However, I noticed some loss of accuracy during the number conversions, and so now I'm thinking on transfering the results as an array of strings into F and then convert them to floats via internal files.

What conversions and why loss of accuracy? This should not happen, and if it did there is probably some error in implementation.

Writing formatted output of numbers and reading them back is both unnecessary and likely to cause more loss in accuracy. Makes no sense to me.

[rudi-gaelzer]

mecej4 wrote:

Quote:

What conversions and why loss of accuracy? This should not happen, and if it did there is probably some error in implementation.

Writing formatted output of numbers and reading them back is both unnecessary and likely to cause more loss in accuracy. Makes no sense to me.

That happened because I'm trying to do something not sanctioned by the standard here.  I want to transfer an array of extended precision (128 bits) numbers. In order to do that, I'm using the quadmath library in C and declaring the dummy arguments in the interface with the kind parameter

integer, parameter :: qp= selected_real_kind(32,4500)

Below a sample:

use, intrinsic :: iso_c_binding
integer, parameter :: qp= SELECTED_REAL_KIND(32,4500)
INTERFACE
   function sqrt_quad(x) bind(c)
   import :: qp
   real(kind= qp) :: sqrt_quad
   real(kind= qp), value :: x
   end function sqrt_quad
END INTERFACE
...
print*, sqrt_quad(2.0_dp)
...

/* function returning square root (quad) */
#include <quadmath.h>
#include <stdio.h>

__float128 sqrt_quad(__float128 x)
{
   __float128 res;
   res= sqrtq(x);
   return res;
}

It works, but I noticed that in some cases there's loss of accuracy.  I guess that happens because real*16 and c-quad numbers are not interoperable kinds, according to the standard.

That's why I'm thinking now on transfering the ouput from the C function as a string and then convert it to real*16 using internal files.

[rudi-gaelzer]

Ok, now I'm somewhat bogged down again...

Trying to implement the routines that transfer a variable array of strings from C to F, I came up with the following samples.

This is the sample C function.

void c_sub (long leng, long *nelm, char **castr)
{
   int i;
   long k;
   char *temp;
   ...
// Generate a result comprised by k strings of length=leng each.
// Each string will be obtained inside a loop.
   ... 
   castr= (char **) malloc(k*leng);
   for (i= 0; i < k; i++)
   {
      ...
      temp= <string obtained>;
      castr= temp;
   }
   *nelm= k;
}

Now the fortran program.

program call_c_sub
use, intrinsic :: iso_c_binding
implicit none
integer, parameter :: mls= 50
integer :: lenstr
integer(c_long) :: ncf, i, leng
type(c_ptr) :: fastr
character(kind=c_char), dimension(:,:), pointer :: fastr_p
character(len= mls) :: str
interface
   subroutine c_sub(leng, nelm, castr) bind(c)
   import :: c_long, c_ptr
   integer(c_long), value       :: leng
   integer(c_long), intent(out) :: nelm
   type(c_ptr), intent(out)     :: castr
   end subroutine c_sub
end interface
!
leng= mls
call c_sub(leng, ncf, fastr)
call c_f_pointer(fastr, fastr_p, [leng, ncf])
!
print*, 'F: # coeff= ', ncf
do i= 1, ncf
   lenstr= cstrlen(fastr(:,i)) !String length up to c_null_char
   str= transfer(fastr(1:lenstr,i), str)
   print*, 'i=', i, '  string=', trim(str)
end do
CONTAINS
   function cstrlen(cstr) result (res) ! Locates c_null_char
   integer :: res
   character(kind=c_char), dimension(mls), intent(in) :: cstr
   integer :: ii
   do ii= mls, 1, -1
      if(cstr(ii) == c_null_char)then
         res= ii - 1
         exit
      end if
   end do
   return
   end function cstrlen
end program call_c_sub

The program compiles, and runs up to the line 

   lenstr= cstrlen(fastr(:,i)) !String length up to c_null_char

when a segmentation fault occurs.

Obviously, I'm doing something wrong with the handling of the string array coming from C.

Could you "point" me to the right direction here?

Thanks.

[mecej4]

Sorry, I don't think that I can help you. I have to admit my misgiving that what you are attempting to do is probably going to fail for a number of reasons.

You seem to be attempting to work in the shadows attempting to get C and Fortran to interoperate with non-interoperable types. There is next to no portability across compilers w.r.t. quadruple and higher precision. Implementing your own high precision arithmetic with support for trancendental functions is no trivial task.

In fact, quad-precision is, rather vaguely, something with more than double precision, and there are multiple implementations of quad-precision. I'd need to see strong justification of the benefits of higher precision before investing time on getting it to work.

[rudi-gaelzer]

mecej4 wrote:

Sorry, I don't think that I can help you. I have to admit my misgiving that what you are attempting to do is probably going to fail for a number of reasons.

You seem to be attempting to work in the shadows attempting to get C and Fortran to interoperate with non-interoperable types. There is next to no portability across compilers w.r.t. quadruple and higher precision. Implementing your own high precision arithmetic with support for trancendental functions is no trivial task.

In fact, quad-precision is, rather vaguely, something with more than double precision, and there are multiple implementations of quad-precision. I'd need to see strong justification of the benefits of higher precision before investing time on getting it to work.

Indeed.  That's why I'm trying to pass the results as strings and only convert to floats at the final destination (i. e., Fortran)...

Thanks a lot for your help.