On 06/03/2012 13:37, Berwin A Turlach wrote:
G'day Berend,
On Tue, 6 Mar 2012 13:06:34 +0100
Berend Hasselman<b...@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, rip...@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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
