This section describes how to code your Fortran program using any of the ESSL run-time libraries.
In Fortran programs, most ESSL subroutines are invoked with the CALL statement:
CALL subroutine-name (argument-1, . . . , argument-n) |
An example of a calling sequence for the SAXPY subroutine might be:
CALL SAXPY (5,A,X,J+INC,Y,1) |
The remaining ESSL subroutines are invoked as functions by coding a function reference. You first declare the type of value returned by the function: short- or long-precision real, short- or long-precision complex, or integer. Then you code the function reference as part of an expression in a statement. An example of declaring and invoking the DASUM function might be:
DOUBLE PRECISION DASUM,SUM,X
.
.
.
SUM = DASUM (N,X,INCX)
|
Values are returned differently for ESSL subroutines and functions. For subroutines, the results of the computation are returned in an argument specified in the calling sequence. In the CALL statement above, the result is returned in argument Y. For functions, the result is returned as the value of the function. In the assignment statement above, the result is assigned to SUM.
See the Fortran publications for details on how to code the CALL statement and a function reference.
Some ESSL numerical quadrature subroutines call a user-supplied subroutine, subf, identified in the ESSL calling sequence. If your program that calls the numerical quadrature subroutines is coded in Fortran, there are some coding rules you must follow:
Table 27 lists the scalar data types in Fortran that are used for
ESSL. Only those types and lengths used by ESSL are listed.
Table 27. Scalar Data Types in Fortran Programs
| Terminology Used by ESSL | Fortran Equivalent |
|---|---|
| Character item1
'N', 'T', 'C' or 'n', 't', 'c' | CHARACTER*1
'N', 'T', 'C' |
| Logical item
.TRUE., .FALSE. | LOGICAL
.TRUE., .FALSE. |
| 32-bit environment integer
12345, -12345 | INTEGER or INTEGER*4
12345, -12345 |
| 64-bit environment integer2
12345, -12345 | INTEGER*83
12345_8, -12345_8 |
| Short-precision real number4
12.345 | REAL or REAL*4
0.12345E2 |
| Long-precision real number4
12.345 | DOUBLE PRECISION or REAL*8
0.12345D2 |
| Short-precision complex number4
(123.45, -54321.0) | COMPLEX or COMPLEX*8
(123.45E0, -543.21E2) |
| Long-precision complex number4
(123.45, -54321.0) | COMPLEX*16
(123.45D0, -543.21D2) |
2 In accordance with the LP64 data model, all ESSL integer arguments remain 32-bits except for the iusadr argument for ERRSET. 3 INTEGER may be used if you specify the compiler option -qintsize=8. 4 Short- and long-precision numbers look the same in this book. | |
Arrays are declared in Fortran by specifying the array name, the number of dimensions, and the range of each dimension in a DIMENSION statement or an explicit data type statement, such as REAL, DOUBLE PRECISION, and so forth.
Each array element can be either a real or complex data item of short or long precision. The type of the array determines the size of the element storage locations. Short-precision data requires 4 bytes, and long-precision data requires 8 bytes. Complex data requires two storage locations of either 4 or 8 bytes each, for short or long precision, respectively, to accommodate the two parts of the complex number: c = a+bi. Therefore, exactly twice as much storage is required for complex data as for real data of the same precision. See How Do You Set Up Your Scalar Data? for a description of real and complex numbers.
Even though complex data items require two storage locations, the same number of elements exist in the array as for real data. A reference to an element--for example, C(3)--in an array containing complex data gives you the whole complex number; that is, it contains both a and b, where the complex number is expressed as follows:
For a one-dimensional array in Fortran 77, you can code:
DIMENSION A(E1:E2)
where A is the name of the array, E1 is the lower bound, and E2 is the upper bound of the single dimension in the array. If the lower bound is not specified, such as in A(E2), the value is assumed to be 1. The upper bound is required.
A one-dimensional array is stored in ascending storage locations (relative to some base storage address) in the following order:
For example, the array A of length 4 specified in the DIMENSION statement as A(0:3) and containing the following elements:
A = (1, 2, 3, 4)
has its elements arranged in storage as follows:
For a two-dimensional array in Fortran 77, you can code:
DIMENSION A(E1:E2,F1:F2)
where A is the name of the array. E1 and F1 are the lower bounds of the first and second dimensions, respectively, and E2 and F2 are the upper bounds of the first and second dimensions, respectively. If either of the lower bounds is not specified, such as in A(E2,F1:F2), the value is assumed to be 1. The upper bounds are always required for each dimension. For examples of Fortran 77 usage, see SGEMV, DGEMV, CGEMV, ZGEMV, SGEMX, DGEMX, SGEMTX, and DGEMTX--Matrix-Vector Product for a General Matrix, Its Transpose, or Its Conjugate Transpose.
The elements of a two-dimensional array are stored in column-major order; that is, they are stored in the following ascending storage locations (relative to some base storage address) with the value of the first (row) subscript expression increasing most rapidly and the value of the second (column) subscript expression increasing least rapidly. Following are the locations of the elements in the array:
For example, the 3 by 4 array A specified in the DIMENSION statement as A(2:4,1:4) and containing the following elements:
* *
| 11 12 13 14 |
A = | 21 22 23 24 |
| 31 32 33 34 |
* *
has its elements arranged in storage as follows:
Each element A(I,J) of the array A, declared A(1:n, 1:m), containing real or complex data, occupies the storage location whose address is given by the following formula:
for:
where:
For a three-dimensional array in Fortran 77, you can code:
DIMENSION A(E1:E2,F1:F2,G1:G2)
where A is the name of the array. E1, F1, and G1 are the lower bounds of the first, second, and third dimensions, respectively, and E2, F2, and G2 are the upper bounds of the first, second, and third dimensions, respectively. If any of the lower bounds are not specified, such as in A(E1:E2,F1:F2,G2), the value is assumed to be 1. The upper bounds are always required for each dimension. For examples of Fortran 77 usage, see SCFT3 and DCFT3--Complex Fourier Transform in Three Dimensions.
The elements of a three-dimensional array can be thought of as a set of two-dimensional arrays, stored sequentially in ascending storage locations in the array. The elements in each two-dimensional array are stored as defined in the previous section. In the three-dimensional array, the value of the first (row) subscript expression increases most rapidly, the second (column) subscript expression increases less rapidly, and the third subscript expression (set of rows and columns) increases least rapidly. Following are the locations of the elements in the array:
* The last set is the G2-G1+1 set.
For example, the 3 by 2 by 4 array A specified in the DIMENSION statement as A(1:3,0:1,2:5) and containing the following sets of rows and columns of elements:
* * * * * * * *
| 111 121 | | 112 122 | | 113 123 | | 114 124 |
A = | 211 221 | | 212 222 | | 213 223 | | 214 224 |
| 311 321 | | 312 322 | | 313 323 | | 314 324 |
* * * * * * * *
has its elements arranged in storage as follows:
Each element A(I,J,K) of the array A, declared A(1:n, 1:m, 1:p), containing real or complex data, occupies the storage location whose address is given by the following formula:
for:
where:
The following example shows how to create up to a maximum of eight threads, where each thread calls the DURAND and DGEICD subroutines.
program matinv_example
implicit none
!
! program to invert m nxn random matrices
!
real(8), allocatable :: A(:,:,:), det(:,:), rcond(:)
real(8) :: dummy_aux, seed=1998, sd
integer :: rc, i, m=8, n=500, iopt=3, naux=0
!
! allocate storage
!
allocate(A(n,n,m),stat=rc)
call error_exit(rc,"Allocation of matrix A")
allocate(det(2,m),stat=rc)
call error_exit(rc,"Allocation of det")
allocate(rcond(m),stat=rc)
call error_exit(rc,"Allocation of rcond")
!
! Calculate inverses in parallel
!
!SMP$ parallel do private(i,sd), schedule(static),
!SMP$& share(n,a,iopt,rcond,det,dummy_aux,naux)
do i=1,m
sd = seed + 100*i
call durand(sd,n*n,A(1,1,i))
call dgeicd(A(1,1,i),n,n,iopt,rcond(i),det(1,i),
& dummy_aux,naux)
enddo
write(*,*)'Reciprocal condition numbers of the matrices are:'
write(*,'(4E12.4)') rcond
!
!
!
deallocate(A,stat=rc)
call error_exit(rc,"Deallocation of matrix A")
deallocate(det,stat=rc)
call error_exit(rc,"Deallocation of det")
deallocate(rcond,stat=rc)
call error_exit(rc,"Deallocation of rcond")
stop
contains
subroutine error_exit(error_code,string)
character(*) :: string
integer :: error_code
if(error_code .eq. 0 ) return
write(0,*)string,": failing return code was ",error_code
stop 1
end subroutine error_exit
end
|
ESSL provides you with flexibilities in handling both input-argument errors and computational errors:
Input-Argument Errors in Fortran and Computational Errors in Fortran explain how to use these facilities by describing the additional statements you must code in your program.
For multithreaded application programs, if you want to initialize the error option table and change the default settings for input-argument and computational errors, you need to implement the steps shown in Input-Argument Errors in Fortran and Computational Errors in Fortran on each thread that calls ESSL. An example is shown in Example of Handling Errors in a Multithreaded Application Program.
To obtain corrected input-argument values in a Fortran program and to avert program termination for the optionally-recoverable input-argument errors 2015, 2030, and 2200 add the statements in the following steps your program. Steps 3 and 7 for ERRSAV and ERRSTR, respectively, are optional. Adding these steps makes the effect of the call to ERRSET temporary.
EXTERNAL ENOTRM |
This declares the ESSL error exit routine ENOTRM as an external reference in your program. This should be coded in the beginning of your program before any of the following statements.
CALL EINFO (0) |
This calls the EINFO subroutine with one argument of value 0 to initialize the ESSL error option table. It is required only if you call ERRSET in your program. It is coded only once in the beginning of your program before any calls to ERRSET. For a description of EINFO, see EINFO--ESSL Error Information-Handler Subroutine.
CALL ERRSAV (ierno,tabent) |
(This is an optional step.) This calls the ERRSAV subroutine, which stores the error option table entry for error number ierno in an 8-byte storage area, tabent, which is accessible to your program. ERRSAV must be called for each entry you want to save. This step is used, along with step 7, for ERRSTR. For information on whether you should use ERRSAV and ERRSTR, see How Can You Control Error Handling in Large Applications by Saving and Restoring Entries in the Error Option Table?. For an example, see Example 3, as the use is the same as for computational errors.
CALL ERRSET (ierno,inoal,inomes,itrace,iusadr,irange) |
This calls the ERRSET subroutine, which allows you to dynamically modify the action taken when an error occurs. For optionally-recoverable ESSL input-argument errors, you need to call ERRSET only if you want to avoid terminating your program and you want the input arguments associated with this error to be assigned correct values in your program when the error occurs. For one error (ierno) or a range of errors (irange), you can specify:
ERRSET must be called for each error code you want to indicate as being recoverable. For ESSL, ierno should have a value of 2015, 2030 or 2200. If you want to eliminate error messages, you should indicate a negative number for inomes; otherwise, you should specify 0 for this argument. All the other ERRSET arguments should be specified as 0.
For a list of the default values set in the ESSL error option table, see Table 26. For a description of the input-argument errors, see Input-Argument Error Messages(2001-2099). For a description of ERRSET, see Chapter 17, Utilities.
CALL name (arg-1,...,arg-n,*yyy,*zzz,...) |
This calls the ESSL subroutine and specifies a branch on one or more return code values, where:
These are the statements at statement number yyy or zzz, shown in the CALL statement in Step 5, and preceded by an "*". The statement to which control is passed corresponds to the return code value for the error.
These statements perform whatever action is desired when the recoverable error occurs. These statements may check the new values set in the input arguments to determine whether adequate program storage is available, and then decide whether to continue or terminate the program. Otherwise, these statements may check that the size of the working storage arrays or the length of the transform agrees with other data in the program. The program may also store this corrected input argument value for future reference.
CALL ERRSTR (ierno,tabent) |
(This is an optional step.) This calls the ERRSTR subroutine, which stores an entry in the error option table for error number ierno from an 8-byte storage area, tabent, which is accessible to your program. ERRSTR must be called for each entry you want to store. This step is used, along with step 3, for ERRSAV. For information on whether you should use ERRSAV and ERRSTR, see How Can You Control Error Handling in Large Applications by Saving and Restoring Entries in the Error Option Table?. For an example, see Example 3, as the use is the same as for computational errors.
This example shows an error code 2015, which resets the size of the work area aux, specified in naux, if the value specified is too small. It also indicates that no error messages should be issued.
.
.
.
C DECLARE ENOTRM AS EXTERNAL
EXTERNAL ENOTRM
.
.
.
C INITIALIZE THE ESSL ERROR
C OPTION TABLE
CALL EINFO(0)
.
.
.
C MAKE ERROR CODE 2015 A RECOVERABLE
C ERROR AND SUPPRESS PRINTING ALL
C ERROR MESSAGES FOR IT
CALL ERRSET(2015,0,-1,0,ENOTRM,2015)
.
.
.
C CALL ESSL ROUTINE SWLEV.
C IF THE NAUX INPUT
C ARGUMENT IS TOO SMALL, ERROR
C 2015 OCCURS. THE MINIMUM VALUE
C REQUIRED IS STORED IN THE NAUX
C INPUT ARGUMENT AND CONTROL GOES
C TO LABEL 400.
CALL SWLEV(X,INCX,U,INCU,Y,INCY,N,AUX,NAUX,*400)
.
.
.
C CHECK THE RESULTING INPUT ARGUMENT
C VALUE IN NAUX AND TAKE THE
C DESIRED ACTION
400 .
.
.
|
To obtain information about an ESSL computational error in a Fortran program, add the statements in the following steps to your program. Steps 2 and 7 for ERRSAV and ERRSTR, respectively, are optional. Adding these steps makes the effect of the call to ERRSET temporary. For a list of those computational errors that return information and to which these steps apply, see EINFO--ESSL Error Information-Handler Subroutine.
CALL EINFO (0) |
This calls the EINFO subroutine with one argument of value 0 to initialize the ESSL error option table. It is required only if you call ERRSET in your program. It is coded only once in the beginning of your program before any calls to ERRSET. For a description of EINFO, see EINFO--ESSL Error Information-Handler Subroutine.
CALL ERRSAV (ierno,tabent) |
(This is an optional step.) This calls the ERRSAV subroutine, which stores the error option table entry for error number ierno in an 8-byte storage area, tabent, which is accessible to your program. ERRSAV must be called for each entry you want to save. This step is used, along with step 7, for ERRSTR. For information on whether you should use ERRSAV and ERRSTR, see How Can You Control Error Handling in Large Applications by Saving and Restoring Entries in the Error Option Table?.
CALL ERRSET (ierno,inoal,inomes,itrace,iusadr,irange) |
This calls the ERRSET subroutine, which allows you to dynamically modify the action taken when an error occurs. For ESSL computational errors, you need to call ERRSET only if you want to change the default values in the ESSL error option table. For one error (ierno) or a range of errors (irange), you can specify:
ERRSET must be called for each error code for which you want to change the default values. For ESSL, ierno should be set to one of the eligible values listed in Table 172. To allow your program to continue after an error in the specified range occurs, inoal must be set to a value greater than 1. For ESSL, iusadr should be specified as either 0 or 1 in a 32-bit environment (0_8 or 1_8 in a 64-bit environment), so a user exit is not taken.
For a list of the default values set in the ESSL error option table, see Table 26. For a description of the computational errors, see Computational Error Messages(2100-2199). For a description of ERRSET, see Chapter 17, Utilities.
CALL name (arg-1,...,arg-n,*yyy,*zzz,...) |
This calls the ESSL subroutine and specifies a branch on one or more return code values, where:
nmbr CALL EINFO (icode,inf1) -or- nmbr CALL EINFO (icode,inf1,inf2) |
This calls the EINFO subroutine, which returns information about certain computational errors, where:
These statements check the values returned in the output argument information receivers, inf1 and inf2, which contain the information about the computational error.
CALL ERRSTR (ierno,tabent) |
(This is an optional step.) This calls the ERRSTR subroutine, which stores an entry in the error option table for error number ierno from an 8-byte storage area, tabent, which is accessible to your program. ERRSTR must be called for each entry you want to store. This step is used, along with step 2, for ERRSAV. For information on whether you should use ERRSAV and ERRSTR, see How Can You Control Error Handling in Large Applications by Saving and Restoring Entries in the Error Option Table?.
This 32-bit environment example shows an error code 2104, which returns one piece of information: the index of the last diagonal with nonpositive value (I1).
.
.
.
C INITIALIZE THE ESSL ERROR
C OPTION TABLE
CALL EINFO(0)
.
.
.
C ALLOW 100 ERRORS FOR CODE 2104
CALL ERRSET(2104,100,0,0,0,2104)
.
.
.
C CALL ESSL ROUTINE DPPF.
C IF THE INPUT MATRIX IS NOT
C POSITIVE DEFINITE, CONTROL GOES TO
C LABEL 400
IOPT=0
CALL DPPF(APP,N,IOPT,*400)
.
.
.
C CALL THE INFORMATION-HANDLER
C ROUTINE FOR ERROR CODE 2104 TO
C RETURN ONE PIECE OF INFORMATION
C IN VARIABLE I1, THE INDEX OF THE
C LAST NONPOSITIVE DIAGONAL FOUND
C BY ROUTINE DPPF
400 CALL EINFO (2104,I1)
.
.
.
|
This 32-bit environment example shows an error code 2103, which returns one piece of information: the index of the zero diagonal (I1) found by DGEF.
.
.
.
C INITIALIZE THE ESSL ERROR
C OPTION TABLE
CALL EINFO(0)
.
.
.
C ALLOW 100 ERRORS FOR CODE 2103
CALL ERRSET(2103,100,0,0,0,2103)
.
.
.
C CALL ESSL SUBROUTINE DGEF.
C IF THE INPUT MATRIX IS
C SINGULAR, CONTROL GOES TO
C LABEL 400
CALL DGEF(A,LDA,N,IPVT,*400)
.
.
.
C CALL THE INFORMATION-HANDLER
C ROUTINE FOR ERROR CODE 2103 TO
C RETURN ONE PIECE OF INFORMATION
C IN VARIABLE I1, THE INDEX OF THE
C LAST ZERO DIAGONAL FOUND BY
C SUBROUTINE DGEF
400 CALL EINFO (2103,I1)
.
.
.
|
This 32-bit environment example shows an error code 2101, which returns two pieces of information: the eigenvalue (I1) that failed to converge after the indicated (I2) number of iterations. It uses ERRSAV and ERRSTR to insulate the effects of the error handling for error 2101 by this program.
.
.
C DECLARE AN AREA TO SAVE THE
C ERROR OPTION TABLE INFORMATION
C FOR ERROR CODE 2101
CHARACTER*8 SAV2101
.
.
C INITIALIZE THE ESSL ERROR
C OPTION TABLE
CALL EINFO(0)
C SAVE THE EXISTING ERROR OPTION
C TABLE ENTRY FOR ERROR CODE 2101
CALL ERRSAV(2101,SAV2101)
.
.
C ALLOW 255 ERRORS FOR CODE 2101
CALL ERRSET(2101,255,0,0,0,2101)
.
.
C CALL ESSL SUBROUTINE DGEEV.
C IF THE EIGENVALUE FAILED TO
C CONVERGE, CONTROL GOES TO LABEL 400
CALL DGEEV(IOPT,A,LDA,W,Z,LDZ,SELECT,N,AUX,NAUX,*400)
.
.
C CALL THE INFORMATION-HANDLER
C ROUTINE FOR ERROR CODE 2101 TO
C RETURN TWO PIECES OF INFORMATION.
C VARIABLE I1 CONTAINS THE EIGENVALUE
C THAT FAILED TO CONVERGE. VARIABLE
C I2 CONTAINS THE NUMBER OF ITERATIONS.
400 CALL EINFO (2101,I1,I2)
.
.
C RESTORE THE PREVIOUS ERROR OPTION
C TABLE ENTRY FOR ERROR CODE 2101.
C ERROR PROCESSING RETURNS TO HOW IT
C WAS BEFORE IT WAS ALTERED BY THE ABOVE
C ERRSET STATEMENT.
CALL ERRSTR(2101,SAV2101)
.
.
|
This 32-bit environment example shows how to modify the MATINV_EXAMPLE program in Creating Multiple Threads and Calling ESSL from Your Fortran Program with calls to the ESSL error handling subroutines. The ESSL error handling subroutines are called from each thread to: initialize the error option table, save the current error option table values for input-argument error 2015 and computational error 2105, change the default values for errors 2015 and 2105, and then restore the original default values for errors 2015 and 2105.
program matinv_example
implicit none
!
! program to invert m nxn random matrices
!
real(8), allocatable :: A(:,:,:), det(:,:), rcond(:)
real(8) :: dummy_aux, seed=1998, sd
integer :: rc, i, m=8, n=500, iopt=3, naux=0
integer :: inf1(8)
character(8) :: sav2015(8)
character(8) :: sav2105(8)
integer :: ENOTRM
!
external ENOTRM
!
! allocate storage
allocate(A(n,n,m),stat=rc)
call error_exit(rc,"Allocation of matrix A")
allocate(det(2,m),stat=rc)
call error_exit(rc,"Allocation of det")
allocate(rcond(m),stat=rc)
call error_exit(rc,"Allocation of rcond")
!
! Calculate inverses in parallel
!
!SMP$ parallel do private(i,sd), schedule(static),
!SMP$& share(n,m,a,iopt,rcond,det,dummy_aux,naux,sav2015,sav2105,inf1)
do i=1,m
!
! initialize error handling
call einfo(0)
!
! Save existing option table values for error 2015
call errsav(2015,sav2015(i))
!
! Set Error 2015 to be non-recoverable so dgeicd will dynamically
! allocate the work area.
call errset(2015,100,100,0,1,2015)
!
! Save existing option table values for error 2105
call errsav(2105,sav2105(i))
!
! Set Error 2105 to be recoverable
call errset(2105,100,100,0,ENOTRM,2105)
!
sd = seed + 100*i
call durand(sd,n*n,A(1,1,i))
call dgeicd(A(1,1,i),n,n,iopt,rcond(i),det(1,i),
& dummy_aux,naux,*10,*20)
10 goto 30
!
! Catch singular matrix returned by dgeicd.
20 CALL EINFO(2105,inf1(i))
WRITE(*,*) 'ERROR: Zero pivot found at location ',inf1(i)
!
! Restore the error option table entries
30 continue
call errstr(2015,SAV2015(i))
call errstr(2105,SAV2105(i))
enddo
|
write(*,*)'Reciprocal condition numbers of the matrices are:'
write(*,'(4E12.4)') rcond
!
deallocate(A,stat=rc)
call error_exit(rc,"Deallocation of matrix A")
deallocate(det,stat=rc)
call error_exit(rc,"Deallocation of det")
deallocate(rcond,stat=rc)
call error_exit(rc,"Deallocation of rcond")
stop
contains
subroutine error_exit(error_code,string)
character(*) :: string
integer :: error_code
if(error_code .eq. 0 ) return
write(0,*)string,": failing return code was ",error_code
stop 1
end subroutine error_exit
end
|