Calling FORTRAN function from R issue?

Section 5.2 of the R manual (on Extending R) says that only
FORTRAN subroutines can be called (not functions), probably
because of the non-standard way the compilers map FORTRAN
function names to symbols in the DLL.

This is consistent with the interface prototype for the BLAS
routine zdotc contained in <R>/include/R_ext/BLAS.h, namely,

BLAS_extern Rcomplex
    F77_NAME(zdotc)(Rcomplex * ret_val, int *n,
		    Rcomplex *zx, int *incx, Rcomplex *zy, int *incy);

[...]

On the other hand, this is not consistent with the standard
FORTRAN definition for zdotc that is contained in
<R>/src/extra/blas/cmplxblas.f, where the first argument is
n, not ret_val. 

Consequently, it is not clear where the wrapper
is defined that is called via the prototype.

Hello,

I am trying to call the BLAS Level1 function zdotc from R via
a .C call like this:

#include "R.h"
#include "R_ext/BLAS.h"

void testzdotc() {
   Rcomplex zx[3], zy[3], ret_val;

   zx[0].r = 1.0; zx[0].i = 0.0;
   zx[1].r = 2.0; zx[0].i = 0.0;
   zx[2].r = 3.0; zx[0].i = 0.0;

   zy[0].r = 1.0; zy[0].i = 0.0;
   zy[1].r = 2.0; zy[0].i = 0.0;
   zy[2].r = 3.0; zy[0].i = 0.0;

   int n=3, incx=1, incy=1;
   F77_CALL(zdotc)(&ret_val, &n, zx, &incx, zy, &incy);
   Rprintf("ret_val = %f, %f\n", ret_val.r, ret_val.i);
}

This does not work. When I run '.C('testzdotc')' there is
typically a delay for a second or so, then I get: 0.0, 0.0
instead of the correct ans: 14.0, 0.0.

On 06-03-2012, at 01:21, Dominick Samperi wrote:

   zx[0].r = 1.0; zx[0].i = 0.0;
   zx[1].r = 2.0; zx[0].i = 0.0;
   zx[2].r = 3.0; zx[0].i = 0.0;

I tried calling zdotc  through an intermediate Fortran routine hoping
it would solve your problem.

Above C routine changed to

The fortran subroutine is

<code>
      subroutine callzdotc(retval,n, zx, incx, zy, incy)
      integer n, incx, incy
      double complex retval, zx(*), zy(*)
      external double complex zdotc

      retval = zdotc(n, zx, incx, zy, incy)

      return
      end
</code>

Made a shared object with

R CMD SHLIB --output=dozdot.so callzdotc.f czdot.c 

and ran

dyn.load("dozdot.so")
.C("testzdotc")

with the result 0.0, 0.0

G'day Berend,

[..]

I tried calling zdotc  through an intermediate Fortran routine hoping
it would solve your problem.

[...] 

Above C routine changed to

[...]

The fortran subroutine is

<code>
     subroutine callzdotc(retval,n, zx, incx, zy, incy)
     integer n, incx, incy
     double complex retval, zx(*), zy(*)
     external double complex zdotc

     retval = zdotc(n, zx, incx, zy, incy)

     return
     end
</code>

Made a shared object with

R CMD SHLIB --output=dozdot.so callzdotc.f czdot.c 

and ran

dyn.load("dozdot.so")
.C("testzdotc")

with the result 0.0, 0.0

Same here.

Once I change the line

	external double complex zdotc

in your fortran subroutine to

	double complex zdotc

everything works fine and I get as result 14.0, 0.0.

It is long time ago that I was taught (and studied) the FORTRAN 77
standard.  But flipping through some books from that time I thing I
have some inkling on what is going on.  The "external" statement is not
needed here (seems to be used in the sense that C is using the
"external" statement).

But I would really like to hear from an Rexpert why you
shouldn't/can't use external here in the Fortran.

G'day Berend,

On Tue, 6 Mar 2012 13:06:34 +0100
Berend Hasselman <bhh at xs4all.nl> wrote:

[... big snip ...]

But I would really like to hear from an Rexpert why you
shouldn't/can't use external here in the Fortran.

Probably less a question for an Rexpert but for a Fortran expert....

If you insert "implicit none" (one of my favourite extensions that I
always use) as the first statement after 
	subroutine callzdotc(retval,n, zx, incx, zy, incy)
you will see what is going on.  The compiler should refuse compilation
and complain that the type of zdotc was not declared (at least my
compiler does).  For FORTRAN to know that zdotc returns a double
complex you need the 
	double complex zdotc
declaration in callzdotc.

An
	external double complex zdotc
would be necessary if you were to call another subroutine/function, say
foo, that accepts functions as formal arguments and you want to pass
zdotc via such an argument.  Something like

	subroutine callzdotc(....)
       ....
       external double complex zdotc
	....
       call foo(a, b, zdotc)
       ....
	return
	end

	subroutine(a, b, h)
	double complex h, t
	....
	t = h(..,..,..,..)
	....
	return
	end

In C, the qualifier (or whatever the technical term is) "external" is
used to indicate that the function/variable/symbol is defined in
another compilation unit.  In FORTRAN77, "external" is used to tell the
compiler that you are passing a function to another
function/subroutine.  At least that is my understanding from what I
re-read in my FORTRAN documentation.

Thus, perhaps strangely, if there is only a 
	external double complex zdotc
declaration in your subroutine, the compiler doesn't know that a call
....

G'day Berend,

On Tue, 6 Mar 2012 13:06:34 +0100
Berend Hasselman<bhh at xs4all.nl>  wrote:

[... big snip ...]

But I would really like to hear from an Rexpert why you
shouldn't/can't use external here in the Fortran.

Probably less a question for an Rexpert but for a Fortran expert....

If you insert "implicit none" (one of my favourite extensions that I
always use) as the first statement after
	subroutine callzdotc(retval,n, zx, incx, zy, incy)
you will see what is going on.  The compiler should refuse compilation
and complain that the type of zdotc was not declared (at least my
compiler does).  For FORTRAN to know that zdotc returns a double
complex you need the
	double complex zdotc
declaration in callzdotc.

An
	external double complex zdotc
would be necessary if you were to call another subroutine/function, say
foo, that accepts functions as formal arguments and you want to pass
zdotc via such an argument.  Something like

	subroutine callzdotc(....)
         ....
         external double complex zdotc
	....
         call foo(a, b, zdotc)
         ....
	return
	end

	subroutine(a, b, h)
	double complex h, t
	....
	t = h(..,..,..,..)
	....
	return
	end

In C, the qualifier (or whatever the technical term is) "external" is

used to indicate that the function/variable/symbol is defined in
another compilation unit.  In FORTRAN77, "external" is used to tell the
compiler that you are passing a function to another
function/subroutine.  At least that is my understanding from what I
re-read in my FORTRAN documentation.

Thus, perhaps strangely, if there is only a
	external double complex zdotc
declaration in your subroutine, the compiler doesn't know that a call

to zdotc will return a double complex but will assume that the result
has the implicitly defined type, a real*8 IIRC.

puts a double complex as result on the stack (heap?), but within
callzdotc a real as return is expected and that is taken from the
stack (heap?), that real is than coerced to a double complex when
assigned to retval.  Note that while I am pretty sure about the above,
this last paragraph is more speculative. :)  But it would explain why
the erroneous code returns 0 on little-endian machines.

[..]
used to indicate that the function/variable/symbol is defined in
another compilation unit.  In FORTRAN77, "external" is used to tell the
compiler that you are passing a function to another
function/subroutine.  At least that is my understanding from what I
re-read in my FORTRAN documentation.

Not quite.  It also means that you really want to externally link and not allow the compiler to find any routine of that name it knows about (e.g. an intrinsic).  See para 8.7 of http://www.fortran.com/F77_std/rjcnf-8.html (although I got this from my Fortran reference on my bookshelf).

Thus, perhaps strangely, if there is only a
	external double complex zdotc
declaration in your subroutine, the compiler doesn't know that a call

The only 'strangely' thing is that is accepted: AFAICS is it not valid according to the link above.

to zdotc will return a double complex but will assume that the result
has the implicitly defined type, a real*8 IIRC.

The Fortran default type for z* is real.

So zdotc is called, and
puts a double complex as result on the stack (heap?), but within
callzdotc a real as return is expected and that is taken from the
stack (heap?), that real is than coerced to a double complex when
assigned to retval.  Note that while I am pretty sure about the above,
this last paragraph is more speculative. :)  But it would explain why
the erroneous code returns 0 on little-endian machines.

We haven't been told which machines, and this actually matters down to the level of optimization used by the compilers.  But these days typically double complex gets passed in sse registers, and double is passed in fpu registers.

Thus, perhaps strangely, if there is only a
	external double complex zdotc
declaration in your subroutine, the compiler doesn't know that a call

The only 'strangely' thing is that is accepted: AFAICS is it not valid according to the link above.

Agree. The fortran 77 standard doesn't allow that syntax and I'm really surprised that no error is flagged.

On 06/03/2012 13:37, Berwin A Turlach wrote:

G'day Berend,

On Tue, 6 Mar 2012 13:06:34 +0100
Berend Hasselman<bhh at xs4all.nl> ?wrote:

[... big snip ...]

But I would really like to hear from an Rexpert why you
shouldn't/can't use external here in the Fortran.

Probably less a question for an Rexpert but for a Fortran expert....

Exactly ....

If you insert "implicit none" (one of my favourite extensions that I
always use) as the first statement after
? ? ? ?subroutine callzdotc(retval,n, zx, incx, zy, incy)
you will see what is going on. ?The compiler should refuse compilation
and complain that the type of zdotc was not declared (at least my
compiler does). ?For FORTRAN to know that zdotc returns a double
complex you need the
? ? ? ?double complex zdotc
declaration in callzdotc.

An
? ? ? ?external double complex zdotc
would be necessary if you were to call another subroutine/function, say
foo, that accepts functions as formal arguments and you want to pass
zdotc via such an argument. ?Something like

? ? ? ?subroutine callzdotc(....)
? ? ? ? ....
? ? ? ? external double complex zdotc
? ? ? ?....
? ? ? ? call foo(a, b, zdotc)
? ? ? ? ....
? ? ? ?return
? ? ? ?end

? ? ? ?subroutine(a, b, h)
? ? ? ?double complex h, t
? ? ? ?....
? ? ? ?t = h(..,..,..,..)
? ? ? ?....
? ? ? ?return
? ? ? ?end

In C, the qualifier (or whatever the technical term is) "external" is

'extern' in C?

used to indicate that the function/variable/symbol is defined in
another compilation unit. ?In FORTRAN77, "external" is used to tell the
compiler that you are passing a function to another
function/subroutine. ?At least that is my understanding from what I
re-read in my FORTRAN documentation.

Not quite. ?It also means that you really want to externally link and not
allow the compiler to find any routine of that name it knows about (e.g. an
intrinsic). ?See para 8.7 of http://www.fortran.com/F77_std/rjcnf-8.html
(although I got this from my Fortran reference on my bookshelf).

Thus, perhaps strangely, if there is only a
? ? ? ?external double complex zdotc
declaration in your subroutine, the compiler doesn't know that a call

The only 'strangely' thing is that is accepted: AFAICS is it not valid
according to the link above.

to zdotc will return a double complex but will assume that the result
has the implicitly defined type, a real*8 IIRC.

The Fortran default type for z* is real.

So zdotc is called, and

puts a double complex as result on the stack (heap?), but within
callzdotc a real as return is expected and that is taken from the
stack (heap?), that real is than coerced to a double complex when
assigned to retval. ?Note that while I am pretty sure about the above,
this last paragraph is more speculative. :) ?But it would explain why
the erroneous code returns 0 on little-endian machines.

We haven't been told which machines, and this actually matters down to the
level of optimization used by the compilers. ?But these days typically
double complex gets passed in sse registers, and double is passed in fpu
registers.

Note that passing double complex/Rcomplex as return values, from C or
Fortran, is rather fragile in that it does depend on a pretty precise match
of compilers (and R's configure does test the constructions it uses, and
from time to time finds combinations which fail -- R 2.12.2 was released to
workaround one of them).

--
Brian D. Ripley, ? ? ? ? ? ? ? ? ?ripley at stats.ox.ac.uk
Professor of Applied Statistics, ?http://www.stats.ox.ac.uk/~ripley/
University of Oxford, ? ? ? ? ? ? Tel: ?+44 1865 272861 (self)
1 South Parks Road, ? ? ? ? ? ? ? ? ? ? +44 1865 272866 (PA)
Oxford OX1 3TG, UK ? ? ? ? ? ? ? ?Fax: ?+44 1865 272595

G'day Dominick,

On Mon, 5 Mar 2012 19:21:01 -0500
Dominick Samperi<djsamperi at gmail.com>  wrote:

This is consistent with the interface prototype for the BLAS
routine zdotc contained in<R>/include/R_ext/BLAS.h, namely,

BLAS_extern Rcomplex
     F77_NAME(zdotc)(Rcomplex * ret_val, int *n,
		    Rcomplex *zx, int *incx, Rcomplex *zy, int *incy);

[...]

On the other hand, this is not consistent with the standard
FORTRAN definition for zdotc that is contained in
<R>/src/extra/blas/cmplxblas.f, where the first argument is
n, not ret_val.

This seems to be indeed inconsistent and, presumably, a bug.  Applying
the attach patch to R's development version (compiles, installs and
passes all checks with this patch), and changing in your code the line

	F77_CALL(zdotc)(&ret_val,&n, zx,&incx, zy,&incy);

to

	ret_val = F77_CALL(zdotc)(&n, zx,&incx, zy,&incy);

produces the expected output.