Initial commit

Lib/Include/CML/matrix/lu.h

@@ -0,0 +1,176 @@
+/* -*- C++ -*- ------------------------------------------------------------
+ 
+Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
+
+The Configurable Math Library (CML) is distributed under the terms of the
+Boost Software License, v1.0 (see cml/LICENSE for details).
+
+ *-----------------------------------------------------------------------*/
+/** @file
+ *  @brief Implements LU decomposition for square matrix expressions.
+ *
+ * @todo The LU implementation does not check for a zero diagonal entry
+ * (implying that the input has no LU factorization).
+ *
+ * @todo Should also have a pivoting implementation.
+ *
+ * @todo need to throw a numeric error if the determinant of the matrix
+ * given to lu(), lu_solve(), or inverse() is 0.
+ *
+ * @internal The implementation is the same for fixed- and dynamic-size
+ * matrices.  It can be sped up for small matrices later.
+ */
+
+#ifndef lu_h
+#define lu_h
+
+#include <cml/et/size_checking.h>
+#include <cml/matrix/matrix_expr.h>
+#include <cml/matvec/matvec_promotions.h>
+
+/* This is used below to create a more meaningful compile-time error when
+ * lu is not provided with a matrix or MatrixExpr argument:
+ */
+struct lu_expects_a_matrix_arg_error;
+
+/* This is used below to create a more meaningful compile-time error when
+ * lu_inplace is not provided with an assignable matrix argument:
+ */
+struct lu_inplace_expects_an_assignable_matrix_arg_error;
+
+namespace cml {
+namespace detail {
+
+/* Compute the LU decomposition in-place: */
+template<class MatT> inline
+void lu_inplace(MatT& A)
+{
+    /* Shorthand: */
+    typedef et::ExprTraits<MatT> arg_traits;
+    typedef typename arg_traits::result_tag arg_result;
+    typedef typename arg_traits::assignable_tag arg_assignment;
+    typedef typename arg_traits::size_tag size_tag;
+    typedef typename arg_traits::value_type value_type;
+
+    /* lu_inplace() requires an assignable matrix expression: */
+    CML_STATIC_REQUIRE_M(
+        (same_type<arg_result, et::matrix_result_tag>::is_true
+         && same_type<arg_assignment, et::assignable_tag>::is_true),
+        lu_inplace_expects_an_assignable_matrix_arg_error);
+    /* Note: parens are required here so that the preprocessor ignores the
+     * commas.
+     */
+
+    /* Verify that the matrix is square, and get the size: */
+    ssize_t N = (ssize_t) cml::et::CheckedSquare(A, size_tag());
+
+
+    for(ssize_t k = 0; k < N-1; ++k) {
+        /* XXX Should check if A(k,k) = 0! */
+        for(ssize_t i = k+1; i < N; ++i) {
+            value_type n = (A(i,k) /= A(k,k));
+            for(ssize_t j = k+1; j < N; ++ j) {
+                A(i,j) -= n*A(k,j);
+            }
+        }
+    }
+}
+
+/* Compute the LU decomposition, and return a copy of the result: */
+template<class MatT>
+inline typename MatT::temporary_type
+lu_copy(const MatT& M)
+{
+    /* Shorthand: */
+    typedef et::ExprTraits<MatT> arg_traits;
+    typedef typename arg_traits::result_tag arg_result;
+
+    /* lu_with_copy() requires a matrix expression: */
+    CML_STATIC_REQUIRE_M(
+        (same_type<arg_result, et::matrix_result_tag>::is_true),
+        lu_expects_a_matrix_arg_error);
+    /* Note: parens are required here so that the preprocessor ignores the
+     * commas.
+     */
+
+    /* Use the in-place LU function, and return the result: */
+    typename MatT::temporary_type A;
+    cml::et::detail::Resize(A,M.rows(),M.cols());
+    A = M;
+    lu_inplace(A);
+    return A;
+}
+
+} // namespace detail
+
+/** LU factorization for a matrix. */
+template<typename E, class AT, typename BO, class L>
+inline typename matrix<E,AT,BO,L>::temporary_type
+lu(const matrix<E,AT,BO,L>& m)
+{
+    return detail::lu_copy(m);
+}
+
+/** LU factorization for a matrix expression. */
+template<typename XprT>
+inline typename et::MatrixXpr<XprT>::temporary_type
+lu(const et::MatrixXpr<XprT>& e)
+{
+    return detail::lu_copy(e);
+}
+
+/** Solve y = LUx for x.
+ *
+ * This solves Lb = y for b by forward substitution, then Ux = b for x by
+ * backward substitution.
+ */
+template<typename MatT, typename VecT> inline
+typename et::MatVecPromote<MatT,VecT>::temporary_type
+lu_solve(const MatT& LU, const VecT& b)
+{
+    /* Shorthand. */
+    typedef et::ExprTraits<MatT> lu_traits;
+    typedef typename et::MatVecPromote<MatT,VecT>::temporary_type vector_type;
+    typedef typename vector_type::value_type value_type;
+
+    /* Verify that the matrix is square, and get the size: */
+    ssize_t N = (ssize_t) cml::et::CheckedSquare(
+            LU, typename lu_traits::size_tag());
+
+    /* Verify that the matrix and vector have compatible sizes: */
+    et::CheckedSize(LU, b, typename vector_type::size_tag());
+
+    /* Solve Ly = b for y by forward substitution.  The entries below the
+     * diagonal of LU correspond to L, understood to be below a diagonal of
+     * 1's:
+     */
+    vector_type y; cml::et::detail::Resize(y,N);
+    for(ssize_t i = 0; i < N; ++i) {
+        y[i] = b[i];
+        for(ssize_t j = 0; j < i; ++j) {
+            y[i] -= LU(i,j)*y[j];
+        }
+    }
+
+    /* Solve Ux = y for x by backward substitution.  The entries at and above
+     * the diagonal of LU correspond to U:
+     */
+    vector_type x; cml::et::detail::Resize(x,N);
+    for(ssize_t i = N-1; i >= 0; --i) {
+        x[i] = y[i];
+        for(ssize_t j = i+1; j < N; ++j) {
+            x[i] -= LU(i,j)*x[j];
+        }
+        x[i] /= LU(i,i);
+    }
+
+    /* Return x: */
+    return x;
+}
+
+} // namespace cml
+
+#endif
+
+// -------------------------------------------------------------------------
+// vim:ft=cpp