We distribute, as source-code, a large package of numerically-focused
programs, which may be compiled
for windows by Lahey compilers LF90, LF95, and Intel IF32 and IF64. We
aim now to add gfortran to that list.
Only windows 32 or 64 (and not Itanium) is of interest at the moment.

Overflow errors are unavoidable, and we try to to inform the user
whence they originated. This is done by
disabling floating point exceptions, and examining the "8087 status
word" after loops (not after every operation).
The different fortrans listed above offer different extensions or
methods to accomplish this.

Is there a similar thing for GFortran ?

I see that we can affect  the "8087 control word" via compile option -
ffpe-trap (so that overflow triggers no exception), but see no way to
detect (via "8087 status word") whether an overflow recently
occurred.

This "8087 status word" is of course an old-fashioned way of speaking.
The facility I seek
(maybe called Get_IEEE_Flags) works for SSE as well, but uses a
structure akin to the old
87 status word.

I see, in  _mingw_float.h mention of "statusfp" and "fpecode".  Are
these accessible to a Gfortran program?
I would prefer if possible not to write and link in a C routine to do
this job.
My desired solution for now can be specific to windows Gfortran -- I
understand that a more general, standards-based, solution may be in
the pipeline, viz: http://gcc.gnu.org/ml/gcc/2008-11/msg00372.html

Mark Horridge

On 11/11/2010 6:12 AM, mark.horridge@buseco.monash.edu.au wrote:
> We distribute, as source-code, a large package of numerically-focused
> programs, which may be compiled
> for windows by Lahey compilers LF90, LF95, and Intel IF32 and IF64. We
> aim now to add gfortran to that list.
> Only windows 32 or 64 (and not Itanium) is of interest at the moment.
>
> Overflow errors are unavoidable, and we try to to inform the user
> whence they originated. This is done by
> disabling floating point exceptions, and examining the "8087 status
> word" after loops (not after every operation).
> The different fortrans listed above offer different extensions or
> methods to accomplish this.
>
> Is there a similar thing for GFortran ?
>
> I see that we can affect  the "8087 control word" via compile option -
> ffpe-trap (so that overflow triggers no exception), but see no way to
> detect (via "8087 status word") whether an overflow recently
> occurred.
>
> This "8087 status word" is of course an old-fashioned way of speaking.
> The facility I seek
> (maybe called Get_IEEE_Flags) works for SSE as well, but uses a
> structure akin to the old
> 87 status word.
>
> I see, in  _mingw_float.h mention of "statusfp" and "fpecode".  Are
> these accessible to a Gfortran program?
> I would prefer if possible not to write and link in a C routine to do
> this job.
> My desired solution for now can be specific to windows Gfortran -- I
> understand that a more general, standards-based, solution may be in
> the pipeline, viz: http://gcc.gnu.org/ml/gcc/2008-11/msg00372.html
>
> Mark Horridge
I link in my own code for support of a portion of ieee_arithmetic, using 
a combination of fenv.h and IA intrinsics, via iso_c_interop.
  Implementing enough of ieee_arithmetic to satisfy more than one or two 
people on more than than the IA architectures is a daunting task. Even 
that might not approach the facility you want.  Compilers with 
ieee_arithmetic are available on free trial, so you can test what can be 
done with it, and make your additions to gfortran work the same.
  If it were implemented in the compiler, you would still be linking in 
C code, but the job would be done for you.  It's certainly feasible to 
append your own code to libgfortran.
I doubt that anyone would want to tackle a Windows-only solution.  My 
partial solution works equally well with gfortran on Windows and linux.

-- 
Tim Prince

In article <8k2d5rFb9eU1@mid.individual.net>,
Tim Prince  <tprince@nospamcomputer.org> wrote:
>>
>I link in my own code for support of a portion of ieee_arithmetic, using 
>a combination of fenv.h and IA intrinsics, via iso_c_interop.

You do like to live dangerously, don't you? :-)  I would expect that
to fail horribly under some circumstances, though you do imply that
you are using very little of the functionality and only a couple
of compilers.

>  Implementing enough of ieee_arithmetic to satisfy more than one or two 
>people on more than than the IA architectures is a daunting task. Even 
>that might not approach the facility you want.  Compilers with 
>ieee_arithmetic are available on free trial, so you can test what can be 
>done with it, and make your additions to gfortran work the same.

Gug.  It's a daunting task on IA architectures if you actually want
to get it right.

>  If it were implemented in the compiler, you would still be linking in 
>C code, but the job would be done for you.  It's certainly feasible to 
>append your own code to libgfortran.

No, you wouldn't.  A Fortran compiler writer would be insane to do
that.  Invoking the hardware facilities directly is FAR safer.

>I doubt that anyone would want to tackle a Windows-only solution.  My 
>partial solution works equally well with gfortran on Windows and linux.

As it's entirely a matter of handshaking between the compiler and
hardware, the operating system is almost irrelevant.


Regards,
Nick Maclaren.

On Nov 11, 3:33=A0pm, Tim Prince <tpri...@computer.org> wrote:
> On 11/11/2010 6:12 AM, mark.horri...@buseco.monash.edu.au wrote:
> I link in my own code for support of a portion of ieee_arithmetic, using
> a combination of fenv.h and IA intrinsics, via iso_c_interop.

Is this module somewhere available? I only found a version which
supports setting/getting the underflow and rounding mode:
http://gcc.gnu.org/ml/fortran/2010-04/msg00144.html


> Implementing enough of ieee_arithmetic to satisfy more than one or two
> people on more than than the IA architectures is a daunting task. Even
> that might not approach the facility you want.

I think work on adding the Fortran 2003 IEEE support to gfortran will
start for GCC/gfortran 4.7. The number of people expecting that this
feature is present [1] is growing -- while the number of other
outstanding Fortran 2003 tasks is shrinking [2]. However, implementing
a fully working Fortran 2003 IEEE support* for all GCC architectures
will be larger job and take a while.

Tobias

[1] That's also due to the relatively large number of other compilers,
which have this feature already longer. Cf.
http://fortranwiki.org/fortran/show/Fortran+2003+status

[2] Cf. [1] and http://gcc.gnu.org/wiki/Fortran2003Status ; while a
lot has been done, a lot still needs to be done: IEEE, UDIO, type
parameters and unlimited polymorphism. Plus F2008's submodules and
coarrays (http://gcc.gnu.org/wiki/Fortran2008Status).

* Whatever one understands under 'fully working'. Like with
interoperability and parallel execution (e.g. using POSIX threads),
there is a certain fuzziness and for things which fail under certain
circumstances.

<mark.horridge@buseco.monash.edu.au> wrote in message 
news:ee320077-1185-43b8-9a3e-c16485abaf8f@b19g2000prj.googlegroups.com...

> We distribute, as source-code, a large package of numerically-focused
> programs, which may be compiled
> for windows by Lahey compilers LF90, LF95, and Intel IF32 and IF64. We
> aim now to add gfortran to that list.
> Only windows 32 or 64 (and not Itanium) is of interest at the moment.

> Overflow errors are unavoidable, and we try to to inform the user
> whence they originated. This is done by
> disabling floating point exceptions, and examining the "8087 status
> word" after loops (not after every operation).
> The different fortrans listed above offer different extensions or
> methods to accomplish this.

> Is there a similar thing for GFortran ?

> I see that we can affect  the "8087 control word" via compile option -
> ffpe-trap (so that overflow triggers no exception), but see no way to
> detect (via "8087 status word") whether an overflow recently
> occurred.

> This "8087 status word" is of course an old-fashioned way of speaking.
> The facility I seek
> (maybe called Get_IEEE_Flags) works for SSE as well, but uses a
> structure akin to the old
> 87 status word.

> I see, in  _mingw_float.h mention of "statusfp" and "fpecode".  Are
> these accessible to a Gfortran program?
> I would prefer if possible not to write and link in a C routine to do
> this job.
> My desired solution for now can be specific to windows Gfortran -- I
> understand that a more general, standards-based, solution may be in
> the pipeline, viz: http://gcc.gnu.org/ml/gcc/2008-11/msg00372.html

In 64-bit Windows you want to read and write MXCSR because it doesn't
use the x87 floating point stack.  In 32-bit Windows, there is a
switch to force SSE instead of x87, so all that is required is a
standard-conforming Fortran program the can read MXCSR (via the
STMXCSR instruction) and write MXCSR (via the LDMXCSR instruction).
This we can do in Windows:

C:\gfortran\clf\mxcsr>type mxcsr.f90
module mxcsr
   use ISO_C_BINDING
   implicit none
   private
   public mxcsr_init, ldmxcsr, stmxcsr
   integer(C_INT8_T) BAD_STUFF(56)
   data BAD_STUFF / &
      Z'31', Z'C0', Z'48', Z'85', Z'C0', Z'75', Z'0A', Z'89', &
      Z'4C', Z'24', Z'08', Z'0F', Z'AE', Z'54', Z'24', Z'08', &
      Z'C3', Z'8B', Z'44', Z'24', Z'04', Z'50', Z'0F', Z'AE', &
      Z'14', Z'24', Z'58', Z'C3', Z'90', Z'90', Z'90', Z'90', &
      Z'31', Z'C0', Z'48', Z'85', Z'C0', Z'75', Z'0A', Z'0F', &
      Z'AE', Z'5C', Z'24', Z'08', Z'8B', Z'44', Z'24', Z'08', &
      Z'C3', Z'50', Z'0F', Z'AE', Z'1C', Z'24', Z'58', Z'C3' /
   abstract interface
      subroutine ldmxcsr_template(source) bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT32_T), value :: source
      end subroutine ldmxcsr_template
      function stmxcsr_template() bind(C)
         use ISO_C_BINDING
         implicit none
         integer(C_INT32_T) :: stmxcsr_template
      end function stmxcsr_template
   end interface
   procedure(ldmxcsr_template), pointer :: ldmxcsr
   procedure(stmxcsr_template), pointer :: stmxcsr
   interface
      function VirtualAlloc(lpAddress, dwSize, flAllocationType, &
            flProtect) bind(C, name = 'VirtualAlloc')
         use ISO_C_BINDING
         implicit none
!GCC$ ATTRIBUTES STDCALL :: VirtualAlloc
         type(C_PTR) VirtualAlloc
         type(C_PTR), value :: lpAddress
         integer(C_SIZE_T), value :: dwSize
         integer(C_LONG), value :: flAllocationType
         integer(C_LONG), value :: flProtect
      end function VirtualAlloc
      function GetLastError() bind(C,name='GetLastError')
         use ISO_C_BINDING
         implicit none
!GCC$ ATTRIBUTES STDCALL :: GetLastError
         integer(C_LONG) GetLastError
      end function GetLastError
   end interface
   contains
      subroutine mxcsr_init()
         type(C_PTR) address
         integer(C_INTPTR_T) temp
         integer(C_LONG) error
         integer(C_INT8_T), pointer :: temp_ptr(:)
         type(C_FUNPTR) fun_address

         address = VirtualAlloc(C_NULL_PTR, &
            size(BAD_STUFF,kind=C_SIZE_T),int(Z'1000',C_LONG), &
            int(Z'40',C_LONG))
         if(.NOT.C_ASSOCIATED(address)) then
            error = GetLastError()
            write(*,'(a,z0,a)') "Error Z'",error,"' allocating memory"
            stop
         end if
         call C_F_POINTER(address, temp_ptr, shape(BAD_STUFF))
         temp_ptr = BAD_STUFF
         temp = transfer(address, temp)
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, ldmxcsr)
         temp = temp+32
         fun_address = transfer(temp, fun_address)
         call C_F_PROCPOINTER(fun_address, stmxcsr)
      end subroutine mxcsr_init
end module mxcsr

program mxcsr_test
   use ISO_C_BINDING
   use mxcsr
   implicit none
   integer(C_INT32_T) reg

   call mxcsr_init
   reg = stmxcsr()
   write(*,'(z0.8)') reg
   reg = ibclr(reg,9)
   write(*,'(z0.8)') reg
   call ldmxcsr(reg)
   reg = stmxcsr()
   write(*,'(z0.8)') reg
end program mxcsr_test

C:\gfortran\clf\mxcsr>gfortran -fno-range-check mxcsr.f90 -omxcsr

C:\gfortran\clf\mxcsr>mxcsr
00001F80
00001D80
00001D80

See how at the start the program read the default value of
Z'00001F80' from MXCSR and then we took the value Z'00001D80'
and wrote it into MXCSR then read it back again?

So just look up the switch to force SSE in 32-bit Windows and
look in 253665.pdf, section 10.2.3 to see what to do with the
MXCSR register and you're good to go!

Oh, I forgot to insert a flag to prevent a second invocation of
mxcsr_init.  Sorry.

-- 
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

On 11/11/2010 8:29 AM, Tobias Burnus wrote:
> On Nov 11, 3:33 pm, Tim Prince<tpri...@computer.org>  wrote:
>> On 11/11/2010 6:12 AM, mark.horri...@buseco.monash.edu.au wrote:
>> I link in my own code for support of a portion of ieee_arithmetic, using
>> a combination of fenv.h and IA intrinsics, via iso_c_binding
>
> Is this module somewhere available? I only found a version which
> supports setting/getting the underflow and rounding mode:
> http://gcc.gnu.org/ml/fortran/2010-04/msg00144.html
>
http://sites.google.com/site/tprincesite/Home/gfortran-ieee-arithmetic
>

-- 
Tim Prince