Linear Equations

We use the standard notation for a system of simultaneous linear equations:

A x = b (2.4)

where A is the coefficient matrix, b is the right hand side, and x is the solution. In (2.4) A is assumed to be a square matrix of order n, but some of the individual routines allow A to be rectangular. If there are several right hand sides, we write

A X = B (2.5)

where the columns of B are the individual right hand sides, and the columns of X are the corresponding solutions. The basic task is to compute X, given A and B.

If A is upper or lower triangular, (2.4) can be solved by a straightforward process of backward or forward substitution. Otherwise, the solution is obtained after first factorizing A as a product of triangular matrices (and possibly also a diagonal matrix or permutation matrix).

The form of the factorization depends on the properties of the matrix A. LAPACK provides routines for the following types of matrices, based on the stated factorizations:

• general matrices (LU factorization with partial pivoting):
  
  A = PLU
  
  
  where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).

• general band matrices including tridiagonal matrices (LU factorization with partial pivoting): If A is m-by-n with kl subdiagonals and ku superdiagonals, the factorization is
  
  A = LU
  
  
  where L is a product of permutation and unit lower triangular matrices with kl subdiagonals, and U is upper triangular with kl+ku superdiagonals.
• symmetric and Hermitian positive definite matrices including band matrices (Cholesky factorization):
  
  $$A = U^{T}U \; \; {\rm or}\; \; A = LL^{T} \mbox{(in the symmetric case)}$$
  
  
  
  
  $$A = U^{H}U \; \; {\rm or}\; \; A = LL^{H} \mbox{(in the Hermitian case)}$$
  
  
  where U is an upper triangular matrix and L is lower triangular.
• symmetric and Hermitian positive definite tridiagonal matrices (L D L^T factorization):
  
  $$A = U D U^{T} \; \; {\rm or}\; \; A = L D L^{T} \mbox{(in the symmetric case)}$$
  
  
  
  
  $$A = U D U^{H} \; \; {\rm or}\; \; A = L D L^{H} \mbox{(in the Hermitian case)}$$
  
  
  where U is a unit upper bidiagonal matrix, L is unit lower bidiagonal, and D is diagonal.

• symmetric and Hermitian indefinite matrices (symmetric indefinite factorization):
  
  $$A = U D U^{T} \; \; {\rm or}\; \; A = L D L^{T} \mbox{(in the symmetric case)}$$
  
  
  
  
  $$A = U D U^{H} \; \; {\rm or}\; \; A = L D L^{H} \mbox{(in the Hermitian case)}$$
  
  
  where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with diagonal blocks of order 1 or 2.

The factorization for a general tridiagonal matrix is like that for a general band matrix with kl = 1 and ku = 1. The factorization for a symmetric positive definite band matrix with k superdiagonals (or subdiagonals) has the same form as for a symmetric positive definite matrix, but the factor U (or L) is a band matrix with k superdiagonals (subdiagonals). Band matrices use a compact band storage scheme described in section 5.3.3. LAPACK routines are also provided for symmetric matrices (whether positive definite or indefinite) using packed storage, as described in section 5.3.2.

While the primary use of a matrix factorization is to solve a system of equations, other related tasks are provided as well. Wherever possible, LAPACK provides routines to perform each of these tasks for each type of matrix and storage scheme (see Tables 2.7 and 2.8). The following list relates the tasks to the last 3 characters of the name of the corresponding computational routine:

xyyTRF:

factorize (obviously not needed for triangular matrices);

xyyTRS:

use the factorization (or the matrix A itself if it is triangular) to solve (2.5) by forward or backward substitution;

xyyCON:

estimate the reciprocal of the condition number $\kappa(A) = \Vert A\Vert . \Vert A^{-1} \Vert$ ; Higham's modification [63] of Hager's method [59] is used to estimate |A^-1|, except for symmetric positive definite tridiagonal matrices for which it is computed directly with comparable efficiency [61];

xyyRFS:

compute bounds on the error in the computed solution (returned by the xyyTRS routine), and refine the solution to reduce the backward error (see below);

xyyTRI:

use the factorization (or the matrix A itself if it is triangular) to compute A^-1 (not provided for band matrices, because the inverse does not in general preserve bandedness);

xyyEQU:

compute scaling factors to equilibrate A (not provided for tridiagonal, symmetric indefinite, or triangular matrices). These routines do not actually scale the matrices: auxiliary routines xLAQyy may be used for that purpose -- see the code of the driver routines xyySVX for sample usage.

Note that some of the above routines depend on the output of others:

xyyTRF:

may work on an equilibrated matrix produced by xyyEQU and xLAQyy, if yy is one of {GE, GB, PO, PP, PB};

xyyTRS:

requires the factorization returned by xyyTRF;

xyyCON:

requires the norm of the original matrix A, and the factorization returned by xyyTRF;

xyyRFS:

requires the original matrices A and B, the factorization returned by xyyTRF, and the solution X returned by xyyTRS;

xyyTRI:

requires the factorization returned by xyyTRF.

The RFS (``refine solution'') routines perform iterative refinement and compute backward and forward error bounds for the solution. Iterative refinement is done in the same precision as the input data. In particular, the residual is not computed with extra precision, as has been traditionally done. The benefit of this procedure is discussed in Section 4.4.

Table 2.7: Computational routines for linear equations

Type of matrix	Operation	Single precision		Double precision
and storage scheme		real	complex	real	complex
general	factorize	SGETRF	CGETRF	DGETRF	ZGETRF
	solve using factorization	SGETRS	CGETRS	DGETRS	ZGETRS
	estimate condition number	SGECON	CGECON	DGECON	ZGECON
	error bounds for solution	SGERFS	CGERFS	DGERFS	ZGERFS
	invert using factorization	SGETRI	CGETRI	DGETRI	ZGETRI
	equilibrate	SGEEQU	CGEEQU	DGEEQU	ZGEEQU
general	factorize	SGBTRF	CGBTRF	DGBTRF	ZGBTRF
band	solve using factorization	SGBTRS	CGBTRS	DGBTRS	ZGBTRS
	estimate condition number	SGBCON	CGBCON	DGBCON	ZGBCON
	error bounds for solution	SGBRFS	CGBRFS	DGBRFS	ZGBRFS
	equilibrate	SGBEQU	CGBEQU	DGBEQU	ZGBEQU
general	factorize	SGTTRF	CGTTRF	DGTTRF	ZGTTRF
tridiagonal	solve using factorization	SGTTRS	CGTTRS	DGTTRS	ZGTTRS
	estimate condition number	SGTCON	CGTCON	DGTCON	ZGTCON
	error bounds for solution	SGTRFS	CGTRFS	DGTRFS	ZGTRFS
symmetric/Hermitian	factorize	SPOTRF	CPOTRF	DPOTRF	ZPOTRF
positive definite	solve using factorization	SPOTRS	CPOTRS	DPOTRS	ZPOTRS
	estimate condition number	SPOCON	CPOCON	DPOCON	ZPOCON
	error bounds for solution	SPORFS	CPORFS	DPORFS	ZPORFS
	invert using factorization	SPOTRI	CPOTRI	DPOTRI	ZPOTRI
	equilibrate	SPOEQU	CPOEQU	DPOEQU	ZPOEQU
symmetric/Hermitian	factorize	SPPTRF	CPPTRF	DPPTRF	ZPPTRF
positive definite	solve using factorization	SPPTRS	CPPTRS	DPPTRS	ZPPTRS
(packed storage)	estimate condition number	SPPCON	CPPCON	DPPCON	ZPPCON
	error bounds for solution	SPPRFS	CPPRFS	DPPRFS	ZPPRFS
	invert using factorization	SPPTRI	CPPTRI	DPPTRI	ZPPTRI
	equilibrate	SPPEQU	CPPEQU	DPPEQU	ZPPEQU
symmetric/Hermitian	factorize	SPBTRF	CPBTRF	DPBTRF	ZPBTRF
positive definite	solve using factorization	SPBTRS	CPBTRS	DPBTRS	ZPBTRS
band	estimate condition number	SPBCON	CPBCON	DPBCON	ZPBCON
	error bounds for solution	SPBRFS	CPBRFS	DPBRFS	ZPBRFS
	equilibrate	SPBEQU	CPBEQU	DPBEQU	ZPBEQU
symmetric/Hermitian	factorize	SPTTRF	CPTTRF	DPTTRF	ZPTTRF
positive definite	solve using factorization	SPTTRS	CPTTRS	DPTTRS	ZPTTRS
tridiagonal	estimate condition number	SPTCON	CPTCON	DPTCON	ZPTCON
	error bounds for solution	SPTRFS	CPTRFS	DPTRFS	ZPTRFS

Table 2.8: Computational routines for linear equations (continued)

Type of matrix	Operation	Single precision		Double precision
and storage scheme		real	complex	real	complex
symmetric/Hermitian	factorize	SSYTRF	CHETRF	DSYTRF	ZHETRF
indefinite	solve using factorization	SSYTRS	CHETRS	DSYTRS	ZHETRS
	estimate condition number	SSYCON	CHECON	DSYCON	ZHECON
	error bounds for solution	SSYRFS	CHERFS	DSYRFS	ZHERFS
	invert using factorization	SSYTRI	CHETRI	DSYTRI	ZHETRI
complex symmetric	factorize		CSYTRF		ZSYTRF
	solve using factorization		CSYTRS		ZSYTRS
	estimate condition number		CSYCON		ZSYCON
	error bounds for solution		CSYRFS		ZSYRFS
	invert using factorization		CSYTRI		ZSYTRI
symmetric/Hermitian	factorize	SSPTRF	CHPTRF	DSPTRF	ZHPTRF
indefinite	solve using factorization	SSPTRS	CHPTRS	DSPTRS	ZHPTRS
(packed storage)	estimate condition number	SSPCON	CHPCON	DSPCON	ZHPCON
	error bounds for solution	SSPRFS	CHPRFS	DSPRFS	ZHPRFS
	invert using factorization	SSPTRI	CHPTRI	DSPTRI	ZHPTRI
complex symmetric	factorize		CSPTRF		ZSPTRF
(packed storage)	solve using factorization		CSPTRS		ZSPTRS
	estimate condition number		CSPCON		ZSPCON
	error bounds for solution		CSPRFS		ZSPRFS
	invert using factorization		CSPTRI		ZSPTRI
triangular	solve	STRTRS	CTRTRS	DTRTRS	ZTRTRS
	estimate condition number	STRCON	CTRCON	DTRCON	ZTRCON
	error bounds for solution	STRRFS	CTRRFS	DTRRFS	ZTRRFS
	invert	STRTRI	CTRTRI	DTRTRI	ZTRTRI
triangular	solve	STPTRS	CTPTRS	DTPTRS	ZTPTRS
(packed storage)	estimate condition number	STPCON	CTPCON	DTPCON	ZTPCON
	error bounds for solution	STPRFS	CTPRFS	DTPRFS	ZTPRFS
	invert	STPTRI	CTPTRI	DTPTRI	ZTPTRI
triangular	solve	STBTRS	CTBTRS	DTBTRS	ZTBTRS
band	estimate condition number	STBCON	CTBCON	DTBCON	ZTBCON
	error bounds for solution	STBRFS	CTBRFS	DTBRFS	ZTBRFS