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/
