linalg.texi

来自「开放gsl矩阵运算」· TEXI 代码 · 共 906 行 · 第 1/3 页

tridiagonal decomposition (@var{A}, @var{tau}) obtained from
@code{gsl_linalg_hermtd_decomp} into the real vectors @var{d} and
@var{sd}.
@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}.
@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 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_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.  The form of @var{A} for the
4-by-4 case is shown below,

@tex
\beforedisplay
$$
A = \pmatrix{d_0&e_0&   &  \cr
             e_0&d_1&e_1&  \cr
                &e_1&d_2&e_2\cr 
                &   &e_2&d_3\cr}
$$
\afterdisplay
@end tex
@ifinfo
@example
A = ( d_0 e_0         )
    ( e_0 d_1 e_1     )
    (     e_1 d_2 e_2 )
    (         e_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.  The form of @var{A}
for the 4-by-4 case is shown below,

@tex
\beforedisplay
$$
A = \pmatrix{d_0&e_0&   &e_3\cr
             e_0&d_1&e_1&   \cr
                &e_1&d_2&e_2\cr 
             e_3&   &e_2&d_3\cr}
$$
\afterdisplay
@end tex
@ifinfo
@example
A = ( d_0 e_0     e_3 )
    ( e_0 d_1 e_1     )
    (     e_1 d_2 e_2 )
    ( e_3     e_2 d_3 )
@end example
@end ifinfo

@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
#include <stdio.h>
#include <gsl/gsl_linalg.h>

int
main (void)
@{
  double a_data[] = @{ 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 @};

  double b_data[] = @{ 1.0, 2.0, 3.0, 4.0 @};

  gsl_matrix_view m 
    = gsl_matrix_view_array(a_data, 4, 4);

  gsl_vector_view b
    = gsl_vector_view_array(b_data, 4);

  gsl_vector *x = gsl_vector_alloc (4);
  
  int s;

  gsl_permutation * p = gsl_permutation_alloc (4);

  gsl_linalg_LU_decomp (&m.matrix, p, &s);

  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);

  printf ("x = \n");
  gsl_vector_fprintf(stdout, x, "%g");

  gsl_permutation_free (p);
  return 0;
@}
@end example
@noindent
Here is the output from the program,

@example
x = -4.05205
-12.6056
1.66091
8.69377
@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
@noindent
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,

@itemize @asis
@item
@cite{LAPACK Users' Guide} (Third Edition, 1999), Published by SIAM,
ISBN 0-89871-447-8.

@url{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, @url{http://www.netlib.org/lapack/} in the @code{lawns} or
@code{lawnspdf} directories.
@end itemize