Paul Hinker's Weblog

Using F95 interfaces to customize access to the Sun Performance Library

Porting dusty deck Fortran source can be an exercise in patience and conditional compilation. An application which needs to run in the ILP-32, LP-64, and ILP-64 models faces the problem of interfacing with external libraries seamlessly.

Using the BLAS AXPY family of routines (caxpy, daxpy, saxpy, and zaxpy) as an example:

ILP-32 interface

SUBROUTINE saxpy(N, ALPHA, X, INCX, Y, INCY)
INTEGER*4 :: N, INCX, INCY
REAL*4    :: ALPHA, X(*), Y(*)

LP-64 interface

SUBROUTINE saxpy(N, ALPHA, X, INCX, Y, INCY)
INTEGER*4 :: N, INCX, INCY
REAL*4    :: ALPHA, X(*), Y(*)

ILP-64 interface 

SUBROUTINE saxpy(N, ALPHA, X, INCX, Y, INCY)
INTEGER*8 :: N, INCX, INCY
REAL*4    :: ALPHA, X(*), Y(*)

ILP-64 interface (strict Fortran type adherence)

SUBROUTINE saxpy(N, ALPHA, X, INCX, Y, INCY)
INTEGER*8 :: N, INCX, INCY
REAL*8    :: ALPHA, X(*), Y(*)

Strict Fortran adherence means that INTEGER and REAL data types have identical bit-width. This flies in the face of the strong implication in the BLAS documentation that single precision (i.e. the 's' prefixed routines) expect 32-bit floating point data.

The ILP-64 interfaces to the Sun Performance Library routines are suffixed by _64 to distinguish them from routines by the same name which expect 32-bit integers but run in an LP-64 model. Some older Cray source codes follow a strict adherence to the Fortran Language specification which requires INTEGER and REAL data types to have the same bit-width and so, expect the floating point data sent to the 's' prefixed routines to be 64-bit.

There are a variety of ways to handle this. First, a purely brute force approach of manually editing the source with conditional compilation code:

#ifdef ILP64
   call saxpy_64(n,alpha,incx,y,incy)
#else
   call saxpy(n,alpha,incx,y,incy)
#endif

Doing this to an application which consists of 1000's of files and millions of lines of source is a waste of engineering time. The same can be done with an awk, sed, or perl script but there are 1700+ routines in the Performance Library and even scripting the process will be time consuming and error prone.

Finally, the Fortran 95 generic interface functionality could be used to allow source to remain virtually unchanged and yet facilitate the use of each of the three programming models (ILP-32, LP-64, and ILP-64).

Let's take a very simple example which calls the BLAS saxpy routine:

% cat tst.f

      program tst saxpy
      implicit none

      integer, parameter :: N = 10, INCX = 1, INCY = 1
      real, parameter    :: ALPHA = 1.0
      real, dimension(N) :: X, Y

      X = 1.0 
      Y = 2.0
   
      call saxpy(N, ALPHA, X, INCX, Y, INCY)
      print *, SUM(Y)
      END

% f90 -o tst tst.f -xlic_lib=sunperf -xarch=[v8|v8plus|v8plusb|sse2]

When compiled with one of the ILP-32 architectures (v8, v8plus, v8plusb, sse2) the saxpy call resolves to the one expecting 32-bit integer and real parameters.

% f90 -o tst tst.f -xlic_lib=sunperf -xarch=[v9|v9b|amd64]

When compiled with one of the ILP-64 architectures (v9, v9b, amd64) the saxpy call resolves to the LP-64 interface which expects 32-bit integer and real parameters. This is for backward compatibility reasons. However, when compiled with one of the ILP-64 libraries, additional entry points are available. That is, our source could look like the following :

      program tst saxpy
      implicit none

      integer, parameter   :: N = 10, INCX = 1, INCY = 1
      real, parameter      :: ALPHA = 1.0
      real, dimension(N)   :: X, Y

      X = 1.0 
      Y = 2.0
   
      call saxpy_64(N, ALPHA, X, INCX, Y, INCY)
      print *, SUM(Y)
      END

% f90 -o tst tst.f -xlic_lib=sunperf -xarch=[v9|v9b|amd64] -xtypemap=integer:64

Note : The integer declarations have been changed to 8 byte integers using the -xtypemap compiler option.

Warning : The xtypemap option only applies to implicitly typed variables. INTEGER :: N gets promoted to INTEGER*8 :: N but INTEGER*4 :: N would not be promoted.

In theory, the xtypemap option is a great way to port to the ILP-64 model. In practice, it's not quite a silver bullet. As long as interfaces are clearly defined and strictly follow typing rules it works well.

A F95 interface can be used to describe the different programming models.

      module sunperf64
         interface saxpy
!
! ILP-32 and LP-64 interface
!
            subroutine saxpy(n,alpha,x,incx,y,incy)
               integer :: n, incx, incy
               real    :: alpha, x(*), y(*)
            end subroutine
!
! ILP-64 interface
!
            subroutine saxpy_64(n,alpha,x,incx,y,incy)
               integer * 8 :: n, incx, incy
               real        :: alpha, x(*), y(*)
            end subroutine
!
! ILP-64 interface w/strict Fortran typing
!
            subroutine daxpy_64(n,alpha,x,incx,y,incy)
               integer * 8 :: n, incx, incy
               real * 8    :: alpha, x(*), y(*)
            end subroutine
         end interface
      end module

This module file will allow a single call to the saxpy routine to be interpreted as any of the ILP-32, LP-64, ILP-64, ILP-64(strict) calling conventions.

      program tst saxpy
      implicit none
      use sunperf64

      integer, parameter   :: N = 10, INCX = 1, INCY = 1
      real, parameter      :: ALPHA = 1.0
      real, dimension(N)   :: X, Y

      X = 1.0 
      Y = 2.0
   
      call saxpy(N, ALPHA, X, INCX, Y, INCY)
      print *, SUM(Y)
      END

% f90 -o tst tst.f -xlic_lib=sunperf -xarch=[v8plus|v8plusb|sse2]
Would call the ILP-32 version 
% f90 -o tst tst.f -xlic_lib=sunperf -xarch=[v9|v9b|amd64] -xtypemap=integer:64
Would call ILP-64 version
% f90 -o tst tst.f -xlic_lib=sunperf -xarch=[v9|v9b|amd64] -xtypemap=integer:64,real:64
Would call the ILP-64(strict) version

Posted at 12:41PM Nov 13, 2006 by hinkthink in General  |  Comments[1]