| @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 |
| |
| |
| |