Most of the techniques in LAPACK are based on a matrix factorization as the first step. There are two main types of factorization forms:

Generally, users do not have to know the details of how the factorization is stored, because other LAPACK routines manipulate the factored form.

Regardless of the form of the factorization, it reduces the solution phase to one that requires only permutations and the solution of triangular systems. For example, the LU factorization of a general matrix, , is used to solve for X in the system of equations by successively applying the inverses of P, L, and U to the right-hand side:

In the last two steps, the inverse of the triangular factors is not computed, but triangular systems of the form and are solved instead.

The following table lists the factorization forms for each of the factorization routines for real matrices. The factorization forms differ for DGETRF and DGBTRF, even though both compute an LU factorization with partial pivoting. You can also obtain the same factorizations through the LAPACK driver routines (for instance, DGESV or DGESVX).

Example 3-1. LU factorization

The DGETRF subroutine performs an LU factorization with partial pivoting () as the first step in solving a general system of linear equations . If DGETRF is called with the following:

Equation 3-18.

details of the factorization are returned, as follows:

Equation 3-19.

Matrices L and U are given explicitly in the lower and upper triangles, respectively, of A:

Equation 3-20.

The IPIV vector specifies the row interchanges that were performed. IPIV(1)= 3 implies that the first and third rows were interchanged when factoring the first column; IPIV(2)= 3 implies that the second and third rows were interchanged when factoring the second column. In this case, IPIV(3) must be 3 because there are only three rows. Thus, the permutation matrix is the following:

Equation 3-21.

Generally, the pivot information is used directly from IPIV without constructing matrix P.

Example 3-2. Symmetric indefinite matrix factorization

DSYTRF factors a symmetric indefinite matrix A into one of the forms or , where L and U are lower triangular and upper triangular matrices, respectively, in factored form, and D is a diagonal matrix with 1-by-1 and 2-by-2 diagonal blocks. To illustrate this factorization, choose a symmetric matrix that requires both 1-by-1 and 2-by-2 pivots:

Equation 3-22.

Only the lower triangle of A is specified because the matrix is symmetric, but you could have specified the upper triangle instead. The output from DSYTRF is the following:

Equation 3-23.

The signs of the indices in the IPIV vector indicate that a 2-by-2 pivot block was used for the first two columns, and 1-by-1 pivots were used for the third and fourth columns. Therefore, D must be the following:

Equation 3-24.

Matrix L is supplied in factored form as , where the parts of each that differ from the identity are stored in A below their corresponding blocks :

Equation 3-25.

Equation 3-26.

The LAPACK routines always check the arguments on entry for incorrect values. If an illegal argument value is detected, the error-handling subroutine XERBLA is called. XERBLA prints a message similar to the following to standard error, and then it aborts:

** On entry to DGETRF parameter number 4 had an illegal value

All other errors in the LAPACK routines are described by error codes returned in info, the last argument. The values returned in info are routine-specific, except for info = 0, which always means that the requested operation completed successfully.

For example, an error code of info > 0 from DGETRF means that one of the diagonal elements of the factor U from the factorization is exactly 0. This indicates that one of the rows of A is a linear combination of the other rows, and the linear system does not have a unique solution.

Example 3-3. Error conditions

If DGETRF is given the matrix

Equation 3-27.

it returns

Equation 3-28.

which corresponds to the factorization

Equation 3-29.

On exit from DGETRF, info = 3, indicating that U(3,3) is exactly 0. This is not an error condition for the factorization because the factors that were computed satisfy , but the factorization cannot be used to solve the system.

A return code of info = 0 from a factorization routine indicates that the triangular factors have nonzero diagonals. The linear system still may be too ill-conditioned to give a meaningful solution.

One indicator that you can examine before computing the solution is the reciprocal condition number, RCOND. The condition number, defined as , tells how much the relative errors in A and b are magnified in the solution x. DGECON and the other condition estimation routines compute by using the exact 1-norm or infinity-norm of A and an estimate for the norm of A ^-1, because computing the inverse explicitly would be very expensive, and the inverse may not even exist. By convention, if A^-1 does not exist, and RCOND should be computed as 0 or a small number on the order of , the machine epsilon.

If the condition number is large (that is, if RCOND is small), small errors in A and b may lead to large errors in the solution x. The rule of thumb is that the solution loses one digit of accuracy for every power of 10 in the condition number, assuming that the elements of A all have about the same magnitude. For example, a condition number of 100 (RCOND = 0.01) implies that the last 2 digits are inaccurate; a condition number of (RCOND < , the machine epsilon which is approximately on Altix systems) implies that all of the digits have been lost. This value for the machine epsilon assumes a model of rounded floating-point arithmetic with base and mantissa digits, and .

The expert driver routine DGESVX uses this rule of thumb to decide whether the solution and error bounds should be computed. If RCOND is less than the machine epsilon, DGESVX returns info = N+1, indicating that the matrix A is singular to working precision, and it does not compute the solution.

Example 3-4. Roundoff errors

The matrix

Equation 3-30.

is singular in exact arithmetic, but on Altix systems, DGETRF returns

Equation 3-31.

where IPIV=[3,3,3] and info = 0. In exact arithmetic, A(3,3) would have been 0, but roundoff error has made this entry instead. The reciprocal condition number computed by DGECON is , which is less than the machine epsilon of . Therefore, DGESVX returns info= 4 and does not try to solve any systems with this A.

You can use the condition number to compute a simple bound on the relative error in the computed solution to a system of equations (see Introduction to Matrix Computations , by Stewart). If x is the exact solution and is the computed solution, let r be the residual . If A is nonsingular:

Equation 3-32.

and

Equation 3-33.

From Equation 3-33 we can already see that if is large, then the error in may be significantly greater than the residual r.

Since , it follows that , therefore

Equation 3-34.

Define the condition number . This gives the bound

Equation 3-35.

For a description of a more precise error bound based on a component-wise error analysis, see “Error Bounds”.

Another application of the condition number is to consider the computed solution to be the exact solution of a slightly perturbed linear system:

Equation 3-36.

where ΔA is small in norm with respect to A, and Δb is small in norm with respect to b. For Gaussian elimination with partial pivoting, it has been shown that and , where ω is the product of ε and a slowly growing function of n (see The Algebraic Eigenvalue Problem, by Wilkinson). Proving that an algorithm has this property is the stuff of backward error analysis, and it ensures that, if the problem is well-conditioned, the computed solution is near the exact solution of the original problem. Because (Equation 3-34) simplifies to

Equation 3-37.

Assuming A is nonsingular,

Equation 3-38.

Taking norms and dividing by ∥x∥,

Equation 3-39.

Using the inequality ,

Equation 3-40.

and substituting ,

Equation 3-41.

provided . In terms of the relative backward error ,

Equation 3-42.

In “Error Bounds”, the backward error is defined slightly differently to obtain a component-wise error bound.

The condition number defined in the last section is sensitive to the scaling of A. For example, the matrix

Equation 3-43.

has as its inverse

Equation 3-44.

and so RCOND= . If this value of RCOND is less than the machine epsilon on a system, DGESVX (with FACT = `N') does not try to solve a system with this matrix. However, A has elements that vary widely in magnitude, so the bounds on the relative error in the solution may be pessimistic. For example, if the right-hand side in is the following:

Equation 3-45.

DGETRF followed by DGETRS produces the exact answer, .

You can improve the condition of this example by a simple row scaling. Scaling a problem before computing its solution is known as equilibration, and it is an option to some of the expert driver routines (those for general or positive definite matrices). Enabling equilibration does not necessarily mean it will be done; the driver routine will choose to do row scaling, column scaling, both row and column scaling, or no scaling, depending on the input data. The usage of this option is as follows:

      CALL DGESVX('E', 'N', N, NRHS, A, LDA, AF,
$         LDAF, IPIV, EQUED, R, C, B, LDB, X, LDX,
$         RCOND, FERR, BERR, WORK, IWORK, INFO)

The 'E' in the first argument enables equilibration. For this example, EQUED = 'R' on return, indicating that only row scaling was done, and the vector R contains the scaling constants:

Equation 3-46.

The form of the equilibrated problem is:

Equation 3-47.

or , where is returned in A:

Equation 3-48.

and is returned in B:

Equation 3-49.

The factored form, AF, returns the same matrix as A in this example, because A is upper triangular, and RCOND = 0.3, which is the estimated reciprocal condition number for the equilibrated matrix. The only output quantity that pertains to the original problem before equilibration is the solution matrix X. In this example, X is also the solution to the equilibrated problem because no column scaling was done, but if EQUED had returned 'C' or 'B' and the solution to the equilibrated system were desired, it could be computed from .

Iterative refinement in LAPACK uses all same-precision arithmetic, a recent innovation, because it was long believed that a successful algorithm would require the residual to be computed in double precision. The following example illustrates the results of iterative refinement; “Error Bounds”, discusses the error bounds computed in the course of the algorithm.

One possible use of iterative refinement is to smooth out numerical differences between floating-point number representations. For example, a result computed on a system, which has about 13 digits of accuracy, may be improved on IEEE system, which has about 15 digits of accuracy, through the 0(n²) process of iterative refinement, instead of the 0(n³) process of recomputing the solution.

Example 3-5. Hilbert matrix

The classic example of an ill-conditioned matrix is the Hilbert matrix, defined by . For example, the 5-by-5 Hilbert matrix is

Equation 3-50.

which has a condition number of . A rule of thumb (“Condition Estimation”) suggests almost 8 digits of accuracy in the solution are possible on systems, because and . If the matrix is factored using DGETRF and DGETRS is used to solve , where .

Error Bounds

In addition to performing iterative refinement on each column of the solution matrix, DGERFS and the other xxxRFS routines also compute error bounds for each column of the solution. These bounds are returned in the real arrays FERR (forward error) and BERR (backward error), both of length NRHS. For a computed solution x^ to the system of equations , the forward error bound ƒ is the following:

Equation 3-51.

and the backward error bound ω bounds the relative errors in each component of A and b in the perturbed equation 2 from “Use in Error Bounds”.

Equation 3-52.

In Example 3-5, the first column of the inverse of the Hilbert matrix of order 5 is computed by solving , and DGERFS computed error bounds of and on systems. This provides direct information about the solution; its relative error is at most 0(10^-8 ); therefore, the largest components of the solution should exhibit about 8 digits of accuracy, and the system for which this solution is an exact solution is within a factor of epsilon from the system whose solution was attempted, so the solution is as good as the data warrants.

The component-wise relative backward error bound (equation 3) is more restrictive than the classical backward error bound , because it assumes has the same sparsity structure as , because if is 0, so must be . The backward error for the solution x^ is computed from the equation

Equation 3-53.

where , and the forward error bound is computed from the equation

Equation 3-54.

where γ is the maximum number of nonzeros in any row of A. To avoid computing A ^-1 an approximation is used for , where g is the nonnegative vector in parentheses (see Arioli, et al., for details).

Subroutines to compute a matrix inverse are provided in LAPACK, but they are not used in the driver routines. The inverse routines sometimes use extra workspace and always require more operations than the solve routines. For example, if there is one right-hand side in the equation , where A is a square general matrix, the solves following the factorization require 2 n² operations, while inverting the matrix A requires 4/3n ³ operations, plus another 2n ² to multiply b by A^-1. The inverse must be computed only once, however, and the cost can be amortized over the number of right-hand sides. Because multiplying by the inverse may be more efficient than doing a triangular solve, the extra cost to compute the inverse may be overcome if the number of right-hand sides is large.

LAPACK provides four different orthogonal factorizations for each data type. For real data, they are as follows:

Each of these factorizations can be done regardless of whether m or n is larger, but the QR and QL factorizations are most often used when , and the LQ and RQ factorizations are most often used when .

The QR factorization of an m-by- n matrix A for has the form

Equation 3-57.

where Q is an m-by- m orthogonal matrix, and R is an n-by-n upper triangular matrix. If m > n, it is convenient to write the factorization as

Equation 3-58.

or simply

Equation 3-59.

where consists of the first n columns of Q, and consists of the remaining m- n columns. The LAPACK routine SGEQRF computes this factorization. See the LAPACK User's Guide for details.

Example 3-6. Orthogonal factorization

The result of calling DGEQRF with the matrix A equal to

Equation 3-60.

is a compact representation of Q and R consisting of

Equation 3-61.

The matrix R appears in the upper triangle of A explicitly:

Equation 3-62.

while the matrix Q is stored in a factored form where each is an elementary Householder transformation of the following form: . Each vector has and is stored below the diagonal in the column of A. Therefore,

Equation 3-63.

Equation 3-64.

Equation 3-65.

Each transformation is orthogonal (to machine precision) because τ is chosen to have the value , so that

Equation 3-66.

In the example of the last subsection, the elementary orthogonal matrices Q_i were expressed in terms of the scalar τ and the vector v that defines them, without multiplying out the expression. This was done to make the point that the elementary transformation is most often used in its factored form. This subsection describes one application in which multiplication by the orthogonal matrix Q is required. An alternative interface for this application is the LAPACK driver routine DGELS.

Given the QR factorization of a matrix A with m>n (Equation 3-57), if R has rank n, the solution to the linear least squares problem (Equation 3-55) is obtained by the following steps:

Equation 3-67.

The LAPACK routine DORMQR is used in step 1 to multiply the right-hand side matrix by the transpose of the orthogonal matrix Q, using Q in its factored form. The triangular system solution in step 2 can be done using the LAPACK routine DTRTRS.

Continuing the example of “Orthogonal Factorizations”, suppose the right-hand side vector is . Multiplying b by by using the LAPACK routine DORMQR, you get

Equation 3-68.

and after solving the triangular system with the 3-by-3 matrix R,

Equation 3-69.

The last m-n elements of can be used to compute the 2-norm of the residual, because . Here, .

The DORMxx routines described in the previous subsection are useful for operating with the orthogonal matrix Q in its factored form, but sometimes it is the matrix Q itself that is of interest.

For example, the first step in the computational procedure for the generalized eigenvalue problem is to find orthogonal matrices U and Z such that is upper Hessenberg and is upper triangular (see Golub and Van Loan for details). First, the matrix B is reduced to upper triangular form by the QR factorization, and A is overwritten by , using the subroutines of the previous section for multiplying by the orthogonal matrix.

Next, the updated A is reduced to Hessenberg form while preserving the triangular structure of B by applying Givens rotations to A and B, alternately from the left and the right. In order to obtain the orthogonal matrices U and Z that reduce the original problem, it is necessary to keep a copy of the matrix Q from and continually update it. This requires that Q be generated after the QR factorization.

The DORGxx routines generate the orthogonal matrix Q as a full matrix by expanding its factored form. Using the example of “Orthogonal Factorizations”, Q can be generated from its factored form by using the following Fortran code:

DO J = 1, N
   DO I = J+1, M
      Q(I,J) = A(I,J)
   END DO
END DO
CALL GORGQR(M,M,N,Q,LDQ,TAU,WORK,LWORK,INFO)

where the data from the QR factorization of A has been copied into a separate matrix Q because DORGQR overwrites its input data with the expanded matrix. The orthogonal matrix is returned in Q:

Equation 3-70.

The matrix Q⁽¹⁾, consisting of only the first n columns of Q, could be generated by specifying N, rather than M, as the second argument in the call to DORGQR.