From djsamperi at gmail.com  Tue Mar  6 16:39:32 2012
From: djsamperi at gmail.com (Dominick Samperi)
Date: Tue, 6 Mar 2012 10:39:32 -0500
Subject: [Rd] Calling FORTRAN function from R issue?
In-Reply-To: <4F561C79.2090002@stats.ox.ac.uk>
References: <CADUbQ5ghmCd-3K4hx=FrzhA8YwfohhUp0-g6rJHHU0s-uEfN+A@mail.gmail.com>
	<B1B70BE3-2D8A-495B-8893-CC67026E1E9F@xs4all.nl>
	<20120306195608.7ec858c0@bossiaea>
	<879BAEC9-1963-4574-B962-CA8CC9D4F13B@xs4all.nl>
	<20120306213703.6a5a091f@absentia>
	<4F561C79.2090002@stats.ox.ac.uk>
Message-ID: <CADUbQ5gmAcy99MS0hGSmmzx0V1bAbDh29_oZLgJVN9OkegHFMw@mail.gmail.com>

On Tue, Mar 6, 2012 at 9:17 AM, Prof Brian Ripley <ripley at stats.ox.ac.uk> 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<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).

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, ? ? ? ? ? ? ? ? ?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
>
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



[Rd] Calling FORTRAN function from R issue?
Dominick Samperi