| Previous | Contents | Index |
This chapter describes:
This entire chapter pertains only to HP Fortran on OpenVMS Alpha systems only. |
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:
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:
The routines are described in Table 15-1.
| Name | Description |
|---|---|
| Basic Linear Algebra | The Basic Linear Algebra Subprograms (BLAS) library includes the industry-standard Basic Linear Algebra Subprograms for Level 1 (vector-vector, BLAS1), Level 2 (matrix-vector, BLAS2), and Level 3 (matrix-matrix, BLAS3). Also included are subprograms for BLAS Level 1 Extensions, and Sparse BLAS Level 1. |
| Signal Processing | The Signal Processing library provides a basic set of signal processing functions. Included are one-, two-, and three-dimensional Fast Fourier Transforms (FFT), group FFTs, Cosine/Sine Transforms (FCT/FST), Convolution, Correlation, and Digital Filters. |
| Sparse Linear System | The Sparse Linear System library provides both direct and iterative sparse linear system solvers. The direct solver package supports both symmetric and symmetrically shaped sparse matrices stored using the compressed row storage scheme. The iterative solver package supports a basic set of storage schemes, preconditioners, and iterative solvers. |
| LAPACK | LAPACK is an industry-standard subprogram package offering an extensive set of linear system and eigenproblem solvers. LAPACK uses blocked algorithms that are better suited to most modern architectures, particularly ones with memory hierarchies. |
| Utility subprograms | Utility subprograms include random number generation, vector math functions, and sorting subprograms. |
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
|
The differences are:
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:
There are a number of implications for this approach:
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:
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:
HP does not encourage users to change the exception-handling or rounding mode within a program. This practice is particularly discouraged for an application using VAX-format floating-point. |
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:
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.
| Example A-1 Changing Floating-Point Exception Mode |
|---|
! 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
|
| Previous | Next | Contents | Index |