- THE PL/I NEWSLETTER -

Two of the great names in numerical methods - James Wilkinson of the National Physical Laboratory, U.K., and C. Reinsch of the Mathematical Institute of TH, Munchen, Germany - combined to assemble a masterful suite of numerical procedures in Handbook for Automatic Computation, Volume II, Linear Algebra, published by Springer-Verlag.

The algorithms cover linear equations, least squares, linear programming, and eigenvalue and eigenvector problems.
Note: in some of the following procedures, the original comment documentation and description of the parameters - which were missing - have been restored. As well, the original procedure names have been restored. The routines have not been re-tested.
Most of these routines have active links; the remainder will be added as soon as they have been converted.

• Cholesky symmetric decomposition of a positive definite matrix.
• Cholesky symmetric decomposition of a positive definite matrix, using compact storage. (Otherwise, it's the same as the previous item.)
• LDLt symmetric decomposition of a positive definite matrix.
• Iterative refinement of the solution of a positive definite system.
• Inversion of positive definite matrices by the Gauss-Jordan method. Two procedures are provided: GJDEF1 takes an n × n matrix as input. GJDEF2 takes the lower triangle stored in compact form as a vector of n(n+1)/2 elements.
• Symmetric decomposition of positive definite band matrices.
• The Conjugate gradient method for solving a system of linear equations that have symmetric positive definite coefficients.
• The solution of symmetric and unsymmetric band equations & the calculation of eigenvectors of band matrices.
• The solution of a real system of linear equations by the Crout method.
• The solution of a complex system of linear equations by the Crout method.
• Linear least squares solutions by Householder transformations. This algorithm can also be used for systems of linear equations of n equations in n unknowns, and for matrix inversion.
• Elimination with weighted row combinations.
  Four procedures are provided for solving least squares problems, and for solving n equations in n unknowns, and for matrix inversion. Choose the one best suited to your application:
  ORTHOLIN1 is for solving for X a least squares problem AX = B, where A is an n × m matrix, X is an m × k matrix, and B is an n × k matrix.
  ORTHOLIN2 is used when B is a single column.
  ORTHO1 is used for matrix inversion.
  ORTHO2 is for matrix inversion, and economises on storage.
• Singular value decomposition and least squares solutions.
• Simplex method based on triangular decompositions.
• The Jacobi method for a real symmetric matrix converts it to tri-diagonal form.
• Householder's tri-diaagonalization of a symmetric matrix. Several algorithms are provided, some requiring the matrix in compressed form, stored as a vector.

  Eigenvalues and eigenvectors

  • The Q R and Q L algorithms for symmetric matrices. TQL1 finds the Eigenvalues of a symmetric tri-diagonal matrix, while TQL2 finds the eigenvalues and eigenvectors of a symmetric tri-diagonal matrix
  • Using the implicit Q L algorithm, IMTQL1 will find all the eigenvalues of a symmetric tri-diagonal matrix, while IMTQL2 will find all the eigenvalues and eigenvectors of a symmetric tri-diagonal matrix.
  • Eigenvalues of a symmetric tri-diagonal matrix. Unsymmetric tri-diagonal matrices will also be handled, provided that they meet certain restrictions.
  • Using the Rational Q R transformation with Newton shift, RATQR finds the m lowest or the m highest eigenvalues of a symmetric tri-diagonal matrix.
  • The Q R algorithm for band symmetric matrices is designed to find one eigenvalue on each call to the procedure.
  • Tri-diagonalization of a symmetric band matrix. Any symmetric band matrix is transformed to tri-diagonal form by means of Jacobi rotations.
  • Simultaneous iteration method for symmetric matrices for obtaining eigenvalues and eigenvectors.
  • Reduction of the symmetric eigenproblem A x = lambda B x to standard form.
  • Balancing a matrix for calculation of eigenvalues and eigenvectors.
  • Solution to the eigenproblem by a norm reducing Jacobi type method.
  • Similarity reduction of a general matrix to Hessenberg form.
  • The Q R algorithm for real Hessenberg matrices produces all the eigenvalues.
  • Finds all the eigenvalues and eigenvectors of a complex matrix.
  • When only a few eigenvectors are required of a complex matrix, this will be found to be more efficient than the previous algorithm.
  • To find all the real and/or complex eigenvalues and eigenvectors of a real matrix A that has been reduced to upper Hessenberg form.
  • The modified L R algorithm for complex Hessenberg matrices finds all the eigenvalues.
  • Solution to the complex eigenproblem by a norm reducing Jacobi type method.
  • Several procedures are available in which specified eigenvectors are computed by inverse iterations.
    • Eigenvalues of a symmetric tridiagonal matrix when up to 25% of the eigenvalues and eigenectors are required. If more are desired, use this algorithm instead.
    • Eigenvalues of a real upper Hessenberg matrix.
    • Eigenvectors of a complex upper Hessenberg matrix from prescribed eigenvectors. When more than 25% of the eigenvectors are required, this procedure is more efficient.

PL/I offers a richer variety of basic data types than most languages. Consider just computational data. C has the floating-point types float, double, long double, integer types short [int] and long [int], and possibly char can be considered a computational data type.

PL/I computational data has base, scale, mode, and precision. Scale can be fixed or float, roughly corresponding to C. Mode can be real or complex; C has no complex data. Base can be binary or decimal; decimal data is unknown in C. Precision is the attribute that provides the most power, but also some degree of confusion. In place of C's short and long, for example, PL/I allows the programmer to specify the precision, the number of bits(binary), or digits(decimal) that are required to store the data, although the hardware can use more if appropriate. For example, a machine with a word size of 32 bits can store a fixed binary value up to FIXED BIN(31) [one bit is reserved for the sign and is not counted].

An often-unnoticed appendage to precision called the scale factor sometimes causes confusion. For fixed data precision is specified as a pair of numbers ( number-of-digits [, scale-factor ] ), where the scale-factor is assumed to be zero if not specified. Float data has no scale-factor, since by definition floating-point data includes its own scale. So far, so good. A number representing the number of items in inventory might be declared FIXED DECIMAL(8), or alternatively FIXED DECIMAL(8,0), which provides storage to handle quantities up to 99,999,999. The count must obviously be a whole number of items, and so the ",0". A number representing a bank balance might be declared FIXED DECIMAL(10,2), which would reserve storage to handle balances up to $99,999,999.99, dollars and cents. Non-US readers feel free to read this using your own currency conventions, PL/I supports them all.

The scale factor can cause some confusion when performing division as part of an arithmetic expression. C INT data is, as the name indicates, always integer. The compiler assures that expressions of type INT will always give integer results: in C 5/3=1. In PL/I, the scale factor is always taken into account, even when it is zero. For addition, subtraction, and multiplication, if all the operands have a scale factor of zero, the result will too. For division, however, PL/I always tries to preserve the maximum possible significance of a fractional quotient. Many books on PL/I fail to point out this fact. In practical terms, the compiler converts the dividend (the number being divided) to a temporary value which has just enough digits to hold its declared value, followed by zeroes to pad out the maximum possible length of that data type. The definition from the standard is:

m = N

n = N-p+q-s

Here "m" is the precision of the temporary result, "N" is the maximum number of bits or digits allowed for the data type, "p" is the number of bits or digits in the original value of the dividend, "q", and "s" are the scale factors of the dividend and the divisor.

This conversion is defined by the language, not the implementation. The only possible differences among machines might be in the maximum sizes of data which can be handled (N). Here are a couple of examples for IBM S/390 architecture:

FIXED BINARY(15,0) divided by FIXED BINARY(15,0): The temporary result is FIXED BINARY(31,15). That's right, the quotient will have fifteen bits of fractional result.

FIXED DECIMAL(5,0) divided by FIXED DECIMAL(5,2): The temporary result is FIXED DECIMAL(15,8).

Most PL/I "Language Reference" or similar manuals will have a table providing formulas for results of expressions involving different types of operands like "Table 16: Results of Arithmetic Operations for Coded Arithmetic Operands" in IBM's SC26-3114.

Having gotten to this point, a natural question is "so what?" In many cases this may not matter. For example 5/3 evaluates to 1.666.... Assigning this result to a FIXED variable with a scale factor of zero will truncate to an integer (1), exactly as in C. Assigning the result to a FLOAT variable, will result in the expected 1.666... to the limits of precision. Try this in C; the result will be 1.00000..., because the expression 5/3 involves integers and evaluates to an INT before assignment. To achieve the same result, C has to resort to casting, an art so arcane that I had to try three times to get it correct for my test. A C cast in this case is a method of forcing the conversion of data during expression evaluation. C has to use a cast to force the expression 5/3 to be evaluated as FLOAT. Guess which one of these is correct C:
f = (float)(5/3); f = (float)5/3; or f = (float)(5)/3;

The other situation in which you have to be aware of what PL/I is doing is when the division is part of a larger expression. Compare the statement i = (5/3) * 3 in PL/I and C. C divides five by three and gets one (integer division), then multiplies one by three and gets the result "three". PL/I divides five by three and gets 1.6666... (scaled integer division), it then multiplies 3.00000 (three scaled to the same number of fractional digits as the result of the division) and gets 4.999.... Assigning this result to i truncates to an integer "four". This may or may not be the desired result, but in any case you'd best be aware of what's happening. Another unobvious point is that PL/I will use decimal operations, since the constants "5" and "3" are decimal constants, the exact equivalent of the C expression would be (101B/11B) * 11B, or (BINARY(5)/BINARY(3))*BINARY(3).

PL/I provides a way for the programmer to exercise further control over division by way of the DIVIDE built-in function. The syntax is DIVIDE( dividend, divisor, precision [, scale] ). Divisor and dividend are self-explanatory. Here "precision" is the required precision of the result of the division, and "scale" is the required scale factor. As usual with PL/I built-in functions, the base, scale, and mode of the result are determined by the attributes of the dividend and the divisor. If the result is FLOAT the scale must be omitted. If the result is FIXED and the scale is omitted, zero is assumed.

The result of DIVIDE(5,3,15) [or DIVIDE(5,3,15,0)] will have a no fractional digits; an integer result of "one", or the same as C's (5/3). To obtain the same result as the C expression (5/3) * 3, code: DIVIDE(5,3,15) * 3

The spring 2000 edition of the COBOL and PL/I newsletter has some topics about PL/I. The specific PL/I topics may be viewed separately, or the entire newsletter in PDF format may be downloaded.
Of particular interest is the use of PL/I as a programming language for the Common Gateway Interface (CGI).

The PL/I Connection newsletter.

The PL/I Connection Newsletter No. 11 , December 1997.

The PL/I Connection Newsletter No. 10 , April 1997.

IBM enhanced VisualAge® PL/I with powerful features to increase development productivity, simplify the maintenance of your legacy code, and provide seamless portability from your host to your workstation. It also provides the tools needed to identify your Year 2000 date-related fields on both OS/2 and Windows NT.

VisualAge PL/I Enterprise Edition Version 2.1 combines the two separate offerings from the previous release of VisualAge PL/I (Standard and Professional) into a single offering available on both the OS/2 and Windows NT platforms. By doing so, it provides these productivity features:

• remote edit/compile
• redevelopment tools
• cross platform development
• year 2000 assessment strategies, find and fix, and test capabilities

Also included as an extra bonus offering is VisualAge CICS® Enterprise Application Development, which enables CICS host application development on the workstation.

• Product number 04L6564: VisualAge PL/1 Enterprise V2.1
• 04L6565: VisualAge PL/I Enterprise V2.1 Upgrade from Prior IBM PL/I Workstation or Competitive Products Program Package CD-ROM
• 04L6567: VisualAge PL/I Enterprise V2.1 Upgrade from Professional V2.0 Program Package CD-ROM
• 04L6566: VisualAge PL/I Enterprise V2.1 Upgrade from Standard V2.0 Package CD-ROM

Built-in functions are generic, and typically operate on any data type and for an array of any number of dimensions (up to the maximum permitted by the compiler).
If a user attempts to write a generic function, it is necessary to provide the procedures to perform the desired operation on each of the data types that he wants.
As well, if the argument(s) can be arrays, the user must also provide procedures to perform the desired operation on arrays of any number of dimensions. Each procedure dealing with an array must also be duplicated to deal with each of the different types that the user wants.
Thus, if there are four data types (integer, decimal, float, float (16)), and the compiler supports arrays having 10 dimensions, then 40 procedures must be provided.
It would be better to have a new attribute GENERIC_ALL which specifies that any procedure in the family is assumed to apply to scalars and arrays of any dimension. e.g.,