Message-ID: <489769E7.8060207@alice.it>
Date: Mon, 04 Aug 2008 22:43:19 +0200
From: Angelo Graziosi <angelo DOT graziosi AT alice DOT it>
User-Agent: Thunderbird 2.0.0.16 (Windows/20080708)
MIME-Version: 1.0
To: cygwin AT cygwin DOT com
Subject: RE: cygwin gcc: Different numerical results in thread vs in main()
Content-Type: text/plain; charset=ISO-8859-15; format=flowed
Content-Transfer-Encoding: 7bit
Mailing-List: contact cygwin-help AT cygwin DOT com; run by ezmlm
Sender: cygwin-owner AT cygwin DOT com
Mail-Followup-To: cygwin AT cygwin DOT com

Dave Korn wrote:

> Take a look at the (legendary) GCC PR323:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
> 
> ... and in particular comment #60:
> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323#c60
> 
> ... which has a bit of code you can adapt (google for the definition of
> _FPU_SETCW, it's an inline asm you can copy from linux/x86 headers); using
> that to set double precision mode in your thread startup code, and using
> double (never float) types throughout your code is probably the best you can
> do.

... and in Fortran?

I have flagged a similar problem on fortran list a few weeeks ago [1]. 
And obviously they pointed out the same solution [2,3].

Now rereading '...323#c60', this simple solution in F2003 comes in mind:

$ cat set_math_double_precision.c
/*
   We cannot use

   #include <fpu_control.h>

   because Cygwin haven't that header. It can be downloaded from

http://boinc.gorlaeus.net/download/DownLoads/Classical/crlibm_libs/freebsd/includes/fpu_control.h
  */
#include "fpu_control.h"

void set_math_double_precision() {
    fpu_control_t fpu_control = 0x027f ;
    _FPU_SETCW(fpu_control);
}

$ cat test_case.f90
!
! gfortran -c set_math_double_precision.c
!
! (or gcc4 -c set_math_double_precision.c)
!
! gfortran test_case.f90 -o test_case set_math_double_precision.o
!
! See: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
! And also: http://cygwin.com/ml/cygwin/2008-08/msg00075.html,
!           http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323#c60
!

module fpu95
   interface
      subroutine set_math_double_precision() &
           bind(C,name='set_math_double_precision')
        use, intrinsic :: iso_c_binding
      end subroutine set_math_double_precision
   end interface
end module fpu95

program test_case
   use fpu95
   implicit none
   integer :: k
   integer, parameter :: DP = kind(1.D0),&
        N = 29
   real(DP), parameter :: A(0:N) = &
        (/1.0000000000000D0,3.8860009363293D0,7.4167881843083D0,&
          9.8599837415463D0,10.431383465276D0,9.4077161304727D0,&
          7.5332956276775D0,5.4995844630326D0,3.7275104474917D0,&
          2.3766149078263D0,1.4394444941337D0,0.8344437543616D0,&
          0.4657058294147D0,0.2514185849207D0,0.1317920740387D0,&
          0.0672995854816D0,0.0335554996094D0,0.0163785206816D0,&
          0.0078368645385D0,0.0036846299276D0,0.0017020281304D0,&
          0.0007756352209D0,0.0003469750783D0,0.0001544224626D0,&
          0.0000666379875D0,0.0000295655909D0,0.0000118821415D0,&
          0.0000057747681D0,0.0000017502652D0,0.0000014682034D0/)
   real(DP), parameter :: BB(0:N-1) = &
        (/(((k+2)*A(k+2)-A(k+1)*(A(1)+(k+1))),k = 0,N-2),&
          (-A(N)*(A(1)+N))/)
   real(DP) :: b(0:N-1),bbb(0:N-1)
   call set_math_double_precision()
   do k = 0,N-2
      b(k) = (k+2)*A(k+2)-A(k+1)*(A(1)+(k+1))
   enddo
   b(N-1) = -A(N)*(A(1)+N)
   bbb(0:N-1) = (/(((k+2)*A(k+2)-A(k+1)*(A(1)+(k+1))),k = 0,N-2),&
        (-A(N)*(A(1)+N))/)
   do k = 0,N-1
      print *, k,b(k)-BB(k),b(k)-bbb(k),BB(k)-bbb(k)
   enddo
end program test_case

$ gfortran -c set_math_double_precision.c
$ gfortran test_case.f90 -o test_case set_math_double_precision.o

Running './test_case' it prints all zeroes!


Cheers,
    Angelo.

---
[1] http://gcc.gnu.org/ml/fortran/2008-06/msg00065.html
[2] http://gcc.gnu.org/ml/fortran/2008-06/msg00118.html
[3] http://gcc.gnu.org/ml/fortran/2008-06/msg00119.html

--
Unsubscribe info:      http://cygwin.com/ml/#unsubscribe-simple
Problem reports:       http://cygwin.com/problems.html
Documentation:         http://cygwin.com/docs.html
FAQ:                   http://cygwin.com/faq/