HP OpenVMS Systems Documentation

HP Fortran for OpenVMS
User Manual

This chapter describes:

• Section 15.1, What Is CXML?
• Section 15.2, CXML Routine Groups
• Section 15.3, Using CXML from Fortran
• Section 15.4, CXML Program Example

The Compaq Extended Math Library (CXML) provides a comprehensive set of mathematical library routines callable from Fortran and other languages. CXML contains a set of over 1500 high-performance mathematical subprograms designed for use in many different types of scientific and engineering applications. It significantly improves the run-time performance of certain HP Fortran programs.

CXML is included with HP Fortran for OpenVMS Alpha Systems and can be installed using the instructions in the HP Fortran Installation Guide for OpenVMS Alpha Systems.

The CXML reference guide is available in both online and hardcopy formats. You can obtain this documentation by accessing the following files:

• SYS$HELP:CXMLREF-VMS.PDF (online)---view with the Adobe Acrobat Reader
• SYS$HELP:CXMLREF-VMS.PS (hardcopy)---print to a PostScript printer

Example programs are also provided with CXML. These programs are located in the following directory:

SYS$COMMON:[SYSHLP.EXAMPLES.CXML]

15.2 CXML Routine Groups

CXML provides a comprehensive set of highly efficient mathematical subroutines designed for use in many different types of scientific and engineering applications. CXML includes the following functional groups of subroutines:

• Basic linear algebra
• Linear system and Eigenproblem solvers
• Sparse linear system solvers
• Signal processing
• Utility subprograms

The routines are described in Table 15-1.

Where appropriate, each subprogram has a version to support each combination of real or complex and single- or double-precision arithmetic. In addition, selected key CXML routines are available in parallel form as well as serial form on HP OpenVMS Alpha systems.

15.3 Using CXML from Fortran

To use CXML, you need to make the CXML routines and their interfaces available to your program and specify the appropriate libraries when linking.

The CXML routines can be called explicitly by your program. There are separate CXML libraries for the IEEE and the VAX floating-point formats. You must compile your program for one of these float formats and then link to the matching CXML library (either IEEE or VAX), depending upon how you compiled the program.

Either the IEEE or VAX CXML library can be established as the systemwide default by the system startup procedure. Individual users can select between the VAX and IEEE version by executing the SYS$LIBRARY:CXML$SET_LIB command procedure. For example, the following command alters the default CXML link library for the current user to the VAX format library:

$ @SYS$LIBRARY:CXML$SET_LIB VAX

For more details, see the section about CXML post-installation startup options in the HP Fortran Installation Guide for OpenVMS Alpha Systems.

If needed, you can instead specify the appropriate CXML library or libraries on the LINK command line (use the /LIBRARY qualifier after each library file name). You must compile your program and then link to the appropriate CXML library (either IEEE or VAX), depending upon how you compiled the program. The following examples show the corresponding CXML commands for compiling and linking for cases where the CXML default library is not used:

$ FORTRAN /FLOAT=IEEE_FLOAT MYPROG.F90 
$ LINK MYPROG.OBJ, SYS$LIBRARY:CXML$IMAGELIB_TS/LIBRARY

The link command uses the name of the CXML library for IEEE floating-point data, cxml$imagelib_ts . To use VAX floating-point data, specify the CXML library name as cxml$imagelib_gs .

If you are using an older version of CXML, use dxml$xxxxx instead of cxml$xxxxx as the library name. For more information on using CXML and specifying the correct object libraries on the LINK command, see the Compaq Extended Mathematical Library Reference Manual.

15.4 CXML Program Example

The free-form Fortran 90 example program below invokes the function SAXPY from the BLAS portion of the CXML libraries. The SAXPY function computes a*x+y .

   PROGRAM example
!
!  This free-form example demonstrates how to call
!  CXML routines from Fortran.
!
   REAL(KIND=4)    :: a(10)
   REAL(KIND=4)    :: x(10)
   REAL(KIND=4)    :: alpha
   INTEGER(KIND=4) :: n
   INTEGER(KIND=4) :: incx
   INTEGER(KIND=4) :: incy
   n = 5 ; incx = 1 ; incy = 1 ; alpha = 3.0
   DO i = 1,n
     a(i) = FLOAT(i)
     x(i) = FLOAT(2*i)
   ENDDO
   PRINT 98, (a(i),i=1,n)
   PRINT 98, (x(i),i=1,n)
98 FORMAT(' Input = ',10F7.3)
   CALL saxpy( n, alpha, a, incx, x, incy )
   PRINT 99, (x(i),I=1,n)
99 FORMAT(/,' Result = ',10F7.3)
   STOP
   END PROGRAM example

This appendix describes:

• Section A.1,HP Fortran Commands on OpenVMS I64 That Are Not Available on OpenVMS Alpha
• Section A.2,HP Fortran Commands on OpenVMS Alpha That Are Not Available on OpenVMS I64
• Section A.3, Differences in Default Values
• Section A.4, Support for VAX-Format Floating-Point
• Section A.5, Changes in Exception Numbers and Places
• Section A.6, Changes in Exception-Mode Selection

• /CHECK=ARG_INFO
• /CHECK=FP_MODE

• /ARCHITECTURE
• /MATH_LIBRARY
• /OLD_F77

• /OPTIMIZE=TUNE
• /SYNCHRONOUS_EXCEPTIONS

The differences are:

• Because the Itanium architecture supports IEEE directly, the default floating-point format on OpenVMS I64 systems is /FLOAT=IEEE_FLOAT. On OpenVMS Alpha systems, the default is /FLOAT=G_FLOAT. See Section 2.3.22.
• The default floating-point exception handling mode on OpenVMS I64 systems is /IEEE_MODE=DENORM_RESULTS). On OpenVMS Alpha systems, the default is /IEEE_MODE=FAST. See Section 2.3.24.

Because there is no direct hardware support for VAX-format floating-point on OpenVMS I64 systems, the VAX-format floating-point formats (F, D, and G) are supported indirectly by a three-step process:

1. Conversion from VAX-format variable to IEEE floating-point temporary (using the denormalized number range)
2. Computation in IEEE floating-point
3. Conversion back to VAX-format floating-point and storage into the target variable

There are a number of implications for this approach:

• Exceptions might move, appear, or disappear, because the calculation is done in a different format with a different range of representable values.
• Because there are very small numbers representable in VAX format that can only be represented in IEEE format using the IEEE denorm range, there is a potential loss of accuracy if an intermediate result of the calculation is one of those numbers.
  At worst, the two bits with the least significance will be set to zero. Note that further calculations might result in magnifying the magnitude of this loss of accuracy.
  Note that this small loss of accuracy does not raise signal FOR$_SIGLOSMAT (or FOR$IOS_SIGLOSMAT).
• Expressions that are used to drive control flow but are not stored back into a variable will not be converted back into VAX format. This can cause exceptions to disappear.
• There can be a significant performance cost for the use of VAX-format floating-point.

Note that floating-point query built-ins (such as TINY and HUGE) will return values appropriate to the floating-point format that you select, despite the fact that all formats are supported by IEEE.

A.5 Changes in Exception Numbers and Places

There will be changes in the number of exceptions raised and in the code location at which they are raised.

This is particularly true for VAX-format floating-point calculations, because many exceptions will only be raised at the point where a result is converted from IEEE format to VAX format. Some valid IEEE-format numbers will be too large or too small to convert and will thus raise underflow or overflow. IEEE exceptional values (such as Infinity and NaN) produced during the evaluation of an expression will not generate exceptions until the final conversion to VAX format is done.

If a VAX-format floating-point calculation has intermediate results (such as the X * Y in the expression (X * Y)/ Z ), and the calculation of that intermediate result raised an exception on OpenVMS Alpha systems, it is not guaranteed that an exception will be raised on OpenVMS I64 systems. An exception will only be raised if the IEEE calculation produces an exception.

A.5.1 Ranges of Representable Values

In general, the range of VAX-format floating-point numbers is the same as the range for IEEE-format. However, the smallest F- or G-format value is one quarter the size of the smallest normal IEEE number, while the largest F- or G-format number is about half that of the largest IEEE number. There are therefore nonexceptional IEEE values that would raise overflows in F- or G-format. There are also nonexceptional F- or G-format values that would raise underflow in IEEE-format in those modes in which denormalized numbers are not supported.

A.5.2 Underflow in VAX Format with /CHECK=UNDERFLOW

OpenVMS Alpha and VAX Fortran applications do not report underflows for VAX-format floating-point operations unless you specifically enable underflow traps by compiling with the /CHECK=UNDERFLOW qualifier (see Section 2.3.11).

The same is true on OpenVMS I64 systems, but with an important caveat: Since all I64 floating-point operations are implemented by means of IEEE-format operations, enabling underflow traps with /CHECK=UNDERFLOW causes exceptions to be raised when values underflow the IEEE-format representation, not the VAX-format one.

This can result in an increased number of underflow exceptions seen with /CHECK=UNDERFLOW when compared with equivalent Alpha or VAX programs, as the computed values may be in the valid VAX-format range, but in the denormalized IEEE-format range.

If your application requires it, a user-written exception handler could catch the IEEE-format underflow exception, inspect the actual value, and determine whether it represented a VAX-format underflow or not.

See Section 8.4 for exact ranges of VAX-format and IEEE-format floating point.

A.6 Changes in Exception-Mode Selection

On OpenVMS Alpha systems, the exception-handling mode and the rounding mode can be chosen on a per-routine basis. This lets you set a desired exception mode and rounding mode using compiler qualifiers. Thus, a single application can have different modes during the execution of different routines.

This is not as easy to do on OpenVMS I64 systems. While the modes can be changed during the execution of a program, there is a significant performance penalty for doing so.

As a result, the HP Fortran compiler and the OpenVMS linker implement a "whole-program" rule for exception handling and rounding modes. The rule says that the whole program is expected to run in the same mode, and that all compilations will have been done using the same mode qualifiers. To assist in enforcing this rule, the compiler, linker and loader work together:

• The compiler tags each compiled object file with a code specifying the modes selected by the user (directly or using the default) for that compilation.
• The linker tags the generated executable file with a code specifying the mode selected by the user.
• The loader initializes the floating-point status register of the process based on the linker code.

If you are using an OpenVMS I64 system and want to change the exception-handling or rounding mode during the execution of a program, use a call to either of the following:

• Fortran routines DFOR$GET_FPE and DFOR$SET_FPE
• System routines SYS$IEEE_SET_FP_CONTROL, SYS$IEEE_SET_PRECISION_MODE, and SYS$IEEE_SET_ROUNDING_MODE

If you call DFOR$GET_FPE and DFOR$SET_FPE, you need to construct a mask using the literals in SYS$LIBRARY:FORDEF.FOR, module FOR_FPE_FLAGS.

The calling format is:

INTEGER*4  OLD_FLAGS, NEW_FLAGS
INTEGER*4  DFOR$GET_FPE, DFOR$GET_FPE

EXTERNAL   DFOR$GET_FPE, DFOR$GET_FPE

! Set the new mask, return old mask.
OLD_FLAGS = DFOR$GET_FPE( NEW_FLAGS )

! Return the current FPE mask.
OLD_FLAGS = DFOR$GET_FPE ()

An example (which does no actual computations) follows. For a more complete example, see Example A-1.

subroutine via_fortran
include 'sys$library:fordef.for'
include '($IEEEDEF)'
integer new_flags, old_flags
old_flags = dfor$get_fpe()
new_flags = FOR_K_FPE_CNT_UNDERFLOW + FPE_M_TRAP_UND
call dfor$set_fpe( new_flags )
!
!  Code here uses new flags
!
call dfor$set_fpe( old_flags )
return
end subroutine

If you call SYS$IEEE_SET_FP_CONTROL, SYS$IEEE_SET_PRECISION_MODE, and SYS$IEEE_SET_ROUNDING_MODE, you need to construct a set of masks using the literals in SYS$LIBRARY:FORSYSDEF.TLB, defined in module IEEEDEF.H (in STARLET). For information about the signature of these routines, see the HP OpenVMS System Services Reference Manual.

An example (which does no actual computations) follows. For a more complete example, see Example A-1.

subroutine via_system( rr )
real rr
include '($IEEEDEF)'
integer*8 clear_flags, set_flags, old_flags, new_flags
clear_flags = IEEE$M_MAP_DNZ + IEEE$M_MAP_UMZ
set_flags   = IEEE$M_TRAP_ENABLE_UNF + IEEE$M_TRAP_ENABLE_DNOE
call sys$ieee_set_fp_control(%ref(clear_flags),
%ref(set_flags),%ref(old_flags))
!
!  Code here uses new flags
!
clear_flags = set_flags
call sys$ieee_set_fp_control(%ref(clear_flags),%ref(old_flags),%ref(new_flags))
return

If you decide to change the exception-handling or rounding mode, be careful to observe the following rules to maintain the "whole-program" rule. Failure to do so might cause unexpected errors in other parts of your program:

• The preexisting mode must be restored at the end of the execution of the section in which the new mode is used. This includes both normal endings, such as leaving a code block, and exits by means of exception handlers.
• It is a good idea to establish a handler to restore the old mode on unwinds, because called code can cause exceptions to be raised (including exceptions not related to floating point).
• The code should be compiled with the same mode qualifiers as the other, "normal" parts of the application, not with the mode that will be set by the call to the special function.
• Be aware that VAX-format expressions are actually calculated in IEEE format, and any change to the modes will impact a calculation in IEEE format, not a calculation in VAX format.
• Consider adding the VOLATILE attribute to the declaration of all the variables used in the calculation in the different mode. This will prevent optimizations that might move all or part of the calculation out of the region in which the different mode is in effect.

System libraries that need to use an alternate mode (for example, the math library) accomplish this by using an architectural mechanism not available to user code: the .sf1 flags of the floating-point status register (user code uses the .sf0 flags).

Therefore, a user's choice of exception-handling or rounding mode will not have an impact on any system library used by the program.

A.6.2 Example of Changing Floating-Point Exception-Handling Mode

Example A-1 shows both methods of changing the floating-point exception-handling mode. However, for a real program, you should pick just one of the two methods.

!  SET_FPE.F90: Change floating-point exception handling mode,
!               and check that it worked.
!
!  Compile and link like this:
!
!  $ f90 set_fpe
!  $ lin set_fpe,SYS$LIBRARY:VMS$VOLATILE_PRIVATE_INTERFACES.OLB/lib
!  $ run set_fpe
!
!  The library is needed to bring in the code for LIB$I64_INS_DECODE,
!  which we call for its side-effect of incrementing the PC in the
!  correct manner for I64.
!

    !  This is a place to save the old FPE flags.
    !
    module saved_flags
    integer saved_old_flags
    end module

    !  Turn on underflow detection for one routine
    !  Using the FORTRAN library function FOR_SET_FPE.
    !
    subroutine via_fortran( rr )
    real rr
    include 'sys$library:fordef.for'
    include '($IEEEDEF)'
    integer new_flags, old_flags
    old_flags   = dfor$get_fpe()
    new_flags   = FPE_M_TRAP_UND
    call dfor$set_fpe( new_flags )
    !
    !  Code here uses new flags
    !
    rr = tiny(rr)
    type *,' Expect a catch #1'
    rr = rr / huge(rr)
    !
    call dfor$set_fpe( old_flags )
    return
    end subroutine

    !  Alternatively, do the same using the system routine.
    !
    subroutine via_system( rr )
    real rr
    include '($IEEEDEF)'
    integer*8 clear_flags, set_flags, old_flags, new_flags
    clear_flags = IEEE$M_MAP_DNZ + IEEE$M_MAP_UMZ
    set_flags   = IEEE$M_TRAP_ENABLE_UNF + IEEE$M_TRAP_ENABLE_DNOE
    call sys$ieee_set_fp_control(%ref(clear_flags), %ref(set_flags), %ref(old_flags))
    !
    !  Code here uses new flags
    !
    rr = tiny(rr)
    type *,' Expect a catch #2'
    rr = rr / huge(rr)
    !
    clear_flags = set_flags
    call sys$ieee_set_fp_control(%ref(clear_flags),%ref(old_flags),%ref(new_flags))
    return
    end subroutine

    !  Main program
    !
    program tester
    use saved_flags
    real, volatile :: r
    !
    !  Establish an exception handler.
    !
    external handler
    call lib$establish( handler )
    !
    !  Save the initial setting of the exception mode flags.
    !
    saved_old_flags = dfor$get_fpe()
    !
    !  This expression underflows, but because this program has
    !  been compiled with /IEEE=DENORM (by default) underflows
    !  do not raise exceptions.
    !
    write (6,100)
100 format(1x,' No catch expected')
    r = tiny(r);
    r = r / huge(r)

    !
    !  Call a routine to turn on underflow and try that expression
    !  again.  After the call, verify that underflow detection has
    !  been turned off.
    !
    call via_fortran( r )
    write (6,100)
    r = tiny(r)
    r = r / huge(r)
    !
    !  Ditto for the other approach
    !
    call via_system( r )
    write (6,100)
    r = tiny(r)
    r = r / huge(r)
    end program


    !  A handler is needed to catch any exceptions.
    !
    integer (kind=4) function handler( sigargs, mechargs )
    use saved_flags
    include '($CHFDEF)'
    include '($SSDEF)'
    integer sigargs(100)
    record /CHFDEF2/ mechargs
    integer  lib$match_cond
    integer  LIB$I64_INS_DECODE
    !
    integer index
    integer status
    integer no_loop      / 20 /
    logical int_over
    logical int_div
    logical float_over
    logical float_div
    logical float_under
    logical float_inval
    logical float_denorm
    logical HP_arith
    logical do_PC_advance
    integer   pc_index
    integer*8 pc_value
    !
    !  Don't loop forever between handler and exception
    !  (in case something goes wrong).
    !
    no_loop = no_loop - 1
    if( no_loop .le. 0 ) then
        handler = ss$_resignal
        return
    endif
    !
    !  We'll need the PC value of the instruction if
    !  this turns out to have been a fault rather than
    !  a trap.
    !
    pc_index = sigargs(1)
    pc_value = sigargs(pc_index)
    !
    !  Figure out what kind of exception we have, and
    !  whether it is a fault and we need to advance the
    !  PC before continuing.
    !
    do_PC_advance = .false.
    int_over      = .false.
    int_div       = .false.
    float_over    = .false.
    float_div     = .false.
    float_under   = .false.
    float_inval   = .false.
    float_denorm  = .false.
    HP_arith      = .false.
    !
    index = lib$match_cond(sigargs(2), SS$_INTOVF)
    if( index .eq. 0 ) then
        int_over = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_INTDIV)
    if( index .eq. 0 ) then
        int_div = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTOVF)
    if( index .eq. 0 ) then
        float_over = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTDIV)
    if( index .eq. 0 ) then
        float_div = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTUND)
    if( index .eq. 0 ) then
        float_under = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTOVF_F)
    if( index .eq. 0 ) then
        float_over = .true.
        do_PC_advance = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTDIV_F)
    if( index .eq. 0 ) then
        float_div = .true.
        do_PC_advance = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTUND_F)
    if( index .eq. 0 ) then
        float_under = .true.
        do_PC_advance = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTINV)
    if( index .eq. 0 ) then
        float_inval = .true.
        do_PC_advance = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_INTOVF_F)
    if( index .eq. 0 ) then
        int_over = .true.
        do_PC_advance = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_FLTDENORMAL)
    if( index .eq. 0 ) then
        float_denorm = .true.
    endif
    !
    index = lib$match_cond(sigargs(2), SS$_HPARITH)
    if( index .eq. 0 ) then
        HP_arith = .true.
    endif
    !
    !  Tell the user what kind of exception this is.
    !
    handler = ss$_continue
    if( float_over) then
        write(6,*) ' - Caught Floating overflow'

    else if ( int_over ) then
        write(6,*) ' - Caught Integer overflow'

    else if ( int_div ) then
        write(6,*) ' - Caught Integer divide by zero'

    else if ( float_div ) then
        write(6,*) ' - Caught Floating divide by zero'

    else if ( float_under ) then
        write(6,*) ' - Caught Floating underflow'

    else if ( float_inval ) then
        write(6,*) ' - Caught Floating invalid'

    else if ( HP_arith ) then
        write(6,*) ' - Caught HP arith error'

    else
        write(6,*) '  - Caught something else: resignal '
        !
        !  Here we have to restore the initial floating-point
        !  exception processing mode in case the exception
        !  happened during one of the times we'd changed it.
        !
        call dfor$set_fpe( saved_old_flags )
        !
        handler = ss$_resignal
    endif
    !
    !  If this was a fault, and we want to continue, then
    !  the PC has to be advanced over the faulting instruction.
    !
    if( do_PC_advance .and. (handler .eq. ss$_continue)) then
        status = lib$i64_ins_decode (pc_value)
        sigargs(pc_index) = pc_value
    endif
    !
    return
    end function handler