Friday, 30 July 2010

Hey, you sass that hoopy Fortran?

We have to work with lots of different Fortran compilers at NAG. So far at Mark 22 of our Fortran Library we have built using products from six distinct sources (GNU, IBM, Intel®, NAG, PGI and Sun). The list of what we build with naturally changes from Mark to Mark as different vendors fall by the wayside or as commercial considerations preclude us from creating a particular Library implementation.

With the six vendors listed above we are blessed that the base Fortran coverage is Fortran 95. At previous Marks we still supported compilers covering only Fortran 77 (g77 and pgf77 for example). This forced us to jump through all manner of hoops so that we could do some of the really useful stuff from Fortran 90—like, woah, dynamically allocating memory—in a Fortran 77 way.

Some hoops still exist though, even with Fortran 90 code.


For example, we recently ran into the following:
    PROGRAM str_read
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: wp = KIND(0.0D0)
!      .. Local Scalars ..
       REAL (kind=wp)                  :: rval
       INTEGER                         :: ioerr
       CHARACTER (200)                 :: my_str
!      .. Executable Statements ..
       my_str = '1.0D400'
       READ (my_str,'(E16.0)',iostat=ioerr) rval
       PRINT *, 'IOERR = ', ioerr
    END PROGRAM str_read
Ignoring the fact that the code fragment is a little contextless, just imagine that we're trying to trap overflowing real values being read from a CHARACTER. On five of the six compilers given above, two return a nonzero ioerr from the READ, two carry on happily and set rval to Infinity, and one core dumps! (At the time of writing, I didn't have access to all six compilers.)

Conferring with the Fortran standard (e.g., the draft Fortran 95) we see the dreaded phrase The set of input/output error conditions is processor dependent. (my italics), so we can't even complain to the compiler vendors about this behaviour! We have good relationships with the vendors, so when our building process reveals a compiler bug we'll report it. But as you can imagine, sometimes it can be a battle to achieve what we want across many Fortran platforms in as simple and maintainable a way as possible.

The next level of complication comes from ensuring that what we do in Fortran is callable from non-Fortran environments. We usually take Microsoft Excel as a sufficiently-far-removed example. I won't go into details of the usual issues raised by cross-language programming; however, now we're at a base level of Fortran 95 we've been looking at pushing our envelope more and trying out some Fortran 2003 features, specifically C Interoperability.

M'colleague Nicolas has been working on a suite of image processing routines. In pure Fortran a NAG_IMAGE TYPE exists to hold the pixel values and other metadata for the image. To make this passable from outside Fortran we're using Fortran 2003. It's understandable with new features in the language that there should be a settling-in period, but it's no exaggeration to say that every compiler we tried did, at some point, fall over on our new code. Thanks to the hard work of the compiler developers on the various teams we're now at a stable state, but there are still some vendors who are lagging behind with Fortran 2003 and for which we have to exclude the image-processing code from our testing.

Looking forward (!) to Fortran 2008, I'm interested to see how the new special-function INTRINSICs (BESSEL_J*, BESSEL_Y*, ERF* and *GAMMA) will behave across different compilers. In theory these new INTRINSICs mean we'd be able to withdraw a sizeable chunk of our S Chapter. But I wonder how easy that will be?


(The images in this post are from now voyager, but are probably copyright PolyGram Filmed Entertainment.)

5 comments:

  1. I think those S chapter functions will be around for a while yet. Remember that the library is called from many non-Fortran environments such as C.

    Cheers,
    Mike

    ReplyDelete
  2. The TR1 extensions to the C++ language (see http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1836.pdf for a draft version) include some challenging special functions such as the hypergeometric functions 2F1 and 1F1 and elliptic integrals.

    Unfortunately there are some areas not specified in detail in the (draft) standard. For example, the 2F1 function (5.2.1.17) for real arguments is required to return a value given by the standard Gauss series (DLMF 15.2.1), but TR1 does not mention that this series is only valid for |z|<1 and the function is defined elsewhere by analytic continuation. In addition there is no information about the numerical accuracy to be required of the implementation!

    So I suspect the S chapter in both C and FORTRAN will have a considerable life and I look forward to finding 2F1 as a member in a future NAG release!

    Alan Bain

    ReplyDelete
  3. Thanks Alan, that's interesting news. I didn't know there was a drive to try to standardize this kind of functionality outside of Fortran. I like how both standards seem to be opting for the 'implementation-dependent approximation' approach (ho-hum)!

    ReplyDelete
  4. I want to print an variable size array as a series of lines, each with a set number of elements, until the last line. The last has a variable number of elements, equal to or less than the previous lines.

    Consider the following array:
    a(1)=1, a(2)=2, a(3)=3, a(4)=4, a(5)=5

    The print appears as:
    1 2 3
    4 5

    Consider the another array:

    a(1)=1, a(2)=2, a(3)=3, a(4)=4, a(5)=5 a(6)=6, a(7)=7

    The print appears as:
    1 2 3
    4 5 6
    7

    Can this be done without programming loops and indices, using Fortran statements?

    Or is it necessary to write special code, say a callable subroutine, that accepts the array, and its length then does the required line printing allocation?

    Thanks.

    ReplyDelete
  5. leomane1 - Fortran formatted output is pretty
    powerful if a little arcane. Here's how to do what you want with integer arrays. Note that
    the string '(3i2)' in the write statement means
    "print up to three integers with field width 2".
    If there are more than three items in the
    array, the last part of the format gets reused
    over and over again. For more information see
    the section on formatted about in any good
    Fortran book.

    integer a(5), b(7)

    a = (/ 1, 2, 3, 4, 5 /)
    b = (/ 1, 2, 3, 4, 5, 6, 7 /)

    write (*,'(3i2)') a
    write (*,*)
    write (*,'(3i2)') b
    end

    Results of running this program:

    1 2 3
    4 5

    1 2 3
    4 5 6
    7

    ReplyDelete

NAG moderates all replies and reserves the right to not publish posts that are deemed inappropriate.