Dogcows Code - chaz/yoink/blob - src/moof/cml/matrix/lu.h

   1 /* -*- C++ -*- ------------------------------------------------------------
   2
   3 Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
   4
   5 The Configurable Math Library (CML) is distributed under the terms of the
   6 Boost Software License, v1.0 (see cml/LICENSE for details).
   7
   8  *-----------------------------------------------------------------------*/
   9 /** @file
  10  *  @brief Implements LU decomposition for square matrix expressions.
  11  *
  12  * @todo The LU implementation does not check for a zero diagonal entry
  13  * (implying that the input has no LU factorization).
  14  *
  15  * @todo Should also have a pivoting implementation.
  16  *
  17  * @todo need to throw a numeric error if the determinant of the matrix
  18  * given to lu(), lu_solve(), or inverse() is 0.
  19  *
  20  * @internal The implementation is the same for fixed- and dynamic-size
  21  * matrices.  It can be sped up for small matrices later.
  22  */
  23
  24 #ifndef lu_h
  25 #define lu_h
  26
  27 #include <cml/et/size_checking.h>
  28 #include <cml/matrix/matrix_expr.h>
  29 #include <cml/matvec/matvec_promotions.h>
  30
  31 /* This is used below to create a more meaningful compile-time error when
  32  * lu is not provided with a matrix or MatrixExpr argument:
  33  */
  34 struct lu_expects_a_matrix_arg_error;
  35
  36 /* This is used below to create a more meaningful compile-time error when
  37  * lu_inplace is not provided with an assignable matrix argument:
  38  */
  39 struct lu_inplace_expects_an_assignable_matrix_arg_error;
  40
  41 namespace cml {
  42 namespace detail {
  43
  44 /* Compute the LU decomposition in-place: */
  45 template<class MatT> inline
  46 void lu_inplace(MatT& A)
  47 {
  48     /* Shorthand: */
  49     typedef et::ExprTraits<MatT> arg_traits;
  50     typedef typename arg_traits::result_tag arg_result;
  51     typedef typename arg_traits::assignable_tag arg_assignment;
  52     typedef typename arg_traits::size_tag size_tag;
  53     typedef typename arg_traits::value_type value_type;
  54
  55     /* lu_inplace() requires an assignable matrix expression: */
  56     CML_STATIC_REQUIRE_M(
  57         (same_type<arg_result, et::matrix_result_tag>::is_true
  58          && same_type<arg_assignment, et::assignable_tag>::is_true),
  59         lu_inplace_expects_an_assignable_matrix_arg_error);
  60     /* Note: parens are required here so that the preprocessor ignores the
  61      * commas.
  62      */
  63
  64     /* Verify that the matrix is square, and get the size: */
  65     ssize_t N = (ssize_t) cml::et::CheckedSquare(A, size_tag());
  66
  67
  68     for(ssize_t k = 0; k < N-1; ++k) {
  69         /* XXX Should check if A(k,k) = 0! */
  70         for(ssize_t i = k+1; i < N; ++i) {
  71             value_type n = (A(i,k) /= A(k,k));
  72             for(ssize_t j = k+1; j < N; ++ j) {
  73                 A(i,j) -= n*A(k,j);
  74             }
  75         }
  76     }
  77 }
  78
  79 /* Compute the LU decomposition, and return a copy of the result: */
  80 template<class MatT>
  81 inline typename MatT::temporary_type
  82 lu_copy(const MatT& M)
  83 {
  84     /* Shorthand: */
  85     typedef et::ExprTraits<MatT> arg_traits;
  86     typedef typename arg_traits::result_tag arg_result;
  87
  88     /* lu_with_copy() requires a matrix expression: */
  89     CML_STATIC_REQUIRE_M(
  90         (same_type<arg_result, et::matrix_result_tag>::is_true),
  91         lu_expects_a_matrix_arg_error);
  92     /* Note: parens are required here so that the preprocessor ignores the
  93      * commas.
  94      */
  95
  96     /* Use the in-place LU function, and return the result: */
  97     typename MatT::temporary_type A;
  98     cml::et::detail::Resize(A,M.rows(),M.cols());
  99     A = M;
 100     lu_inplace(A);
 101     return A;
 102 }
 103
 104 } // namespace detail
 105
 106 /** LU factorization for a matrix. */
 107 template<typename E, class AT, typename BO, class L>
 108 inline typename matrix<E,AT,BO,L>::temporary_type
 109 lu(const matrix<E,AT,BO,L>& m)
 110 {
 111     return detail::lu_copy(m);
 112 }
 113
 114 /** LU factorization for a matrix expression. */
 115 template<typename XprT>
 116 inline typename et::MatrixXpr<XprT>::temporary_type
 117 lu(const et::MatrixXpr<XprT>& e)
 118 {
 119     return detail::lu_copy(e);
 120 }
 121
 122 /** Solve y = LUx for x.
 123  *
 124  * This solves Lb = y for b by forward substitution, then Ux = b for x by
 125  * backward substitution.
 126  */
 127 template<typename MatT, typename VecT> inline
 128 typename et::MatVecPromote<MatT,VecT>::temporary_type
 129 lu_solve(const MatT& LU, const VecT& b)
 130 {
 131     /* Shorthand. */
 132     typedef et::ExprTraits<MatT> lu_traits;
 133     typedef typename et::MatVecPromote<MatT,VecT>::temporary_type vector_type;
 134     typedef typename vector_type::value_type value_type;
 135
 136     /* Verify that the matrix is square, and get the size: */
 137     ssize_t N = (ssize_t) cml::et::CheckedSquare(
 138             LU, typename lu_traits::size_tag());
 139
 140     /* Verify that the matrix and vector have compatible sizes: */
 141     et::CheckedSize(LU, b, typename vector_type::size_tag());
 142
 143     /* Solve Ly = b for y by forward substitution.  The entries below the
 144      * diagonal of LU correspond to L, understood to be below a diagonal of
 145      * 1's:
 146      */
 147     vector_type y; cml::et::detail::Resize(y,N);
 148     for(ssize_t i = 0; i < N; ++i) {
 149         y[i] = b[i];
 150         for(ssize_t j = 0; j < i; ++j) {
 151             y[i] -= LU(i,j)*y[j];
 152         }
 153     }
 154
 155     /* Solve Ux = y for x by backward substitution.  The entries at and above
 156      * the diagonal of LU correspond to U:
 157      */
 158     vector_type x; cml::et::detail::Resize(x,N);
 159     for(ssize_t i = N-1; i >= 0; --i) {
 160         x[i] = y[i];
 161         for(ssize_t j = i+1; j < N; ++j) {
 162             x[i] -= LU(i,j)*x[j];
 163         }
 164         x[i] /= LU(i,i);
 165     }
 166
 167     /* Return x: */
 168     return x;
 169 }
 170
 171 } // namespace cml
 172
 173 #endif
 174
 175 // -------------------------------------------------------------------------
 176 // vim:ft=cpp