--- /dev/null
+/* -*- 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
projects / chaz / yoink / blobdiff