Organization of the Documentation

This manual contains a concise description of each routine, with at least one demonstrated example of each routine, including sample input and results. You will find all information pertaining to the MATH/LIBRARY in this manual. Moreover, all information pertaining to a particular routine is in one place within a chapter.

Each chapter begins with an introduction followed by a table of contents that lists the routines included in the chapter. Documentation of the routines consists of the following information:

Required Arguments: a description of the required arguments in the order of their occurrence. Input arguments usually occur first, followed by input/output arguments, with output arguments described last. Futhermore, the following terms apply to arguments:

Input[/Output] Argument must be initialized; the routine may return output through this argument based on other optional data the user may choose to pass to this routine; cannot be a constant or an expression.

The names of the routines are mnemonic and unique. Most routines are available in both a single precision and a double precision version, with names of the two versions sharing a common root. The root name is also the generic interface name. The name of the double precision specific version begins with a “D_” and the single precision specific version begins with an “S_”. For example, the following pairs are precision specific names of routines in the two different precisions: S_GQRUL/D_GQRUL (the root is “GQRUL ,” for “Gauss quadrature rule”) and S_RECCF/D_RECCF (the root is “RECCF,” for “recurrence coefficient”). The precision specific names of the IMSL routines that return or accept the type complex data begin with the letter “C_” or “Z_” for complex or double complex, respectively. Of course, the generic name can be used as an entry point for all precisions supported.

When this convention is not followed the generic and specific interfaces are noted in the documentation. For example, in the case of the BLAS and trigonometric intrinsic functions where standard names are already established, the standard names are used as the precision specific names. There may also be other interfaces supplied to the routine to provide for backwards compatibility with previous versions of the IMSL Fortran Numerical Library. These alternate interfaces are noted in the documentation when they are available.

Except when expressly stated otherwise, the names of the variables in the argument lists follow the Fortran default type for integer and floating point. In other words, a variable whose name begins with one of the letters “I” through “N” is of type INTEGER, and otherwise is of type REAL or DOUBLE PRECISION , depending on the precision of the routine.

An assumed-size array with more than one dimension that is used as a Fortran argument can have an assumed-size declarator for the last dimension only. In the MATH/LIBRARY routines, the information about the first dimension is passed by a variable with the prefix “LD” and with the array name as the root. For example, the argument LDA contains the leading dimension of array A. In most cases, information about the dimensions of arrays is obtained from the array through the use of Fortran 90’s size function. Therefore, arguments carrying this type of information are usually defined as optional arguments.

Where appropriate, the same variable name is used consistently throughout a chapter in the MATH/LIBRARY. For example, in the routines for random number generation, NR denotes the number of random numbers to be generated, and R or IR denotes the array that stores the numbers.

When writing programs accessing the MATH/LIBRARY, the user should choose Fortran names that do not conflict with names of IMSL subroutines, functions, or named common blocks. The careful user can avoid any conflicts with IMSL names if, in choosing names, the following rules are observed:

For further details, see the section on Reserved Names in the Reference Material.

The documentation for the routines uses the generic name and omits the prefix, and hence the entire suite of routines for that subject is documented under the generic name.

Examples that appear in the documentation also use the generic name. To further illustrate this principle, note the LIN_SOL_GEN documentation (see Chapter 1, “Linear Systems”), for solving general systems of linear algebraic equations. A description is provided for just one data type. There are four documented routines in this subject area: s_lin_sol_gen, d_lin_sol_gen, c_lin_sol_gen, and z_lin_sol_gen.

These routines constitute single-precision, double-precision, complex, and double-complex precision versions of the code.

The Fortran 90 compiler identifies the appropriate routine. Use of a module is required with the routines. The naming convention for modules joins the suffix “_int” to the generic routine name. Thus, the line “use lin_sol_gen_int” is inserted near the top of any routine that calls the subprogram “lin_sol_gen”. More inclusive modules are also available, such as imsl_libraries and numerical libraries. To avoid name conflicts, Fortran 90 permits re-labeling names defined in modules so they do not conflict with names of routines or variables in the user’s program. The user can also restrict access to names defined in IMSL Library modules by use of the “: ONLY, <list of names>” qualifier.

where the overstrike denotes complex conjugation. IMSL Fortran Numerical Library linear algebra software uses this convention to conserve the utility of generic documentation for that code subject. All references to orthogonal matrices are to be replaced by their complex counterparts, unitary matrices. Thus, an n × n orthogonal matrix Q satisfies the condition

. An n × n unitary matrix V satisfies the analogous condition for complex matrices,

In general, the IMSL MATH/LIBRARY codes are written so that computations are not affected by underflow, provided the system (hardware or software) places a zero value in the register. In this case, system error messages indicating underflow should be ignored.

IMSL codes are also written to avoid overflow. A program that produces system error messages indicating overflow should be examined for programming errors such as incorrect input data, mismatch of argument types, or improper dimensioning.

In many cases, the documentation for a routine points out common pitfalls that can lead to failure of the algorithm.

Library routines detect error conditions, classify them as to severity, and treat them accordingly. This error-handling capability provides automatic protection for the user without requiring the user to make any specific provisions for the treatment of error conditions. See the section on User Errors in the Reference Material for further details.

Users are required to incorporate a “use” statement near the top of their program for the IMSL routine being called when writing new code that uses this library. However, legacy code which calls routines in the previous version of the library without the use of a “use” statement will continue to work as before. Also, code that employed the “use numerical_libraries” statement from the previous version of the library will continue to work properly with this version of the library.

Users wishing to update existing programs so as to call other routines from this library should incorporate a use statement for the specific new routine being called. (Here, the term “new routine” implies any routine in the library, only “new” to the user’s program.) Use of the more encompassing “imsl_libraries” module in this case could result in argument mismatches for the “old” routine(s) being called. (The compiler would catch this.)

Users wishing to update existing programs to call the new generic versions of the routines must change their calls to the existing routines to match the new calling sequences and use either the routine specific interface modules or the all-encompassing “imsl_libraries” module.

Users of the IMSL Fortran Numerical Library benefit by having a standard (MPI) Message Passing Interface environment. This is needed to accomplish parallel computing within parts of the Library. Either of the icons above clues the reader when this is the case. If parallel computing is not required, then the IMSL Library suite of dummy MPI routines can be substituted for standard MPI routines. All requested MPI routines called by the IMSL Library are in this dummy suite. Warning messages will appear if a code or example requires more than one process to execute. Typically users need not be aware of the parallel codes.

NOTE: that a standard MPI environment is not part of the IMSL Fortran Numerical Library. The standard includes a library of MPI Fortran and C routines, MPI “include” files, usage documentation, and other run-time utilities.

NOTE: Details on linking to the appropriate libraries are explained in the online README file of the product distribution.

There are three situations of MPI usage in the IMSL Fortran Numerical Library:

1. There are some computations that are performed with the ‘box’ data type that benefit from the use of parallel processing. For computations involving a single array or a single problem, there is no IMSL use of parallel processing or MPI codes. The box type data type implies that several problems of the same size and type are to be computed and solved. Each rack of the box is an independent problem. This means that each problem could potentially be solved in parallel. The default for computing a box data type calculation is that a single processor will do all of the problems, one after the other. If this is acceptable there should be no further concern about which version of the libraries is used for linking. If the problems are to be solved in parallel, then the user must link with a working version of an MPI Library and the appropriate IMSL Library. Examples demonstrating the use of box type data may be found in Chapter 10, “Linear Algebra Operators and Generic Functions”.

2. Various routines in Chapter 1, “Linear Systems” allow the user to interface with the ScaLAPACK Library routines. If the user chooses to run on only one processor then these routines will utilize either IMSL Library code or LAPACK Library code based on the libraries the user chooses to use during linking. If the user chooses to run on multiple processors then working versions of MPI, ScaLAPACK, PBLAS, and Blacs will need to be present. These routines are marked with the MPI Capable icon.

3. There are some routines or operators in the Library that require that a working MPI Library be present in order for them to run. Examples are the large-scale parallel solvers and the ScaLAPACK utilities. Routines of this type are marked with the MPI Required icon. For these routines, the user must link with a working version of an MPI Library and the appropriate IMSL Library.

In all cases described above it is the user’s responsibility to supply working versions of the aforementioned third party libraries when those libraries are required.

Table 1 below lists the chapters and IMSL routines calling MPI routines or the replacement non-parallel package.

Table 1 — IMSL Routines Calling MPI Routines or Replacement Non-Parallel Package
Chapter Name and Number	Routine with MPI Utilized
Linear Systems, 1	PARALLEL_NONNEGATIVE_LSQ
Linear Systems, 1	PARALLEL_BOUNDED_LSQ
Linear Systems, 1	Those routines which utilize ScaLAPACK listed in Table D below.
Linear Algebra and Generic Functions, 10	See entire following Table 2, “Defined Operators and Generic Functions for Dense Arrays.”
Utilities, 11	ScaLAPACK_SETUP
Utilities, 11	ScaLAPACK_GETDIM
Utilities, 11	ScaLAPACK_READ
Utilities, 11	ScaLAPACK_WRITE
Utilities, 11	ScaLAPACK_MAP
Utilities, 11	ScaLAPACK_UNMAP
Utilities, 11	ScaLAPACK_EXIT
Reference Material	Entire Error Processor Package for IMSL Library, if MPI is utilized

Each subject routine called or otherwise referenced requires the “use” statement for an interface block designed for that subject routine. The contents of this interface block are the interfaces to the separate routines available for that subject. Packaged descriptive names for option numbers that modify documented optional data or internal parameters might also be provided in the interface block. Although this seems like an additional complication, many errors are avoided at an early stage in development through the use of these interface blocks. The “use” statement is required for each routine called in the user’s program. As illustrated in Examples 3 and 4 in routine lin_geig_gen, the “use” statement is required for defining the secondary option flags.

The function subprogram for s_NaN() or d_NaN() does not require an interface block because it has only a single “required” dummy argument. Also, if one is only using the Fortran 77 interfaces supplied for backwards compatibility then the “use” statements are not required.

IMSL Fortran Numerical Library routines have required arguments and may have optional arguments. All arguments are documented for each routine. For example, consider the routine lin_sol_gen that solves the linear algebraic matrix equation Ax = b. The required arguments are three rank-2 Fortran 90 arrays: A, b, and x. The input data for the problem are the A and b arrays; the solution output is the x array. Often there are other arguments for this linear solver that are closely connected with the computation but are not as compelling as the primary problem. The inverse matrix A-1 may be needed as part of a larger application. To output this parameter, use the optional argument given by the “ainv=” keyword. The rank-2 output array argument used on the right-hand side of the equal sign contains the inverse matrix. See Example 2 of LIN_SOL_GEN in Chapter 1, “Linear Systems” for an example of computing the inverse matrix.

For compatibility with previous versions of the IMSL Libraries, the NUMERICAL_LIBRARIES interface module includes backwards-compatible positional argument interfaces to all routines that existed in the Fortran 77 version of the Library. Note that it is not necessary to include “use” statements when calling these routines by themselves. Existing programs that called these routines will continue to work in the same manner as before.

Some of the primary routines have arguments “epack=” and “iopt=”. As noted the “epack=” argument is of derived type s_error or d_error. The prefix “s_” or “d_” is chosen depending on the precision of the data type for that routine. These optional arguments are part of the interface to certain routines, and are used to modify internal algorithm choices or other parameters.

This additional optional argument (available for some routines) is further distinguished—a derived type array that contains a number of parameters to modify the internal algorithm of a routine. This derived type has the name ?_options, where “?_” is either “s_” or “d_”. The choice depends on the precision of the data type. The declaration of this derived type is packaged within the modules for these codes.

The definition of the derived types is:

where the “?_” is either “s_” or “d_”, and the kind value matches the desired data type indicated by the choice of “s” or “d”.

Example 3 of LIN_SOL_GEN in Chapter 1, “Linear Systems” illustrates the use of iterative refinement to compute a double-precision solution based on a single-precision factorization of the matrix. This is communicated to the routine using an optional argument with optional data. For efficiency of iterative refinement, perform the factorization step once, and then save the factored matrix in the array A and the pivoting information in the rank-1 integer array, ipivots. By default, the factorization is normally discarded. To enable the routine to be re-entered with a previously computed factorization of the matrix, optional data are used as array entries in the “iopt=” optional argument. The packaging of LIN_SOL_GEN includes the definitions of the self-documenting integer parameters lin_sol_gen_save_LU and lin_sol_gen_solve_A. These parameters have the values 2 and 3, but the programmer usually does not need to be aware of it.
The following rules apply to the “iopt=iopt” optional argument:

2. The data structure for the optional data array has the following form:
iopt (IO) = ?_options (Option_number, Optional_data)
[iopt (IO + 1) =?_options (Option_number, Optional_data)]
The length of the data set is specified by the documentation for an individual routine. (The Optional_data is output in some cases and may not be used in other cases.) The square braces […] denote optional items.

Illustration: Example 3 of LIN_EIG_SELF in Chapter 2, “Singular Value and Eigenvalue Decomposition”, a new definition for a small diagonal term is passed to lin_sol_self. There is one line of code required for the change and the new tolerance:

3. The internal processing of option numbers stops when Option_number == 0 or when IO > SIZE(iopt). This signals each routine having this optional argument that all desired changes to default values of internal parameters have been made. This implies that the last option number is the value zero or the value of SIZE(iopt) matches the last optional value changed.

Option numbers can be written in any order, and any selected set of options can be changed from the defaults. They may be repeated. Example 3 in of LIN_SOL_SELF in Chapter 1, “Linear Systems” uses three and then four option numbers for purposes of computing an eigenvector associated with a known eigenvalue.

To assist users in writing compact and readable code, the IMSL Fortran Numerical Library provides overloaded assignment and logical operations for the derived types s_options, d_options, s_error, and d_error. Each of these derived types has an individual record consisting of an integer and a floating-point number. The components of the derived types, in all cases, are named idummy followed by rdummy. In many cases, the item referenced is the component idummy. This integer value can be used exactly as any integer by use of the component selector character (%). Thus, a program could assign a value and test after calling a routine:

Using the overloaded assignment and logical operations, this code fragment can be written in the equivalent and more readable form:

Generally the assignments and logical operations refer only to component idummy. The assignment “s_epack(1)=0” is equivalent to “s_epack(1)=s_error(0,0E0)”. Thus, the floating-point component rdummy is assigned the value 0E0. The assignment statement “I=s_epack(1)”, for I an integer type, is equivalent to “I=s_epack(1)%idummy”. The value of component rdummy is ignored in this assignment. For the logical operators, a single element of any of the IMSL Fortran Numerical Library derived types can be in either the first or second operand.

Derived Type	Overloaded Assignments and Tests
s_options	I=s_options(1);s_options(1)=I	= =	/=	<	<=	>	>=
d_options	I=d_options(1);d_options(1)=I	= =	/=	<	<=	>	>=
s_epack	I=s_epack(1);s_epack(1)=I	= =	/=	<	<=	>	>=
d_epack	I=d_epack(1);d_epack(1)=I	= =	/=	<	<=	>	>=

In the examples, operator_ex01, … , _ex37, the overloaded assignments and tests have been used whenever they improve the readability of the code.

The routines in the IMSL MATH/LIBRARY attempt to detect and report errors and invalid input. Errors are classified and are assigned a code number. By default, errors of moderate or worse severity result in messages being automatically printed by the routine. Moreover, errors of worse severity cause program execution to stop. The severity level and the general nature of the error are designated by an “error type” ranging from 0 to 5. An error type 0 is no error; types 1 through 5 are progressively more severe. In most cases, you need not be concerned with our method of handling errors. For those interested, a complete description of the error-handling system is given in the Reference Material, which also describes how you can change the default actions and access the error code numbers.

A separate error handler is provided to allow users to handle errors of differing types being reported from several nodes without danger of “jumbling” or mixing error messages. The design of this error handler is described more fully in Hanson (1992). The primary feature of the design is the use of a separate array for each parallel call to a routine. This allows the user to summarize errors using the routine error_post in a non-parallel part of an application. For a more detailed discussion of the use of this error handler in applications which use MPI for distributed computing, see the Reference Material.

Most of the routines in the IMSL MATH/LIBRARY (except the line printer routines and special utility routines) do not print any of the results. The output is returned in Fortran variables, and you can print these yourself. See Chapter 11, “Utilities” for detailed descriptions of these routines.

A commonly used routine in the examples is the IMSL routine UMACH (see the Reference Material), which retrieves the Fortran device unit number for printing the results. Because this routine obtains device unit numbers, it can be used to redirect the input or output. The section on Machine-Dependent Constants in the Reference Material contains a description of the routine UMACH.

The IMSL Fortran Numerical Library contains routines which take advantage of Fortran 90 language constructs, including Fortran 90 array data types. One feature of the design is that the default use may be as simple as the problem statement. Complicated, professional-quality mathematical software is hidden from the casual or beginning user.

In addition, high-level operators and functions are provided in the Library. They are described in Chapter 10, “Linear Algebra Operators and Generic Functions”.

The IMSL Fortran Numerical Library allows users to leverage the high-performance technology of shared memory parallelism (SMP) when their environment supports it. Support for SMP systems within the IMSL Library is delivered through various means, depending upon the availability of technologies such as OpenMP, high performance LAPACK and BLAS, and hardware-specific IMSL algorithms. Use of the IMSL Fortran Numerical Library on SMP systems can be achieved by using the appropriate link environment variable when building your application. Details on the available link environment variables for your installation of the IMSL Fortran Numerical Library can be found in the online README file of the product distribution.

The IMSL Fortran Numerical Library is thread-safe in those environments that support OpenMP. This was achieved by using OpenMP directives that define global variables located in the code so they are private to the individual threads. Thread safety allows users to create instances of routines running on multiple threads and to include any routine in the IMSL Fortran Numerical Library in these threads.

For users who are primarily interested in easy-to-use software for numerical linear algebra, see Chapter 10, “Linear Algebra Operators and Generic Functions”. This compact notation for writing Fortran 90 programs, when it applies, results in code that is easier to read and maintain than traditional subprogram usage.

Users may begin their code development using operators and generic functions. If a more efficient executable code is required, a user may need to switch to equivalent subroutine calls using IMSL Fortran Numerical Library routines.

Table 2 and Table 3 contain lists of the defined operators and some of their generic functions.

Table 2 — Defined Operators and Generic Functions for Dense Arrays
Defined Array Operation	Matrix Operation
A .x. B	AB
.i. A	A-1
.t. A, .h. A	AT,A*
A .ix. B	A-1B
B .xi. A	BA-1
A .tx. B, or (.t. A) .x. B A .hx. B, or (.h. A) .x. B	ATB,A*B
B .xt. A, or B .x. (.t. A) B .xh. A, or B .x. (.h. A)	BAT,BA*
S=SVD(A [,U=U, V=V])	A = USVT
E=EIG(A [[,B=B, D=D], V=V, W=W])	(AV = VE), AVD = BVE, (AW = WE), AWD = BWE
R=CHOL(A)	A = RTR
Q=ORTH(A [,R=R])	(A = QR),QTQ = I
U=UNIT(A)	[u1,…] = [a1/∥a1∥,…]
F=DET(A)	det(A) = determinant
K=RANK(A)	rank(A) = rank
P=NORM(A[,[type=]i])
C=COND(A)
Z=EYE(N)	Z = IN
A=DIAG(X)	A = diag(x1,…)
X=DIAGONALS(A)	x = (x11,…)
W=FFT(Z); Z=IFFT(W)	Discrete Fourier Transform, Inverse
A=RAND(A)	random numbers, 0 < A < 1
L=isNaN(A)	test for NaN, if (l) then…

Table 3 — Defined Operators and Generic Functions for Harwell-Boeing Sparse Matrices
Defined Operation	Matrix Operation
Data Management	Define entries of sparse matrices
A .x. B	AB
.t. A, .h. A	AT,A*
A .ix. B	A-1B
B .xi. A	BA-1
A .tx. B, or (.t. A) .x. B A .hx. B, or (.h. A) .x. B	ATB,A*B
B .xt. A, or B .x. (.t. A) B .xh. A, or B .x. (.h. A)	BAT,BA*
A+B	Sum of two sparse matrices
C=COND(A)