This is the mail archive of the cygwin mailing list for the Cygwin project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

RE: cygwin gcc: Different numerical results in thread vs in main()


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/


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]