ARSC HPC Users' Newsletter 263, February 7, 2003

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.

The meetings will also be available via the Access Grid:

For more information contact:

Paul Mercer of ARSC.

Email: , ph: 907-474-3785.

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:

  program namelist_tester
    implicit none
    integer               :: i = -1
    character (len=20)    :: text = ""
    logical               :: L = .false.
    real,dimension(10)    :: x = 10*(-1)
    namelist /nml_a/ i, text
    namelist /nml_b/ L
    namelist /nml_c/ x
    open (77, file="test.nml")
    read (77, nml=nml_a)
    close (77)
    print*, "nml_a: ", i, text
    open (77, file="test.nml")
    read (77, nml=nml_c)
    close (77)
    print*, "nml_c: ", x

Here's a file, "test.nml," which can serve as input to the program. It provides data for all three groups, even though "nml_b" is unused:

    i = 77, 
    text = "Some Text" /
      L = .true.  /
      x=1000.1,1002.1,1003.1,6*1004.1,999.9 /

Here's a run on an IBM SP:

  ICEHAWK2$ xlf90 -qsuffix=f=f90 -o test test.f90      
  ** namelist_tester   === End of Compilation 1 ===
  1501-510  Compilation successful for file test.f90.
  ICEHAWK2$ ./test
   nml_a:  77 Some Text           
   nml_c:  1000.099976 1002.099976 1003.099976 1004.099976 1004.099976
           1004.099976 1004.099976 1004.099976 1004.099976 999.9000244

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'.
   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.

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.

And a test...

  CHILKOOT$ eval `assign -Y on u:77`
  CHILKOOT$ ./test
   nml_a:  77 Some Text           
   nml_c:  1000.099999999999,  1002.099999999999,  1003.099999999999,
           6*1004.099999999999,  999.9000000000015

An additional problem is possible if the namelists are read out of order. The test program, above, avoids this problem completely by opening and closing the file for each READ. Like this:

    open (77, file="test.nml")
    read (77, nml=nml_a)
    close (77)

    open (77, file="test.nml")
    read (77, nml=nml_c)
    close (77)

Since the READs are in order, it would also work with just one open:

    open (77, file="test.nml")
    read (77, nml=nml_a)
    read (77, nml=nml_c)
    close (77)

The following would fail on any system, however, because the file position indicator would be at the end of the file after reading nml_c, and thus, scanning forward would not find nml_b:

    open (77, file="test.nml")
    read (77, nml=nml_a)
    read (77, nml=nml_c)
    read (77, nml=nml_b)
    close (77)

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:

      integer, parameter :: r8  = selected_real_kind(P1,R1)
      integer, parameter :: r16 = selected_real_kind(P2,R2)

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:

      PROGRAM ranger

      real*4  :: b04
      real*8  :: b08
      real*16 :: b16

      print *, 'real*4  ', precision(b04), range(b04)
      print *, 'real*8  ', precision(b08), range(b08)
      print *, 'real*16 ', precision(b16), range(b16)

      END PROGRAM ranger

We can get some surprisingly different values:

 real*4   6,  37
 real*8   15,  307
 real*16  31,  275

 real*4   6 37
 real*8   15 307
 real*16  31 291

 real*4   6 37
 real*8   15 307
 real*16  33 4931

Cray T3E:
f90-391 f90: WARNING RANGER, File = range.f90, Line = 5, Column = 12 
  Type REAL*16 will be mapped to REAL*8.

 real*4   6,  37
 real*8   15,  307
 real*16  15,  307

Cray SV1:
 real*4   6,  2465
 real*8   13,  2465
 real*16  28,  2465

Cray SX-6:
 real*4    6  37
 real*8    15  307
 real*16   31  307

Linux Portland Group Fortran:
PGF90-W-0031-Illegal data type length specifier for real (range.f90: 5)
PGF90-W-0031-Illegal data type length specifier for b16 (range.f90: 5)

 real*4              6           37
 real*8             15          307
 real*16             6           37

So what does this mean for the problem at hand? We need to use precision and range values at least as small as the smallest of our results here:

      integer, parameter :: r4  = selected_real_kind(6,37)
      integer, parameter :: r8  = selected_real_kind(13,307)
      integer, parameter :: r16 = selected_real_kind(28,275)

We will get an error when using r16 on the T3E and Linux box, but it should work on everything else. The return value of selected_real_kind on an error is:

   -1 if precision unavailable
   -2 if range unavailable
   -3 if neither is available

We should make sure that these values work as expected. Let's test it by printing out pi at each precision. Here is the test code:

      PROGRAM pi

      integer, parameter :: r4  = selected_real_kind(6,37)
      integer, parameter :: r8  = selected_real_kind(13,307)
      integer, parameter :: r16 = selected_real_kind(28,275)

      real(r4)  :: a04
      real(r8)  :: a08
      real(r16) :: a16

      a04 = 4.0_r4*atan(1.0_r4)
      a08 = 4.0_r8*atan(1.0_r8)
      a16 = 4.0_r16*atan(1.0_r16)

      print *, a04
      print *, a08
      print *, a16

      END PROGRAM pi


  Cray SV1:
  Cray T3E:
   -1                         result of selected_real_kind(28,275)
  Cray SX-6:
             -1               result of selected_real_kind(28,275)

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:

f90-130 f90: ERROR PI, File = pi.f90, Line = 9, Column = 12 
  The kind type parameter value -1 is not valid for type REAL.

Any readers with a suggestion other than using the C preprocessor?

============== Editor's Note: ==============

Similar functions exist for integer variables:

  RANGE ()

The big surprise is that 8-byte integers on the SX-6 have a smaller range than is possible from 64 bits.

From the SX-6 manual ("float0" format only):

  Values that are represented are integers within the range shown

    2-byte integer
      Maximum 32767
      Minimum -32768

    4-byte integer
      Maximum 2147483647
      Minimum -2147483648

    8-byte integer
      Maximum 2^52 -1
      Minimum -2^52 +1


        When only add and subtract instructions are
        related, 8-byte integer type is as follows. 
          Maximum: 2^63 -1
          Minimum: -2^63 

Minor changes to Kate's program:

      PROGRAM rangei

      integer*4  :: b04
      integer*8  :: b08

      print *, 'integer*4  ', range(b04)
      print *, 'integer*8  ', range(b08)

      END PROGRAM ranger

Gives these results:

   integer*4    9
   integer*8    15
   integer*4   9
   integer*8   18
   integer*4   9
   integer*8   18
   integer*4   9
   integer*8   18
   integer*4   9
   integer*8   18

Training Opportunities

All local ARSC users are invited to participate in the following classes:

  User's Introduction to ARSC Supercomputers,
    Date:       Weds., Feb 12th, 2pm
    Location:   Butrovich 109
    Instructor: Kate Hedstrom, ARSC
  Introduction to the ARSC IBM Regatta System,
    Date:       Weds., Feb 19th, 2pm
    Location:   Butrovich 007
    Instructor: Jim Long, ARSC
  Introduction to ARSC Visualization Software and Facilities
    Date:       Weds., March 19th, 2pm
    Location:   Butrovich 109/007
    Instructor: TBD

For registration and details, see:

Also available:

  Intermediate MPI Programming
    February 6, 13, 20, and 27, 2003 (Every Thursday in February)
    9am - 12noon, Alaska Time, each day.
    Butrovich 109, on the UAF campus 

  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:

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

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
    #define M_PI        3.14159265358979323846264338327950288
    #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) : 

     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?

[[ Answers, Questions, and Tips Graciously Accepted ]]

Current Editors:
Ed Kornkven ARSC HPC Specialist ph: 907-450-8669
Kate Hedstrom ARSC Oceanographic Specialist ph: 907-450-8678
Arctic Region Supercomputing Center
University of Alaska Fairbanks
PO Box 756020
Fairbanks AK 99775-6020
E-mail Subscriptions: Archives:
    Back issues of the ASCII e-mail edition of the ARSC T3D/T3E/HPC Users' Newsletter are available by request. Please contact the editors.

Back to Top