Quantcast

Calling FORTRAN function from R issue?

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Calling FORTRAN function from R issue?

Dominick Samperi-3
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.

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);

But this seems to BOTH return a result, and pass the result
as the first argument?

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.

My search is complicated by the fact that the libraries
libR.so, libRblas.so, libRlapack.so are stripped.

When I install the standard (FORTRAN-based) BLAS on my
system (Fedora 16), I find that zdotc is defined in the traditional
way (without adjustment to the first argument). This is probably
irrelevant because R does not use it in my configuration.

I found some documentation on Intel FORTRAN that seems to
suggest that the first argument on the C side is always the
same as (a pointer to) the FORTRAN function return value, but
this is not so if I use the standard definition of zdotc.f with
gfortran.

Any ideas?

Thanks,
Dominick

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berwin A Turlach-3
G'day Dominick,

On Mon, 5 Mar 2012 19:21:01 -0500
Dominick Samperi <[hidden email]> wrote:

> 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.

Section 5.2 deals with calling C/FORTRAN code from R via .C()
or .Fortran(), and is not directly relevant to the question on how to
call FORTRAN code from C code. :)
 

> 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.

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

The F77_xxxx macros seem to be defined in <R>/include/R_ext/RS.h, and
their sole purpose seems to be to append a _ to the argument if
needed.  

Cheers,

        Berwin

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

R-Patch (1K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berend Hasselman
In reply to this post by Dominick Samperi-3

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

> 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.

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

Above C routine changed to

<code>
#include "R.h"

void F77_NAME(callzdotc)(Rcomplex *, int *, Rcomplex *, int *, Rcomplex *, int *);

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(callzdotc)(&ret_val, &n, zx, &incx, zy, &incy);
   Rprintf("ret_val = %f, %f\n", ret_val.r, ret_val.i);
}
</code>

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

I would've expected this to give the correct result.

Berend

Mac OS X 10.6.8
R2.14.2
Using reference Rblas.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berwin A Turlach-3
G'day Berend,

On Tue, 6 Mar 2012 11:19:07 +0100
Berend Hasselman <[hidden email]> wrote:

> 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;

Just noticing that it is always zx[0].i, same with the imaginary part
of zy.  But this is probably not of importance. :)

> 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).

Cheers,

        Berwin

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berend Hasselman

On 06-03-2012, at 12:56, Berwin A Turlach wrote:

> 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).


Thanks.
I should have tried that too.

This implies that Dominick's original issue can be avoided by using an intermediate Fortran routine.

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

Berend

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berwin A Turlach-3
G'day Berend,

On Tue, 6 Mar 2012 13:06:34 +0100
Berend Hasselman <[hidden email]> 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.  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.

HTH.

Cheers,
       
        Berwin




Cheers,

        Berwin

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berend Hasselman

On 06-03-2012, at 14:37, Berwin A Turlach wrote:

> G'day Berend,
>
> On Tue, 6 Mar 2012 13:06:34 +0100
> Berend Hasselman <[hidden email]> 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.
>

I usually use -fimplicit-none or a similar option as argument for a fortran compiler.
I didn't use it in a syntax checking pre compiling run.

> 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.

Not quite true.
It is required to use external if you are passing the name of a function/subroutine to another routine.
Otherwise it is just a statement telling the compile that the <external-name> is external to the function/subroutine/.. being compiled.

See http://www.fortran.com/F77_std/rjcnf0001-sh-8.html#sh-8.7


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

If I am reading the above reference correctly, the syntax

        external double complex zdotc

is simply incorrect (correct syntax uses commas to separate names e.g. A,B,C).
The compiler (gfortran) simply seems to ignore the "double" and the "complex" if you don't use implicit none.

Berend

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Prof Brian Ripley
In reply to this post by Berwin A Turlach-3
On 06/03/2012 13:37, Berwin A Turlach wrote:

> G'day Berend,
>
> On Tue, 6 Mar 2012 13:06:34 +0100
> Berend Hasselman<[hidden email]>  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,                  [hidden email]
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berend Hasselman

On 06-03-2012, at 15:17, Prof Brian Ripley wrote:

>> [..]
>> 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.
>

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

>> 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.
>

Mac OS X 10.6.8, Mini Core 2 Duo
Apple gcc (i686-apple-darwin10-gcc-4.2.1 (GCC) 4.2.1 (Apple Inc. build 5664))
and gfortran (GNU Fortran (GCC) 4.2.1 (Apple Inc. build 5664)) from r.research.att.com.

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

gfortran -arch x86_64   -fPIC  -g -O2 -c callzdotc.f -o callzdotc.o
gcc -arch x86_64 -std=gnu99 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/x86_64  -I/usr/local/include    -fPIC  -g -O2 -c czdot.c -o czdot.o
gcc -arch x86_64 -std=gnu99 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/lib -o dozdot.so callzdotc.o czdot.o -lgfortran -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation

Same thing happens with gcc (Ubuntu/Linaro 4.6.1-9ubuntu3) 4.6.1 and gfortran on Xubuntu 11.10 64-bit.

Berend

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Berend Hasselman

On 06-03-2012, at 15:46, Berend Hasselman wrote:

>
>>
>>> 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.
>

Got it.
See http://www.fortran.com/F77_std/rjcnf0001-sh-3.html#sh-3.1.6

Compiler ignores the blanks and

        external double complex zdotc

becomes

        external doublecomplexzdotc

And that is legal.

And how long have I been using Fortran? I won't tell.

Berend

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Dominick Samperi-3
In reply to this post by Prof Brian Ripley
On Tue, Mar 6, 2012 at 9:17 AM, Prof Brian Ripley <[hidden email]> wrote:

> On 06/03/2012 13:37, Berwin A Turlach wrote:
>>
>> G'day Berend,
>>
>> On Tue, 6 Mar 2012 13:06:34 +0100
>> Berend Hasselman<[hidden email]>  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).

It turns the problem reported above was discovered in a completely
different context, unrelated to R, and when I looked at how R handles
the problem I found similar issues. As Brian says, the C/Fortran
interface can be fragile and machine/compiler-dependent.

The code appended below may shed some light on the issues. The
Fortran code provides three interfaces to zdotc, two function interfaces
and one subroutine interface. The test driver exercises each interface
and the results are what you would expect.

This was tested using gfortran/gcc under Fedora 16 (64bit).

This depends on the particular way the gfortran and gcc
interact, specifically, the -std=gnu flag is important (it is
also the default). Also, it is important to remember that
REAL in Fortran maps to float in C, REAL*16 maps to
double, DOUBLE COMPLEX maps to complex, etc.

To run the executable you may have to use:
export LD_LIBRARY_PATH=.

*** Makefile ***
testzdotc:
        gfortran -std=gnu -fPIC -shared -o libmyblas.so zdotc.f
        gcc -o testzdotc testzdotc.c -L/home/dsamperi -lmyblas -lm

*** testzdotc.c ***

#include <stdio.h>
#include <complex.h>

extern double complex zdotc1_(int*,complex*,int*,complex*,int*);
extern void zdotc2_(complex*,int*,complex*,int*,complex*,int*);
extern float zdotc3_(int*,complex*,int*,complex*,int*);

int main() {
    int n=3;
    int incx = 1;
    int incy = 1;
    complex zx[3], zy[3], zdotc,zdotc2;
    double z;
    printf("Size complex: %d\n", sizeof(complex));
    printf("Size double = %d\n", sizeof(double));
    printf("Size float:   %d\n", sizeof(float));
    zx[0] = 1+0*I; zx[1] = 2+0*I; zx[2] = 3+0*I;
    zy[0] = 1+0*I; zy[1] = 2+0*I; zy[2] = 3+0*I;
    zdotc = zdotc1_(&n,zx,&incx,zy,&incy);
    zdotc2_(&zdotc2,&n,zx,&incx,zy,&incy);
    z = zdotc3_(&n,zx,&incx,zy,&incy);
    printf("Using zdotc1: %g,%g\n", creal(zdotc),cimag(zdotc));
    printf("Using zdotc2: %g,%g\n", creal(zdotc2),cimag(zdotc2));
    printf("Using zdotc3: %f\n", z);
}

*** zdotc.f ***

** Three variants of zdotc based on the BLAS code of jack dongarra (3/11/78).
      DOUBLE COMPLEX FUNCTION ZDOTC1(N,ZX,INCX,ZY,INCY)
      INTEGER INCX,INCY,N
      DOUBLE COMPLEX ZX(*),ZY(*)
      DOUBLE COMPLEX ZTEMP
      INTEGER I,IX,IY
      INTRINSIC DCONJG
      ZTEMP = (0.0d0,0.0d0)
      ZDOTC = (0.0d0,0.0d0)
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
         DO I = 1,N
            ZTEMP = ZTEMP + DCONJG(ZX(I))*ZY(I)
         END DO
      ELSE
         IX = 1
         IY = 1
         IF (INCX.LT.0) IX = (-N+1)*INCX + 1
         IF (INCY.LT.0) IY = (-N+1)*INCY + 1
         DO I = 1,N
            ZTEMP = ZTEMP + DCONJG(ZX(IX))*ZY(IY)
            IX = IX + INCX
            IY = IY + INCY
         END DO
      END IF
      ZDOTC1 = ZTEMP
      RETURN
      END

      SUBROUTINE ZDOTC2(ZDOTC,N,ZX,INCX,ZY,INCY)
      EXTERNAL ZDOTC1
      INTEGER INCX,INCY,N
      DOUBLE COMPLEX ZX(*),ZY(*),ZDOTC, ZDOTC1
      ZDOTC = ZDOTC1(N,ZX,INCX,ZY,INCY)
      RETURN
      END

      REAL FUNCTION ZDOTC3(N,ZX,INCX,ZY,INCY)
      EXTERNAL ZDOTC1
      INTEGER INCX,INCY,N
      DOUBLE COMPLEX ZX(*),ZY(*), ZDOTC1
      ZDOTC3 = ZDOTC1(N,ZX,INCX,ZY,INCY)
      RETURN
      END

>
> --
> Brian D. Ripley,                  [hidden email]
> 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
>
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Calling FORTRAN function from R issue?

Prof Brian Ripley
In reply to this post by Berwin A Turlach-3
On 06/03/2012 06:28, Berwin A Turlach wrote:
> G'day Dominick,
>
> On Mon, 5 Mar 2012 19:21:01 -0500
> Dominick Samperi<[hidden email]>  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

As I said elsewhere in this thread, this really is very
compiler-specific, and rather than being a bug, that header is not
appropriate to the compilers used (it came from the days of f2c and
Fortran compilers based on it such as g77).

I'll change the sources to follow your patch as I think it is much more
likely to be correct these days, but also add a warning in the header.
I don't think it is safe to call these functions from C without
configure-testing the effect.

>
> 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.


--
Brian D. Ripley,                  [hidden email]
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Loading...