diff options
Diffstat (limited to 'gsl-1.9/doc/linalg.texi')
| -rw-r--r-- | gsl-1.9/doc/linalg.texi | 1117 |
1 files changed, 1117 insertions, 0 deletions
diff --git a/gsl-1.9/doc/linalg.texi b/gsl-1.9/doc/linalg.texi new file mode 100644 index 0000000..839f4d5 --- /dev/null +++ b/gsl-1.9/doc/linalg.texi @@ -0,0 +1,1117 @@ +@cindex linear algebra +@cindex solution of linear systems, Ax=b +@cindex matrix factorization +@cindex factorization of matrices + +This chapter describes functions for solving linear systems. The +library provides linear algebra operations which operate directly on +the @code{gsl_vector} and @code{gsl_matrix} objects. These routines +use the standard algorithms from Golub & Van Loan's @cite{Matrix +Computations}. + +@cindex LAPACK, recommended for linear algebra +When dealing with very large systems the routines found in @sc{lapack} +should be considered. These support specialized data representations +and other optimizations. + +The functions described in this chapter are declared in the header file +@file{gsl_linalg.h}. + + +@menu +* LU Decomposition:: +* QR Decomposition:: +* QR Decomposition with Column Pivoting:: +* Singular Value Decomposition:: +* Cholesky Decomposition:: +* Tridiagonal Decomposition of Real Symmetric Matrices:: +* Tridiagonal Decomposition of Hermitian Matrices:: +* Hessenberg Decomposition of Real Matrices:: +* Bidiagonalization:: +* Householder Transformations:: +* Householder solver for linear systems:: +* Tridiagonal Systems:: +* Balancing:: +* Linear Algebra Examples:: +* Linear Algebra References and Further Reading:: +@end menu + +@node LU Decomposition +@section LU Decomposition +@cindex LU decomposition + +A general square matrix @math{A} has an @math{LU} decomposition into +upper and lower triangular matrices, +@tex +\beforedisplay +$$ +P A = L U +$$ +\afterdisplay +@end tex +@ifinfo + +@example +P A = L U +@end example + +@end ifinfo +@noindent +where @math{P} is a permutation matrix, @math{L} is unit lower +triangular matrix and @math{U} is upper triangular matrix. For square +matrices this decomposition can be used to convert the linear system +@math{A x = b} into a pair of triangular systems (@math{L y = P b}, +@math{U x = y}), which can be solved by forward and back-substitution. +Note that the @math{LU} decomposition is valid for singular matrices. + +@deftypefun int gsl_linalg_LU_decomp (gsl_matrix * @var{A}, gsl_permutation * @var{p}, int * @var{signum}) +@deftypefunx int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * @var{A}, gsl_permutation * @var{p}, int * @var{signum}) +These functions factorize the square matrix @var{A} into the @math{LU} +decomposition @math{PA = LU}. On output the diagonal and upper +triangular part of the input matrix @var{A} contain the matrix +@math{U}. The lower triangular part of the input matrix (excluding the +diagonal) contains @math{L}. The diagonal elements of @math{L} are +unity, and are not stored. + +The permutation matrix @math{P} is encoded in the permutation +@var{p}. The @math{j}-th column of the matrix @math{P} is given by the +@math{k}-th column of the identity matrix, where @math{k = p_j} the +@math{j}-th element of the permutation vector. The sign of the +permutation is given by @var{signum}. It has the value @math{(-1)^n}, +where @math{n} is the number of interchanges in the permutation. + +The algorithm used in the decomposition is Gaussian Elimination with +partial pivoting (Golub & Van Loan, @cite{Matrix Computations}, +Algorithm 3.4.1). +@end deftypefun + +@cindex linear systems, solution of +@deftypefun int gsl_linalg_LU_solve (const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +@deftypefunx int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector_complex * @var{b}, gsl_vector_complex * @var{x}) +These functions solve the square system @math{A x = b} using the @math{LU} +decomposition of @math{A} into (@var{LU}, @var{p}) given by +@code{gsl_linalg_LU_decomp} or @code{gsl_linalg_complex_LU_decomp}. +@end deftypefun + +@deftypefun int gsl_linalg_LU_svx (const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, gsl_vector * @var{x}) +@deftypefunx int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, gsl_vector_complex * @var{x}) +These functions solve the square system @math{A x = b} in-place using the +@math{LU} decomposition of @math{A} into (@var{LU},@var{p}). On input +@var{x} should contain the right-hand side @math{b}, which is replaced +by the solution on output. +@end deftypefun + +@cindex refinement of solutions in linear systems +@cindex iterative refinement of solutions in linear systems +@cindex linear systems, refinement of solutions +@deftypefun int gsl_linalg_LU_refine (const gsl_matrix * @var{A}, const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x}, gsl_vector * @var{residual}) +@deftypefunx int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector_complex * @var{b}, gsl_vector_complex * @var{x}, gsl_vector_complex * @var{residual}) +These functions apply an iterative improvement to @var{x}, the solution +of @math{A x = b}, using the @math{LU} decomposition of @math{A} into +(@var{LU},@var{p}). The initial residual @math{r = A x - b} is also +computed and stored in @var{residual}. +@end deftypefun + +@cindex inverse of a matrix, by LU decomposition +@cindex matrix inverse +@deftypefun int gsl_linalg_LU_invert (const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, gsl_matrix * @var{inverse}) +@deftypefunx int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, gsl_matrix_complex * @var{inverse}) +These functions compute the inverse of a matrix @math{A} from its +@math{LU} decomposition (@var{LU},@var{p}), storing the result in the +matrix @var{inverse}. The inverse is computed by solving the system +@math{A x = b} for each column of the identity matrix. It is preferable +to avoid direct use of the inverse whenever possible, as the linear +solver functions can obtain the same result more efficiently and +reliably (consult any introductory textbook on numerical linear algebra +for details). +@end deftypefun + +@cindex determinant of a matrix, by LU decomposition +@cindex matrix determinant +@deftypefun double gsl_linalg_LU_det (gsl_matrix * @var{LU}, int @var{signum}) +@deftypefunx gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * @var{LU}, int @var{signum}) +These functions compute the determinant of a matrix @math{A} from its +@math{LU} decomposition, @var{LU}. The determinant is computed as the +product of the diagonal elements of @math{U} and the sign of the row +permutation @var{signum}. +@end deftypefun + +@cindex logarithm of the determinant of a matrix +@deftypefun double gsl_linalg_LU_lndet (gsl_matrix * @var{LU}) +@deftypefunx double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * @var{LU}) +These functions compute the logarithm of the absolute value of the +determinant of a matrix @math{A}, @math{\ln|\det(A)|}, from its @math{LU} +decomposition, @var{LU}. This function may be useful if the direct +computation of the determinant would overflow or underflow. +@end deftypefun + +@cindex sign of the determinant of a matrix +@deftypefun int gsl_linalg_LU_sgndet (gsl_matrix * @var{LU}, int @var{signum}) +@deftypefunx gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * @var{LU}, int @var{signum}) +These functions compute the sign or phase factor of the determinant of a +matrix @math{A}, @math{\det(A)/|\det(A)|}, from its @math{LU} decomposition, +@var{LU}. +@end deftypefun + +@node QR Decomposition +@section QR Decomposition +@cindex QR decomposition + +A general rectangular @math{M}-by-@math{N} matrix @math{A} has a +@math{QR} decomposition into the product of an orthogonal +@math{M}-by-@math{M} square matrix @math{Q} (where @math{Q^T Q = I}) and +an @math{M}-by-@math{N} right-triangular matrix @math{R}, +@tex +\beforedisplay +$$ +A = Q R +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = Q R +@end example + +@end ifinfo +@noindent +This decomposition can be used to convert the linear system @math{A x = +b} into the triangular system @math{R x = Q^T b}, which can be solved by +back-substitution. Another use of the @math{QR} decomposition is to +compute an orthonormal basis for a set of vectors. The first @math{N} +columns of @math{Q} form an orthonormal basis for the range of @math{A}, +@math{ran(A)}, when @math{A} has full column rank. + +@deftypefun int gsl_linalg_QR_decomp (gsl_matrix * @var{A}, gsl_vector * @var{tau}) +This function factorizes the @math{M}-by-@math{N} matrix @var{A} into +the @math{QR} decomposition @math{A = Q R}. On output the diagonal and +upper triangular part of the input matrix contain the matrix +@math{R}. The vector @var{tau} and the columns of the lower triangular +part of the matrix @var{A} contain the Householder coefficients and +Householder vectors which encode the orthogonal matrix @var{Q}. The +vector @var{tau} must be of length @math{k=\min(M,N)}. The matrix +@math{Q} is related to these components by, @math{Q = Q_k ... Q_2 Q_1} +where @math{Q_i = I - \tau_i v_i v_i^T} and @math{v_i} is the +Householder vector @math{v_i = +(0,...,1,A(i+1,i),A(i+2,i),...,A(m,i))}. This is the same storage scheme +as used by @sc{lapack}. + +The algorithm used to perform the decomposition is Householder QR (Golub +& Van Loan, @cite{Matrix Computations}, Algorithm 5.2.1). +@end deftypefun + +@deftypefun int gsl_linalg_QR_solve (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the square system @math{A x = b} using the @math{QR} +decomposition of @math{A} into (@var{QR}, @var{tau}) given by +@code{gsl_linalg_QR_decomp}. The least-squares solution for rectangular systems can +be found using @code{gsl_linalg_QR_lssolve}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_svx (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_vector * @var{x}) +This function solves the square system @math{A x = b} in-place using the +@math{QR} decomposition of @math{A} into (@var{QR},@var{tau}) given by +@code{gsl_linalg_QR_decomp}. On input @var{x} should contain the +right-hand side @math{b}, which is replaced by the solution on output. +@end deftypefun + +@deftypefun int gsl_linalg_QR_lssolve (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, const gsl_vector * @var{b}, gsl_vector * @var{x}, gsl_vector * @var{residual}) +This function finds the least squares solution to the overdetermined +system @math{A x = b} where the matrix @var{A} has more rows than +columns. The least squares solution minimizes the Euclidean norm of the +residual, @math{||Ax - b||}.The routine uses the @math{QR} decomposition +of @math{A} into (@var{QR}, @var{tau}) given by +@code{gsl_linalg_QR_decomp}. The solution is returned in @var{x}. The +residual is computed as a by-product and stored in @var{residual}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_QTvec (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_vector * @var{v}) +This function applies the matrix @math{Q^T} encoded in the decomposition +(@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q^T +v} in @var{v}. The matrix multiplication is carried out directly using +the encoding of the Householder vectors without needing to form the full +matrix @math{Q^T}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_Qvec (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_vector * @var{v}) +This function applies the matrix @math{Q} encoded in the decomposition +(@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q +v} in @var{v}. The matrix multiplication is carried out directly using +the encoding of the Householder vectors without needing to form the full +matrix @math{Q}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_Rsolve (const gsl_matrix * @var{QR}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the triangular system @math{R x = b} for +@var{x}. It may be useful if the product @math{b' = Q^T b} has already +been computed using @code{gsl_linalg_QR_QTvec}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_Rsvx (const gsl_matrix * @var{QR}, gsl_vector * @var{x}) +This function solves the triangular system @math{R x = b} for @var{x} +in-place. On input @var{x} should contain the right-hand side @math{b} +and is replaced by the solution on output. This function may be useful if +the product @math{b' = Q^T b} has already been computed using +@code{gsl_linalg_QR_QTvec}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_unpack (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_matrix * @var{Q}, gsl_matrix * @var{R}) +This function unpacks the encoded @math{QR} decomposition +(@var{QR},@var{tau}) into the matrices @var{Q} and @var{R}, where +@var{Q} is @math{M}-by-@math{M} and @var{R} is @math{M}-by-@math{N}. +@end deftypefun + +@deftypefun int gsl_linalg_QR_QRsolve (gsl_matrix * @var{Q}, gsl_matrix * @var{R}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the system @math{R x = Q^T b} for @var{x}. It can +be used when the @math{QR} decomposition of a matrix is available in +unpacked form as (@var{Q}, @var{R}). +@end deftypefun + +@deftypefun int gsl_linalg_QR_update (gsl_matrix * @var{Q}, gsl_matrix * @var{R}, gsl_vector * @var{w}, const gsl_vector * @var{v}) +This function performs a rank-1 update @math{w v^T} of the @math{QR} +decomposition (@var{Q}, @var{R}). The update is given by @math{Q'R' = Q +R + w v^T} where the output matrices @math{Q'} and @math{R'} are also +orthogonal and right triangular. Note that @var{w} is destroyed by the +update. +@end deftypefun + +@deftypefun int gsl_linalg_R_solve (const gsl_matrix * @var{R}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the triangular system @math{R x = b} for the +@math{N}-by-@math{N} matrix @var{R}. +@end deftypefun + +@deftypefun int gsl_linalg_R_svx (const gsl_matrix * @var{R}, gsl_vector * @var{x}) +This function solves the triangular system @math{R x = b} in-place. On +input @var{x} should contain the right-hand side @math{b}, which is +replaced by the solution on output. +@end deftypefun + +@node QR Decomposition with Column Pivoting +@section QR Decomposition with Column Pivoting +@cindex QR decomposition with column pivoting + +The @math{QR} decomposition can be extended to the rank deficient case +by introducing a column permutation @math{P}, +@tex +\beforedisplay +$$ +A P = Q R +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A P = Q R +@end example + +@end ifinfo +@noindent +The first @math{r} columns of @math{Q} form an orthonormal basis +for the range of @math{A} for a matrix with column rank @math{r}. This +decomposition can also be used to convert the linear system @math{A x = +b} into the triangular system @math{R y = Q^T b, x = P y}, which can be +solved by back-substitution and permutation. We denote the @math{QR} +decomposition with column pivoting by @math{QRP^T} since @math{A = Q R +P^T}. + +@deftypefun int gsl_linalg_QRPT_decomp (gsl_matrix * @var{A}, gsl_vector * @var{tau}, gsl_permutation * @var{p}, int * @var{signum}, gsl_vector * @var{norm}) +This function factorizes the @math{M}-by-@math{N} matrix @var{A} into +the @math{QRP^T} decomposition @math{A = Q R P^T}. On output the +diagonal and upper triangular part of the input matrix contain the +matrix @math{R}. The permutation matrix @math{P} is stored in the +permutation @var{p}. The sign of the permutation is given by +@var{signum}. It has the value @math{(-1)^n}, where @math{n} is the +number of interchanges in the permutation. The vector @var{tau} and the +columns of the lower triangular part of the matrix @var{A} contain the +Householder coefficients and vectors which encode the orthogonal matrix +@var{Q}. The vector @var{tau} must be of length @math{k=\min(M,N)}. The +matrix @math{Q} is related to these components by, @math{Q = Q_k ... Q_2 +Q_1} where @math{Q_i = I - \tau_i v_i v_i^T} and @math{v_i} is the +Householder vector @math{v_i = +(0,...,1,A(i+1,i),A(i+2,i),...,A(m,i))}. This is the same storage scheme +as used by @sc{lapack}. The vector @var{norm} is a workspace of length +@var{N} used for column pivoting. + +The algorithm used to perform the decomposition is Householder QR with +column pivoting (Golub & Van Loan, @cite{Matrix Computations}, Algorithm +5.4.1). +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_decomp2 (const gsl_matrix * @var{A}, gsl_matrix * @var{q}, gsl_matrix * @var{r}, gsl_vector * @var{tau}, gsl_permutation * @var{p}, int * @var{signum}, gsl_vector * @var{norm}) +This function factorizes the matrix @var{A} into the decomposition +@math{A = Q R P^T} without modifying @var{A} itself and storing the +output in the separate matrices @var{q} and @var{r}. +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_solve (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the square system @math{A x = b} using the @math{QRP^T} +decomposition of @math{A} into (@var{QR}, @var{tau}, @var{p}) given by +@code{gsl_linalg_QRPT_decomp}. +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_svx (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, const gsl_permutation * @var{p}, gsl_vector * @var{x}) +This function solves the square system @math{A x = b} in-place using the +@math{QRP^T} decomposition of @math{A} into +(@var{QR},@var{tau},@var{p}). On input @var{x} should contain the +right-hand side @math{b}, which is replaced by the solution on output. +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_QRsolve (const gsl_matrix * @var{Q}, const gsl_matrix * @var{R}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the square system @math{R P^T x = Q^T b} for +@var{x}. It can be used when the @math{QR} decomposition of a matrix is +available in unpacked form as (@var{Q}, @var{R}). +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_update (gsl_matrix * @var{Q}, gsl_matrix * @var{R}, const gsl_permutation * @var{p}, gsl_vector * @var{u}, const gsl_vector * @var{v}) +This function performs a rank-1 update @math{w v^T} of the @math{QRP^T} +decomposition (@var{Q}, @var{R}, @var{p}). The update is given by +@math{Q'R' = Q R + w v^T} where the output matrices @math{Q'} and +@math{R'} are also orthogonal and right triangular. Note that @var{w} is +destroyed by the update. The permutation @var{p} is not changed. +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_Rsolve (const gsl_matrix * @var{QR}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the triangular system @math{R P^T x = b} for the +@math{N}-by-@math{N} matrix @math{R} contained in @var{QR}. +@end deftypefun + +@deftypefun int gsl_linalg_QRPT_Rsvx (const gsl_matrix * @var{QR}, const gsl_permutation * @var{p}, gsl_vector * @var{x}) +This function solves the triangular system @math{R P^T x = b} in-place +for the @math{N}-by-@math{N} matrix @math{R} contained in @var{QR}. On +input @var{x} should contain the right-hand side @math{b}, which is +replaced by the solution on output. +@end deftypefun + +@node Singular Value Decomposition +@section Singular Value Decomposition +@cindex SVD +@cindex singular value decomposition + +A general rectangular @math{M}-by-@math{N} matrix @math{A} has a +singular value decomposition (@sc{svd}) into the product of an +@math{M}-by-@math{N} orthogonal matrix @math{U}, an @math{N}-by-@math{N} +diagonal matrix of singular values @math{S} and the transpose of an +@math{N}-by-@math{N} orthogonal square matrix @math{V}, +@tex +\beforedisplay +$$ +A = U S V^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = U S V^T +@end example + +@end ifinfo +@noindent +The singular values +@c{$\sigma_i = S_{ii}$} +@math{\sigma_i = S_@{ii@}} are all non-negative and are +generally chosen to form a non-increasing sequence +@c{$\sigma_1 \ge \sigma_2 \ge ... \ge \sigma_N \ge 0$} +@math{\sigma_1 >= \sigma_2 >= ... >= \sigma_N >= 0}. + +The singular value decomposition of a matrix has many practical uses. +The condition number of the matrix is given by the ratio of the largest +singular value to the smallest singular value. The presence of a zero +singular value indicates that the matrix is singular. The number of +non-zero singular values indicates the rank of the matrix. In practice +singular value decomposition of a rank-deficient matrix will not produce +exact zeroes for singular values, due to finite numerical +precision. Small singular values should be edited by choosing a suitable +tolerance. + +For a rank-deficient matrix, the null space of @math{A} is given by +the columns of @math{V} corresponding to the zero singular values. +Similarly, the range of @math{A} is given by columns of @math{U} +corresponding to the non-zero singular values. + +@deftypefun int gsl_linalg_SV_decomp (gsl_matrix * @var{A}, gsl_matrix * @var{V}, gsl_vector * @var{S}, gsl_vector * @var{work}) +This function factorizes the @math{M}-by-@math{N} matrix @var{A} into +the singular value decomposition @math{A = U S V^T} for @c{$M \ge N$} +@math{M >= N}. On output the matrix @var{A} is replaced by +@math{U}. The diagonal elements of the singular value matrix @math{S} +are stored in the vector @var{S}. The singular values are non-negative +and form a non-increasing sequence from @math{S_1} to @math{S_N}. The +matrix @var{V} contains the elements of @math{V} in untransposed +form. To form the product @math{U S V^T} it is necessary to take the +transpose of @var{V}. A workspace of length @var{N} is required in +@var{work}. + +This routine uses the Golub-Reinsch SVD algorithm. +@end deftypefun + +@deftypefun int gsl_linalg_SV_decomp_mod (gsl_matrix * @var{A}, gsl_matrix * @var{X}, gsl_matrix * @var{V}, gsl_vector * @var{S}, gsl_vector * @var{work}) +This function computes the SVD using the modified Golub-Reinsch +algorithm, which is faster for @c{$M \gg N$} +@math{M>>N}. It requires the vector @var{work} of length @var{N} and the +@math{N}-by-@math{N} matrix @var{X} as additional working space. +@end deftypefun + +@deftypefun int gsl_linalg_SV_decomp_jacobi (gsl_matrix * @var{A}, gsl_matrix * @var{V}, gsl_vector * @var{S}) +This function computes the SVD of the @math{M}-by-@math{N} matrix @var{A} +using one-sided Jacobi orthogonalization for @c{$M \ge N$} +@math{M >= N}. The Jacobi method can compute singular values to higher +relative accuracy than Golub-Reinsch algorithms (see references for +details). +@end deftypefun + +@deftypefun int gsl_linalg_SV_solve (gsl_matrix * @var{U}, gsl_matrix * @var{V}, gsl_vector * @var{S}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the system @math{A x = b} using the singular value +decomposition (@var{U}, @var{S}, @var{V}) of @math{A} given by +@code{gsl_linalg_SV_decomp}. + +Only non-zero singular values are used in computing the solution. The +parts of the solution corresponding to singular values of zero are +ignored. Other singular values can be edited out by setting them to +zero before calling this function. + +In the over-determined case where @var{A} has more rows than columns the +system is solved in the least squares sense, returning the solution +@var{x} which minimizes @math{||A x - b||_2}. +@end deftypefun + +@node Cholesky Decomposition +@section Cholesky Decomposition +@cindex Cholesky decomposition +@cindex square root of a matrix, Cholesky decomposition +@cindex matrix square root, Cholesky decomposition + +A symmetric, positive definite square matrix @math{A} has a Cholesky +decomposition into a product of a lower triangular matrix @math{L} and +its transpose @math{L^T}, +@tex +\beforedisplay +$$ +A = L L^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = L L^T +@end example + +@end ifinfo +@noindent +This is sometimes referred to as taking the square-root of a matrix. The +Cholesky decomposition can only be carried out when all the eigenvalues +of the matrix are positive. This decomposition can be used to convert +the linear system @math{A x = b} into a pair of triangular systems +(@math{L y = b}, @math{L^T x = y}), which can be solved by forward and +back-substitution. + +@deftypefun int gsl_linalg_cholesky_decomp (gsl_matrix * @var{A}) +This function factorizes the positive-definite symmetric square matrix +@var{A} into the Cholesky decomposition @math{A = L L^T}. On input +only the diagonal and lower-triangular part of the matrix @var{A} are +needed. On output the diagonal and lower triangular part of the input +matrix @var{A} contain the matrix @math{L}. The upper triangular part +of the input matrix contains @math{L^T}, the diagonal terms being +identical for both @math{L} and @math{L^T}. If the matrix is not +positive-definite then the decomposition will fail, returning the +error code @code{GSL_EDOM}. +@end deftypefun + +@deftypefun int gsl_linalg_cholesky_solve (const gsl_matrix * @var{cholesky}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the system @math{A x = b} using the Cholesky +decomposition of @math{A} into the matrix @var{cholesky} given by +@code{gsl_linalg_cholesky_decomp}. +@end deftypefun + +@deftypefun int gsl_linalg_cholesky_svx (const gsl_matrix * @var{cholesky}, gsl_vector * @var{x}) +This function solves the system @math{A x = b} in-place using the +Cholesky decomposition of @math{A} into the matrix @var{cholesky} given +by @code{gsl_linalg_cholesky_decomp}. On input @var{x} should contain +the right-hand side @math{b}, which is replaced by the solution on +output. +@end deftypefun + +@node Tridiagonal Decomposition of Real Symmetric Matrices +@section Tridiagonal Decomposition of Real Symmetric Matrices +@cindex tridiagonal decomposition + +A symmetric matrix @math{A} can be factorized by similarity +transformations into the form, +@tex +\beforedisplay +$$ +A = Q T Q^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = Q T Q^T +@end example + +@end ifinfo +@noindent +where @math{Q} is an orthogonal matrix and @math{T} is a symmetric +tridiagonal matrix. + +@deftypefun int gsl_linalg_symmtd_decomp (gsl_matrix * @var{A}, gsl_vector * @var{tau}) +This function factorizes the symmetric square matrix @var{A} into the +symmetric tridiagonal decomposition @math{Q T Q^T}. On output the +diagonal and subdiagonal part of the input matrix @var{A} contain the +tridiagonal matrix @math{T}. The remaining lower triangular part of the +input matrix contains the Householder vectors which, together with the +Householder coefficients @var{tau}, encode the orthogonal matrix +@math{Q}. This storage scheme is the same as used by @sc{lapack}. The +upper triangular part of @var{A} is not referenced. +@end deftypefun + +@deftypefun int gsl_linalg_symmtd_unpack (const gsl_matrix * @var{A}, const gsl_vector * @var{tau}, gsl_matrix * @var{Q}, gsl_vector * @var{diag}, gsl_vector * @var{subdiag}) +This function unpacks the encoded symmetric tridiagonal decomposition +(@var{A}, @var{tau}) obtained from @code{gsl_linalg_symmtd_decomp} into +the orthogonal matrix @var{Q}, the vector of diagonal elements @var{diag} +and the vector of subdiagonal elements @var{subdiag}. +@end deftypefun + +@deftypefun int gsl_linalg_symmtd_unpack_T (const gsl_matrix * @var{A}, gsl_vector * @var{diag}, gsl_vector * @var{subdiag}) +This function unpacks the diagonal and subdiagonal of the encoded +symmetric tridiagonal decomposition (@var{A}, @var{tau}) obtained from +@code{gsl_linalg_symmtd_decomp} into the vectors @var{diag} and @var{subdiag}. +@end deftypefun + +@node Tridiagonal Decomposition of Hermitian Matrices +@section Tridiagonal Decomposition of Hermitian Matrices +@cindex tridiagonal decomposition + +A hermitian matrix @math{A} can be factorized by similarity +transformations into the form, +@tex +\beforedisplay +$$ +A = U T U^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = U T U^T +@end example + +@end ifinfo +@noindent +where @math{U} is a unitary matrix and @math{T} is a real symmetric +tridiagonal matrix. + + +@deftypefun int gsl_linalg_hermtd_decomp (gsl_matrix_complex * @var{A}, gsl_vector_complex * @var{tau}) +This function factorizes the hermitian matrix @var{A} into the symmetric +tridiagonal decomposition @math{U T U^T}. On output the real parts of +the diagonal and subdiagonal part of the input matrix @var{A} contain +the tridiagonal matrix @math{T}. The remaining lower triangular part of +the input matrix contains the Householder vectors which, together with +the Householder coefficients @var{tau}, encode the orthogonal matrix +@math{Q}. This storage scheme is the same as used by @sc{lapack}. The +upper triangular part of @var{A} and imaginary parts of the diagonal are +not referenced. +@end deftypefun + +@deftypefun int gsl_linalg_hermtd_unpack (const gsl_matrix_complex * @var{A}, const gsl_vector_complex * @var{tau}, gsl_matrix_complex * @var{Q}, gsl_vector * @var{diag}, gsl_vector * @var{subdiag}) +This function unpacks the encoded tridiagonal decomposition (@var{A}, +@var{tau}) obtained from @code{gsl_linalg_hermtd_decomp} into the +unitary matrix @var{U}, the real vector of diagonal elements @var{diag} and +the real vector of subdiagonal elements @var{subdiag}. +@end deftypefun + +@deftypefun int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * @var{A}, gsl_vector * @var{diag}, gsl_vector * @var{subdiag}) +This function unpacks the diagonal and subdiagonal of the encoded +tridiagonal decomposition (@var{A}, @var{tau}) obtained from the +@code{gsl_linalg_hermtd_decomp} into the real vectors +@var{diag} and @var{subdiag}. +@end deftypefun + +@node Hessenberg Decomposition of Real Matrices +@section Hessenberg Decomposition of Real Matrices +@cindex hessenberg decomposition + +A general matrix @math{A} can be decomposed by orthogonal +similarity transformations into the form +@tex +\beforedisplay +$$ +A = U H U^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = U H U^T +@end example + +@end ifinfo +where @math{U} is orthogonal and @math{H} is an upper Hessenberg matrix, +meaning that it has zeros below the first subdiagonal. The +Hessenberg reduction is the first step in the Schur decomposition +for the nonsymmetric eigenvalue problem, but has applications in +other areas as well. + +@deftypefun int gsl_linalg_hessenberg (gsl_matrix * @var{A}, gsl_vector * @var{tau}) +This function computes the Hessenberg decomposition of the matrix +@var{A} by applying the similarity transformation @math{H = U^T A U}. +On output, @math{H} is stored in the upper portion of @var{A}. The +information required to construct the matrix @math{U} is stored in +the lower triangular portion of @var{A}. @math{U} is a product +of @math{N - 2} Householder matrices. The Householder vectors +are stored in the lower portion of @var{A} (below the subdiagonal) +and the Householder coefficients are stored in the vector @var{tau}. +@var{tau} must be of length @var{N}. +@end deftypefun + +@deftypefun int gsl_linalg_hessenberg_unpack (gsl_matrix * @var{H}, gsl_vector * @var{tau}, gsl_matrix * @var{U}) +This function constructs the orthogonal matrix @math{U} from the +information stored in the Hessenberg matrix @var{H} along with the +vector @var{tau}. @var{H} and @var{tau} are outputs from +@code{gsl_linalg_hessenberg}. +@end deftypefun + +@deftypefun int gsl_linalg_hessenberg_unpack_accum (gsl_matrix * @var{H}, gsl_vector * @var{tau}, gsl_matrix * @var{V}) +This function is similar to @code{gsl_linalg_hessenberg_unpack}, except +it accumulates the matrix @var{U} into @var{V}, so that @math{V' = VU}. +The matrix @var{V} must be initialized prior to calling this function. +Setting @var{V} to the identity matrix provides the same result as +@code{gsl_linalg_hessenberg_unpack}. If @var{H} is order @var{N}, then +@var{V} must have @var{N} columns but may have any number of rows. +@end deftypefun + +@deftypefun void gsl_linalg_hessenberg_set_zero (gsl_matrix * @var{H}) +This function sets the lower triangular portion of @var{H}, below +the subdiagonal, to zero. It is useful for clearing out the +Householder vectors after calling @code{gsl_linalg_hessenberg}. +@end deftypefun + +@node Bidiagonalization +@section Bidiagonalization +@cindex bidiagonalization of real matrices + +A general matrix @math{A} can be factorized by similarity +transformations into the form, +@tex +\beforedisplay +$$ +A = U B V^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = U B V^T +@end example + +@end ifinfo +@noindent +where @math{U} and @math{V} are orthogonal matrices and @math{B} is a +@math{N}-by-@math{N} bidiagonal matrix with non-zero entries only on the +diagonal and superdiagonal. The size of @var{U} is @math{M}-by-@math{N} +and the size of @var{V} is @math{N}-by-@math{N}. + +@deftypefun int gsl_linalg_bidiag_decomp (gsl_matrix * @var{A}, gsl_vector * @var{tau_U}, gsl_vector * @var{tau_V}) +This function factorizes the @math{M}-by-@math{N} matrix @var{A} into +bidiagonal form @math{U B V^T}. The diagonal and superdiagonal of the +matrix @math{B} are stored in the diagonal and superdiagonal of @var{A}. +The orthogonal matrices @math{U} and @var{V} are stored as compressed +Householder vectors in the remaining elements of @var{A}. The +Householder coefficients are stored in the vectors @var{tau_U} and +@var{tau_V}. The length of @var{tau_U} must equal the number of +elements in the diagonal of @var{A} and the length of @var{tau_V} should +be one element shorter. +@end deftypefun + +@deftypefun int gsl_linalg_bidiag_unpack (const gsl_matrix * @var{A}, const gsl_vector * @var{tau_U}, gsl_matrix * @var{U}, const gsl_vector * @var{tau_V}, gsl_matrix * @var{V}, gsl_vector * @var{diag}, gsl_vector * @var{superdiag}) +This function unpacks the bidiagonal decomposition of @var{A} given by +@code{gsl_linalg_bidiag_decomp}, (@var{A}, @var{tau_U}, @var{tau_V}) +into the separate orthogonal matrices @var{U}, @var{V} and the diagonal +vector @var{diag} and superdiagonal @var{superdiag}. Note that @var{U} +is stored as a compact @math{M}-by-@math{N} orthogonal matrix satisfying +@math{U^T U = I} for efficiency. +@end deftypefun + +@deftypefun int gsl_linalg_bidiag_unpack2 (gsl_matrix * @var{A}, gsl_vector * @var{tau_U}, gsl_vector * @var{tau_V}, gsl_matrix * @var{V}) +This function unpacks the bidiagonal decomposition of @var{A} given by +@code{gsl_linalg_bidiag_decomp}, (@var{A}, @var{tau_U}, @var{tau_V}) +into the separate orthogonal matrices @var{U}, @var{V} and the diagonal +vector @var{diag} and superdiagonal @var{superdiag}. The matrix @var{U} +is stored in-place in @var{A}. +@end deftypefun + +@deftypefun int gsl_linalg_bidiag_unpack_B (const gsl_matrix * @var{A}, gsl_vector * @var{diag}, gsl_vector * @var{superdiag}) +This function unpacks the diagonal and superdiagonal of the bidiagonal +decomposition of @var{A} given by @code{gsl_linalg_bidiag_decomp}, into +the diagonal vector @var{diag} and superdiagonal vector @var{superdiag}. +@end deftypefun + +@node Householder Transformations +@section Householder Transformations +@cindex Householder matrix +@cindex Householder transformation +@cindex transformation, Householder + +A Householder transformation is a rank-1 modification of the identity +matrix which can be used to zero out selected elements of a vector. A +Householder matrix @math{P} takes the form, +@tex +\beforedisplay +$$ +P = I - \tau v v^T +$$ +\afterdisplay +@end tex +@ifinfo + +@example +P = I - \tau v v^T +@end example + +@end ifinfo +@noindent +where @math{v} is a vector (called the @dfn{Householder vector}) and +@math{\tau = 2/(v^T v)}. The functions described in this section use the +rank-1 structure of the Householder matrix to create and apply +Householder transformations efficiently. + +@deftypefun double gsl_linalg_householder_transform (gsl_vector * @var{v}) +This function prepares a Householder transformation @math{P = I - \tau v +v^T} which can be used to zero all the elements of the input vector except +the first. On output the transformation is stored in the vector @var{v} +and the scalar @math{\tau} is returned. +@end deftypefun + +@deftypefun int gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A) +This function applies the Householder matrix @math{P} defined by the +scalar @var{tau} and the vector @var{v} to the left-hand side of the +matrix @var{A}. On output the result @math{P A} is stored in @var{A}. +@end deftypefun + +@deftypefun int gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A) +This function applies the Householder matrix @math{P} defined by the +scalar @var{tau} and the vector @var{v} to the right-hand side of the +matrix @var{A}. On output the result @math{A P} is stored in @var{A}. +@end deftypefun + +@deftypefun int gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w) +This function applies the Householder transformation @math{P} defined by +the scalar @var{tau} and the vector @var{v} to the vector @var{w}. On +output the result @math{P w} is stored in @var{w}. +@end deftypefun + +@comment @deftypefun int gsl_linalg_householder_hm1 (double tau, gsl_matrix * A) +@comment This function applies the Householder transform, defined by the scalar +@comment @var{tau} and the vector @var{v}, to a matrix being build up from the +@comment identity matrix, using the first column of @var{A} as a householder vector. +@comment @end deftypefun + +@node Householder solver for linear systems +@section Householder solver for linear systems +@cindex solution of linear system by Householder transformations +@cindex Householder linear solver + +@deftypefun int gsl_linalg_HH_solve (gsl_matrix * @var{A}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the system @math{A x = b} directly using +Householder transformations. On output the solution is stored in @var{x} +and @var{b} is not modified. The matrix @var{A} is destroyed by the +Householder transformations. +@end deftypefun + +@deftypefun int gsl_linalg_HH_svx (gsl_matrix * @var{A}, gsl_vector * @var{x}) +This function solves the system @math{A x = b} in-place using +Householder transformations. On input @var{x} should contain the +right-hand side @math{b}, which is replaced by the solution on output. The +matrix @var{A} is destroyed by the Householder transformations. +@end deftypefun + +@node Tridiagonal Systems +@section Tridiagonal Systems +@cindex tridiagonal systems + +@deftypefun int gsl_linalg_solve_tridiag (const gsl_vector * @var{diag}, const gsl_vector * @var{e}, const gsl_vector * @var{f}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the general @math{N}-by-@math{N} system @math{A x = +b} where @var{A} is tridiagonal (@c{$N\geq 2$} +@math{N >= 2}). The super-diagonal and +sub-diagonal vectors @var{e} and @var{f} must be one element shorter +than the diagonal vector @var{diag}. The form of @var{A} for the 4-by-4 +case is shown below, +@tex +\beforedisplay +$$ +A = \pmatrix{d_0&e_0& 0& 0\cr + f_0&d_1&e_1& 0\cr + 0 &f_1&d_2&e_2\cr + 0 &0 &f_2&d_3\cr} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = ( d_0 e_0 0 0 ) + ( f_0 d_1 e_1 0 ) + ( 0 f_1 d_2 e_2 ) + ( 0 0 f_2 d_3 ) +@end example +@end ifinfo +@end deftypefun + +@deftypefun int gsl_linalg_solve_symm_tridiag (const gsl_vector * @var{diag}, const gsl_vector * @var{e}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the general @math{N}-by-@math{N} system @math{A x = +b} where @var{A} is symmetric tridiagonal (@c{$N\geq 2$} +@math{N >= 2}). The off-diagonal vector +@var{e} must be one element shorter than the diagonal vector @var{diag}. +The form of @var{A} for the 4-by-4 case is shown below, +@tex +\beforedisplay +$$ +A = \pmatrix{d_0&e_0& 0& 0\cr + e_0&d_1&e_1& 0\cr + 0 &e_1&d_2&e_2\cr + 0 &0 &e_2&d_3\cr} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = ( d_0 e_0 0 0 ) + ( e_0 d_1 e_1 0 ) + ( 0 e_1 d_2 e_2 ) + ( 0 0 e_2 d_3 ) +@end example +@end ifinfo +The current implementation uses a variant of Cholesky decomposition +which can cause division by zero if the matrix is not positive definite. +@end deftypefun + +@deftypefun int gsl_linalg_solve_cyc_tridiag (const gsl_vector * @var{diag}, const gsl_vector * @var{e}, const gsl_vector * @var{f}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the general @math{N}-by-@math{N} system @math{A x = +b} where @var{A} is cyclic tridiagonal (@c{$N\geq 3$} +@math{N >= 3}). The cyclic super-diagonal and +sub-diagonal vectors @var{e} and @var{f} must have the same number of +elements as the diagonal vector @var{diag}. The form of @var{A} for the +4-by-4 case is shown below, +@tex +\beforedisplay +$$ +A = \pmatrix{d_0&e_0& 0 &f_3\cr + f_0&d_1&e_1& 0 \cr + 0 &f_1&d_2&e_2\cr + e_3& 0 &f_2&d_3\cr} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = ( d_0 e_0 0 f_3 ) + ( f_0 d_1 e_1 0 ) + ( 0 f_1 d_2 e_2 ) + ( e_3 0 f_2 d_3 ) +@end example +@end ifinfo +@end deftypefun + + +@deftypefun int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * @var{diag}, const gsl_vector * @var{e}, const gsl_vector * @var{b}, gsl_vector * @var{x}) +This function solves the general @math{N}-by-@math{N} system @math{A x = +b} where @var{A} is symmetric cyclic tridiagonal (@c{$N\geq 3$} +@math{N >= 3}). The cyclic +off-diagonal vector @var{e} must have the same number of elements as the +diagonal vector @var{diag}. The form of @var{A} for the 4-by-4 case is +shown below, +@tex +\beforedisplay +$$ +A = \pmatrix{d_0&e_0& 0 &e_3\cr + e_0&d_1&e_1& 0 \cr + 0 &e_1&d_2&e_2\cr + e_3& 0 &e_2&d_3\cr} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A = ( d_0 e_0 0 e_3 ) + ( e_0 d_1 e_1 0 ) + ( 0 e_1 d_2 e_2 ) + ( e_3 0 e_2 d_3 ) +@end example +@end ifinfo +@end deftypefun + +@node Balancing +@section Balancing +@cindex balancing matrices + +The process of balancing a matrix applies similarity transformations +to make the rows and columns have comparable norms. This is +useful, for example, to reduce roundoff errors in the solution +of eigenvalue problems. Balancing a matrix @math{A} consists +of replacing @math{A} with a similar matrix +@tex +\beforedisplay +$$ +A' = D^{-1} A D +$$ +\afterdisplay +@end tex +@ifinfo + +@example +A' = D^(-1) A D +@end example + +@end ifinfo +where @math{D} is a diagonal matrix whose entries are powers +of the floating point radix. + +@deftypefun int gsl_linalg_balance_matrix (gsl_matrix * @var{A}, gsl_vector * @var{D}) +This function replaces the matrix @var{A} with its balanced counterpart +and stores the diagonal elements of the similarity transformation +into the vector @var{D}. +@end deftypefun + +@node Linear Algebra Examples +@section Examples + +The following program solves the linear system @math{A x = b}. The +system to be solved is, +@tex +\beforedisplay +$$ +\left( +\matrix{0.18& 0.60& 0.57& 0.96\cr +0.41& 0.24& 0.99& 0.58\cr +0.14& 0.30& 0.97& 0.66\cr +0.51& 0.13& 0.19& 0.85} +\right) +\left( +\matrix{x_0\cr +x_1\cr +x_2\cr +x_3} +\right) += +\left( +\matrix{1.0\cr +2.0\cr +3.0\cr +4.0} +\right) +$$ +\afterdisplay +@end tex +@ifinfo + +@example +[ 0.18 0.60 0.57 0.96 ] [x0] [1.0] +[ 0.41 0.24 0.99 0.58 ] [x1] = [2.0] +[ 0.14 0.30 0.97 0.66 ] [x2] [3.0] +[ 0.51 0.13 0.19 0.85 ] [x3] [4.0] +@end example + +@end ifinfo +@noindent +and the solution is found using LU decomposition of the matrix @math{A}. + +@example +@verbatiminclude examples/linalglu.c +@end example + +@noindent +Here is the output from the program, + +@example +@verbatiminclude examples/linalglu.out +@end example + +@noindent +This can be verified by multiplying the solution @math{x} by the +original matrix @math{A} using @sc{gnu octave}, + +@example +octave> A = [ 0.18, 0.60, 0.57, 0.96; + 0.41, 0.24, 0.99, 0.58; + 0.14, 0.30, 0.97, 0.66; + 0.51, 0.13, 0.19, 0.85 ]; + +octave> x = [ -4.05205; -12.6056; 1.66091; 8.69377]; + +octave> A * x +ans = + 1.0000 + 2.0000 + 3.0000 + 4.0000 +@end example + +@noindent +This reproduces the original right-hand side vector, @math{b}, in +accordance with the equation @math{A x = b}. + +@node Linear Algebra References and Further Reading +@section References and Further Reading + +Further information on the algorithms described in this section can be +found in the following book, + +@itemize @asis +@item +G. H. Golub, C. F. Van Loan, @cite{Matrix Computations} (3rd Ed, 1996), +Johns Hopkins University Press, ISBN 0-8018-5414-8. +@end itemize + +@noindent +The @sc{lapack} library is described in the following manual, + +@itemize @asis +@item +@cite{LAPACK Users' Guide} (Third Edition, 1999), Published by SIAM, +ISBN 0-89871-447-8. + +@uref{http://www.netlib.org/lapack} +@end itemize + +@noindent +The @sc{lapack} source code can be found at the website above, along +with an online copy of the users guide. + +@noindent +The Modified Golub-Reinsch algorithm is described in the following paper, + +@itemize @asis +@item +T.F. Chan, ``An Improved Algorithm for Computing the Singular Value +Decomposition'', @cite{ACM Transactions on Mathematical Software}, 8 +(1982), pp 72--83. +@end itemize + +@noindent +The Jacobi algorithm for singular value decomposition is described in +the following papers, + +@itemize @asis +@item +J.C. Nash, ``A one-sided transformation method for the singular value +decomposition and algebraic eigenproblem'', @cite{Computer Journal}, +Volume 18, Number 1 (1973), p 74--76 + +@item +James Demmel, Kresimir Veselic, ``Jacobi's Method is more accurate than +QR'', @cite{Lapack Working Note 15} (LAWN-15), October 1989. Available +from netlib, @uref{http://www.netlib.org/lapack/} in the @code{lawns} or +@code{lawnspdf} directories. +@end itemize + + + |