ARSC HPC Users' Newsletter 263, February 7, 2003

Contents

UAF Bioinformatics Journal Club

If you're trying to learn a little more from the world of bioinformatics, please join informal, bi-monthly discussions of relevant topics and journal articles, to be lead by Nat Goodman.

The first meeting will be on Thursday, Feb 20, 12-1:30 Alaska Time, Butrovich 109. The focus, for at least the opening few months, will be on classic techniques and applications. For the first session, the topic will be the Big Three of Sequence Analysis: Smith-Waterman, Fasta, and BLAST.

Fortran 90 Namelist Behavior

A user had difficulties last week porting an f90 code to the Cray SV1ex. The problem was in reading a namelist from a file containing several. ("Namelist I/O" is a handy Fortran facility for reading/writing groups of variables, annotated by name. It's often used for initialization.)

On the Crays, the default behavior for a READ statement is to not scan ahead for a missing namelist. To enable such scanning, the user had to apply "assign -Y on" to the file. This reproduced the behavior which is default on the SGI's, IBMs, and SX-6.

Here's the unavoidable "small test program." It defines three namelist groups, nml_a, nml_b, and nml_c, but only reads two of them, nml_a and nml_c:

On the SV1ex, it reads "nml_a", but crashes on the attempt to read "nml_c":

CHILKOOT$ f90 -o test test.f90
CHILKOOT$ ./test
nml_a: 77 Some Text
lib-1304 ./test: UNRECOVERABLE library error
Namelist input group name "NML_A" does not
match READ group name.
Encountered during a namelist READ from unit 77
Fortran unit 77 is connected to a sequential formatted text file: "test.nml"
Error initiated at line 403 in routine '_FRN'.
Abort
Beginning of Traceback:
Interrupt at address 67347a in routine '_lwp_killm'.
Called from line 32 (address 64036a) in routine 'raise'.
Called from line 127 (address 1106d) in routine 'abort'.
Called from line 59 (address 371247a) in routine '_ferr'.
Called from line 403 (address 435246d) in routine '_FRN'.
Called from line 21 (address 613d) in routine 'NAMELIST_TESTER'.
Called from line 350 (address 22715c) in routine '$START$'.
End of Traceback.
Abort(coredump)

The solution is the assign statement. From "man assign":

-Y setting Skip unmatched namelist group in the namelist input
record. setting can be either on or off. The default
setting on UNICOS and UNICOS/mk systems is off. The
default setting on IRIX systems is on.

While the mistake is obvious in this example, it might be obscured in a real code which could have READ statements scattered about or even in multiple subroutines.

Getting the Right Type

[ Thanks to Kate Hedstrom of ARSC for this article. See her related article, "The Size of IBM XLF Fortran Reals," in:
issue 257
As well as: "SELECTED_REAL_KIND, Caution for Portability," also in issue 257.]

In the old days, we would write our Fortran code with REAL for the floating point, getting 64-bit values on the Cray. We have talked about the various IBM options for promoting REAL in another article. This time we want to talk about using Fortran 90 kinds to specify the real type used. Let's say that we want most of our calculations to happen in 64-bit precision, but there is one delicate operation we want in 128-bits if available.

The recommended technique is to use selected_real_kind to create the type variables:

The two arguments to selected_real_kind are the precision (number of digits) and range (largest and smallest exponent). What values should we use for these arguments? After some wrong initial guesses, we decided that it would be prudent to write a little test script to find out what types are available:

It seems to be working! The mistakes we had made previously were of two sorts: accidentally getting 128-bit reals on the SV1 by asking for a 64-bit precision of 15, and getting errors on SGI and IBM by asking for a 128-bit range of 300. We have fixed those problems.

The remaining question is what to do about systems that don't support 128-bit reals. If you ask for real*16 in the old way, you get a default real (could be 32-bit!) and a compiler warning. With the new style, you get a compiler error such as:

This course is being delivered over the access grid. It is taught by
Dr. David Ennis of the Ohio Supercomputing Center.
The first meeting was yesterday, but the materials are available if
you'd like to join for the remaining three sessions. Send inquires
to:
training@arsc.edu

Quick-Tip Q & A

A:[[ What's an easier way to get the value of pi into my C/C++ or Fortran
[[ programs?
#
# Thanks to six readers for responding. Here's a sample:
#
A standard way of getting Pi into any program, getting full machine
precision, is
Pi = 4*arctan(1.0)
We note that tan(Pi/4) = 1, which implies that arctan(1) = Pi/4,
giving us the above formula. Since the arctan function (typically
something like atan()) is written to full machine precision, you can
expect your resulting Pi to have full machine precision.
---
Easier than what???!!!!!
In Fortran, I use
pi = 4.0*atan(1.)
---
For C and C++, you can use pi=4.0*atan(1.0). Remember to include
math.h and compile with -lm.
---
Easier than what? In C/C++, I use:
#include <math.h>
After which point the name M_PI will refer to the value of pi (I
assume to double floating point precision, though I've never
checked...).
Editor's note:
==============
Inspired by the reader's comment, we thought we'd check. So, here's
M_PI from /usr/include/math.h on several systems. If you're using 128
bit precision, looks like you'd be safer to use the "atan" method:
SGI, Cray, SX-6, Linux
#define M_PI 3.14159265358979323846
IBM:
#define M_PI 3.14159265358979323846264338327950288
Linux
#define M_PIl 3.1415926535897932384626433832795029L
Q: A C style question. Don't worry what any of this does, the goal is
just to compare a couple things and choose one:
Option 1:
---------
vel = (envel / scaleFct) / ((*outv)[n] + posStrt) < minvel ?
(envel / scaleFct) / ((*outv)[n] + posStrt) :
minvel;
Option 2:
---------
vel = (envel / scaleFct) / ((*outv)[n] + posStrt);
if ( ! (vel < minvel) )
vel = minvel;
Isn't there a cleaner way? If not, which should I prefer?

The University of Alaska Fairbanks is an affirmative action/equal
opportunity employer and educational institution and is a part of the University
of Alaska system.
Arctic Region Supercomputing Center (ARSC) |PO Box 756020, Fairbanks, AK 99775 | voice: 907-450-8602 | fax: 907-450-8601 | Supporting high performance computational research in science and engineering with emphasis on high latitudes and the arctic.
For questions or comments regarding this website, contact info@arsc.edu