[Rd] Fortran BLAS giving bad results

[Emmanuel Sharef]

Hello r-devel,

When compiling Fortran code containing BLAS functions and calling it using
dyn.load, I am getting incorrect results. A small example with which I can
reproduce the problem is below.
I am running on OSX Mavericks (upgraded R, Xcode, etc per instructions on
this list), but I do not think it is related since the check on my
"blupsurv" package on r-forge seems to be showing similar behavior, and the
problem predates Mavericks. The package used to build and check without
errors when I originally put it on r-forge. Would appreciate suggestions.

Thanks,
Emmanuel

Example: this program creates two double vectors, takes the dot product
with ddot, and prints the result:

test.f:
  subroutine testddot(n,x,y)
  integer i,n
  double precision x,y,d1(n),d2(n)
  do 100 i=1,n
d1(i)=x
d2(i)=y
 100  continue
  print *,ddot(n,d1,1,d2,1)
  end

Compiling
R CMD SHLIB test.f

Running:
> dyn.load("test.so")
> test<-.Fortran("testddot",as.integer(3),as.double(1),as.double(1.2))
 -1.0737418E+08

The result should be 3.6, not -1.07e8. Most other inputs give incorrect
results too.

Here is my session info. I upgraded to the latest devel to see if it would
fix the problem.

> sessionInfo()
R Under development (unstable) (2014-01-10 r64727)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Fortran BLAS giving bad results

[Berwin A Turlach]

G'day Emmanuel,

On Sat, 11 Jan 2014 14:08:41 -0800
Emmanuel Sharef  wrote:

> Example: this program creates two double vectors, takes the dot
> product with ddot, and prints the result:
> 
> test.f:
>   subroutine testddot(n,x,y)
>   integer i,n
>   double precision x,y,d1(n),d2(n)
>   do 100 i=1,n
> d1(i)=x
> d2(i)=y
>  100  continue
>   print *,ddot(n,d1,1,d2,1)
>   end
> 
> Compiling
> R CMD SHLIB test.f

Mmh, FORTRAN77 does not have dynamical memory allocation (though I have
been told that many compilers implement it as an extension to the
standard).  So I am surprised that this code works.  It is certainly
not portable (i.e. a compiler that implements the FORTRAN77 standard
strictly would barf on this code).

Also, output to * in FORTRAN code dynamically linked to R is
discouraged, see chapter 5.7 (Fortran I/O) in the "Writing R Extension"
manual.  From memory, many years ago one could do so (on linux and
SUN machines atleast), then it stopped working.  I am surprised that it
works again.

But neither of these two issues are the reason for your strange
results.  You would have pinpointed your actual problem easily if you
would have used my favourite none-FORTRAN77 feature (which is
implemented by most compilers as an extension), namely having "implicit
none" as the first line in each function/subroutine. (The FORTRAN77
compliant version, AFAIK, was to write "IMPLICIT CHARACTER (A-Z)" into
the first line of all functions/subroutines.)

Essentially, your routine testddot does not know what kind of result
ddot returns and by the implicit type definition rules of FORTRAN77 it
expects an INTEGER back.  The fact that a DOUBLE PRECISION value is
returned when an INTEGER is expected creates some confusion.  Somewhere
at the start of your routine you have to add a "DOUBLE PRECISION ddot"
statement.

HTH.

Cheers,

Berwin

== Full address 
A/Prof Berwin A Turlach   Tel.: +61 (8) 6488 3338 (secr)
School of Maths and Stats (M019)+61 (8) 6488 3383 (self)
The University of Western Australia   FAX : +61 (8) 6488 1028
35 Stirling Highway   
Crawley WA 6009 e-mail: berwin.turl...@gmail.com
Australiahttp://www.maths.uwa.edu.au/~berwin
 http://www.researcherid.com/rid/A-4995-2008

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel