Friday, 30 July 2010

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.

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_IMAGETYPE 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.)

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!

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)!

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 34 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 34 5 67

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?

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

About NAG

The Numerical Algorithms Group (NAG) is dedicated to applying its unique expertise in numerical engineering to delivering high-quality computational software and high performance computing services. For 40 years NAG experts have worked closely with world-leading researchers in academia and industry to create powerful, reliable and flexible software which today is relied on by tens of thousands of individual users, as well as numerous independent software vendors. NAG serves its customers from offices in Oxford, Manchester, Chicago, Tokyo and Taipei, through local staff in France and Germany, as well as via a global network of distributors.