[Rd] Fortran BLAS giving bad results

Berend Hasselman bhh at xs4all.nl
Sun Jan 12 11:18:05 CET 2014


On 12-01-2014, at 08:26, Berwin A Turlach <Berwin.Turlach at gmail.com> wrote:

> G'day Emmanuel,
> 
> On Sat, 11 Jan 2014 14:08:41 -0800
> Emmanuel Sharef <esharef at gmail.com> 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.


No. It is expecting ddot to return a REAL result (single precision).

Berend

>  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.Turlach at gmail.com
> Australia                http://www.maths.uwa.edu.au/~berwin
>                         http://www.researcherid.com/rid/A-4995-2008
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



More information about the R-devel mailing list