What is the current status of Quad-float support in Fortran?

[FX]

Dear Alan,

I am sending this reply to the gfortran developers and users mailing-list (as you suggested in your email). Thus, I keep your email in full below so people understand the context of the reply.

Quad-float (aka IEEE binary128) support was merged into mainline GCC thanks to the efforts of Tobias Burnus and the gfortran maintainers, and is part of the 4.6.x release branch (the latest release of this branch is 4.6.1, released on 2011-06-??). This is not so hard to build, and if you get into trouble there you can ask for help on the gcc-help@gcc.gnu.org mailing-list.

Regarding your -fdefault-real-8 questions, I will reply with a small code snippet (run here on x86_64 Mac OS):

 real a
 double precision b
 real*4 c
 real*8 d
 print *, sizeof(a), sizeof(b)
 print *, sizeof(2.), sizeof(2.d0)
 print *, sizeof(c), sizeof(d)
 print *, sizeof(2._4), sizeof(2._8)
 end


$ gfortran-4.5 a.f90 && ./a.out                 
                   4                    8
                   4                    8
                   4                    8
                   4                    8
$ gfortran-4.5 a.f90 -fdefault-real-8 && ./a.out
                   8                    8
                   8                    8
                   4                    8
                   4                    8
$ gfortran-4.6 a.f90 -fdefault-real-8 && ./a.out
                   8                   16
                   8                   16
                   4                    8
                   4                    8

So:
 -- with GCC < 4.6.0, -fdefault-real-8 is promotes "real" but not "double precision" (unless on a platform with native support for 128-bit long double, such as Mac OS on powerpc)
 -- with quad-float support, "double precision" is promoted to 128-bit
 -- literal constants without specific kind are promoted the same as "real" and "double precision"
 -- variables specified with specific width (real*8) or kind (real(kind=8)) are not promoted

I understand the last item on the list will disappoint you. I dispute your statement that "uniform promotion is the only version of that option that makes sense to the user", but I don't think it's worth arguing as the existing (and documented) behaviour is clear.

I will argue, though, that existing good programming practices recommend defining your kind parameters in a single place (module or include file), by means of selected_real_kind:

 integer, parameter :: sp = selected_real_kind (5)
 integer, parameter :: dp = selected_real_kind (15)

and then using "real(kind=sp)" and "real(kind=dp)" everywhere. This makes it really easy to recompile code with doubled precision by adjusting just two constants. That doesn't work, of course, for older codebases.


Hope this helps clarify the situation,
FX




Le 7 juil. 2011 à 04:33, Alan W. Irwin a écrit :

> Dear Dr. Coudert:
> 
> I noticed with considerable interest the October 2010 thread at
> http://old.nabble.com/-RFC--Quad-float-support-in-Fortran-td30071592.html.
> 
> What is the status now of that quadruple-precision gfortran effort? That thread refers to the "old" build method, but do you have a URL
> for the latest build method and latest revised patch (assuming your
> revised patch for 128-bit floating point for gfortran has not yet made
> it into an official gcc release)?  If you want to answer that question
> (and the one about the -fdefault-real-8 option below) on some list so
> it is propagated widely, please CC your answer to me in the list you
> prefer, and I will subscribe to it to follow up if there are any
> further questions.
> 
> My platform is Debian stable (on AMD64 hardware) with gfortran version
> "GNU Fortran (Debian 4.4.5-8) 4.4.5", and I do have a lot of recent
> experience building software packages (although the last time I was
> forced to build gcc was 1996 so I could get access to g77 back then).
> So I should feel reasonably comfortable building gcc (patched or new
> official version) that gives me access to a quadruple-precision
> capable gfortran.  I am probably "preaching to the choir" here, but I
> would not be so comfortable wasting a lot of my time and effort
> learning how to use a commercial Fortran compiler such as ifort
> because I far prefer to learn about and use free software tools such
> as gfortran for my scientific research.
> 
> The reason I am highly interested in your gfortran quadruple precision
> work is I have been frustrated for years in my research with
> scientific/math problems that required large least squares solutions
> that were impossible due to ill-conditioning problems.
> 
> To take but one example, there exist precise time ephemeris
> computations I have done a decade ago (using g77) based on the JPL
> planet ephemerides which stretch for many thousands of years before
> and after the year 2000. Those calculated time ephemeris data can be
> represented as a straight sloped line (the secular term that accounts
> for earth-based clocks running at a different mean rate because of the
> mean gravitational field they are in) + linear combinations of
> sinusoids with many different integer combinations of planetary
> orbital frequencies (including some frequencies that correspond to
> periods exceeding 1000 years). That secular term is of theoretical
> interest because its value affects many known physical constants and
> important timing corrections of astronomical and space craft data.
> Thus, it is essential to use as reliable a method as possible to
> separate out that secular term from the rest of my numerically
> calculated time ephemeris data.  (My paper on this subject can be
> found at http://adsabs.harvard.edu/full/1999A&A...348..642I).
> 
> My calculated time ephemeris results have good numerical precision
> (3.d-13 seconds out of a function [if the secular term is ignored]
> whose maximum deviation from zero is ~ 1.7d-3 seconds), and the
> calculations are tabulated every half day for several thousands of
> years.  Those excellent numerical qualities of the time ephemeris
> function data being fitted mean that double-precision least squares
> with singular value decomposition will probably give good results up
> to several hundreds of sinusoidal components (+ secular term) being
> fitted.  Of course, the reliability of the determined secular term
> depends on the number of sinusoids used to fit the data, but once that
> number gets too large, then the whole fit will break down due to
> ill-conditioning problems in the linear system of equations that must
> be solved.
> 
> My experience with other large least squares fits is that inevitable
> ill-conditioning breakdown tends to occur substantially before the fit
> is a good representation of the original high-quality data.
> Nevertheless, I intend to go ahead and programme that least squares
> fit in double precision to see where the breakdown occurs, and how bad
> the fit is at that point, but I would definitely like to have the
> option to switch to quadruple precision to get substantially beyond
> that double-precision ill-conditioning breakdown point.
> 
> My present version of gfortran has an option -fdefault-real-8 which
> purports to "promote the default width of 'DOUBLE PRECISION' to 16
> bytes if possible".  Currently that appears to be a no-op.  Will that
> option work correctly for your patched version?  If working, will it
> be specifically limited to just types specified by the "DOUBLE
> PRECISION" string or will all the other possibilities for specifying a
> type corresponding to 64-bit floating point (e.g., real*8, or the
> various KIND possibilities) also be promoted to 128-bit floating point
> by that option?  Uniform promotion is the only version of that option
> that makes sense to the user, but I wanted to confirm that with you.
> 
> If that option for your patched gfortran uniformly promotes all 64-bit
> floating-point varieties of constants, variables, and functions to
> 128-bit floating point, then I could use that convenient switch to
> compile the bog-standard lapack/blas double-precision routines that I
> currently use for my fits to obtain quadruple-precision results.  If
> you don't expect -fdefault-real-8 to work in that uniform way (or at
> all), I will dig further to find a non-standard quadruple precision
> lapack/blas or roll my own quadruple-precision lapack/blas version
> using a sed script to convert double-precision type statements,
> constants and math library calls, etc., in that source code to
> quadruple precision.  And similarly for my own double precision source
> code that calls the lapack/blas routines.
> 
> Alan
> __________________________
> Alan W. Irwin
> 
> Astronomical research affiliation with Department of Physics and Astronomy,
> University of Victoria (astrowww.phys.uvic.ca).
> 
> Programming affiliations with the FreeEOS equation-of-state implementation
> for stellar interiors (freeeos.sf.net); PLplot scientific plotting software
> package (plplot.org); the libLASi project (unifont.org/lasi); the Loads of
> Linux Links project (loll.sf.net); and the Linux Brochure Project
> (lbproject.sf.net).
> __________________________
> 
> Linux-powered Science
> __________________________