X-Git-Url: https://git.dogcows.com/gitweb?p=chaz%2Fyoink;a=blobdiff_plain;f=src%2Fcml%2Fmatrix%2Fdeterminant.h;fp=src%2Fcml%2Fmatrix%2Fdeterminant.h;h=fa2cba300443841cb39fc32190b65a599b61f315;hp=0000000000000000000000000000000000000000;hb=0fffd0097d7b496454413e57b398c903ecc252e4;hpb=79becf045222f385da5a1b9eb79081f6f5266c86 diff --git a/src/cml/matrix/determinant.h b/src/cml/matrix/determinant.h new file mode 100644 index 0000000..fa2cba3 --- /dev/null +++ b/src/cml/matrix/determinant.h @@ -0,0 +1,192 @@ +/* -*- 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 Compute the determinant of a square matrix using LU factorization. + * + * @todo This should be specialized on the matrix size for small matrices. + */ + +#ifndef determinant_h +#define determinant_h + +#include + +namespace cml { +namespace detail { + +/* Need to use a functional, since template functions cannot be + * specialized. N is used to differentiate dimension only, so this can be + * used for any matrix size type: + */ +template struct determinant_f; + +/* 2x2 determinant. Despite being marked for fixed_size matrices, this can + * be used for dynamic-sized ones also: + */ +template +struct determinant_f +{ + typename MatT::value_type operator()(const MatT& M) const + { + return M(0,0)*M(1,1) - M(1,0)*M(0,1); + } + +}; + +/* 3x3 determinant. Despite being marked for fixed_size matrices, this can + * be used for dynamic-sized ones also: + */ +template +struct determinant_f +{ + /* [00 01 02] + * M = [10 11 12] + * [20 21 22] + */ + typename MatT::value_type operator()(const MatT& M) const + { + return M(0,0)*(M(1,1)*M(2,2) - M(1,2)*M(2,1)) + + M(0,1)*(M(1,2)*M(2,0) - M(1,0)*M(2,2)) + + M(0,2)*(M(1,0)*M(2,1) - M(1,1)*M(2,0)); + } + +}; + +/* 4x4 determinant. Despite being marked for fixed_size matrices, this can + * be used for dynamic-sized ones also: + */ +template +struct determinant_f +{ + /* [00 01 02 03] + * M = [10 11 12 13] + * [20 21 22 23] + * [30 31 32 33] + * + * |11 12 13| |10 12 13| + * C00 = |21 22 23| C01 = |20 22 23| + * |31 32 33| |30 32 33| + * + * |10 11 13| |10 11 12| + * C02 = |20 21 23| C03 = |20 21 22| + * |30 31 33| |30 31 32| + * + * d00 = 11 * (22*33 - 23*32) d01 = 10 * (22*33 - 23*32) + * + 12 * (23*31 - 21*33) + 12 * (23*30 - 20*33) + * + 13 * (21*32 - 22*31) + 13 * (20*32 - 22*30) + * + * d02 = 10 * (21*33 - 23*31) d03 = 10 * (21*32 - 22*31) + * + 11 * (23*30 - 20*33) + 11 * (22*30 - 20*32) + * + 13 * (20*31 - 21*30) + 12 * (20*31 - 21*30) + */ + typename MatT::value_type operator()(const MatT& M) const + { + /* Shorthand. */ + typedef typename MatT::value_type value_type; + + /* Common cofactors: */ + value_type m_22_33_23_32 = M(2,2)*M(3,3) - M(2,3)*M(3,2); + value_type m_23_30_20_33 = M(2,3)*M(3,0) - M(2,0)*M(3,3); + value_type m_20_31_21_30 = M(2,0)*M(3,1) - M(2,1)*M(3,0); + value_type m_21_32_22_31 = M(2,1)*M(3,2) - M(2,2)*M(3,1); + value_type m_23_31_21_33 = M(2,3)*M(3,1) - M(2,1)*M(3,3); + value_type m_20_32_22_30 = M(2,0)*M(3,2) - M(2,2)*M(3,0); + + value_type d00 = M(0,0)*( + M(1,1) * m_22_33_23_32 + + M(1,2) * m_23_31_21_33 + + M(1,3) * m_21_32_22_31); + + value_type d01 = M(0,1)*( + M(1,0) * m_22_33_23_32 + + M(1,2) * m_23_30_20_33 + + M(1,3) * m_20_32_22_30); + + value_type d02 = M(0,2)*( + M(1,0) * - m_23_31_21_33 + + M(1,1) * m_23_30_20_33 + + M(1,3) * m_20_31_21_30); + + value_type d03 = M(0,3)*( + M(1,0) * m_21_32_22_31 + + M(1,1) * - m_20_32_22_30 + + M(1,2) * m_20_31_21_30); + + return d00 - d01 + d02 - d03; + } + +}; + +/* General NxN determinant by LU factorization: */ +template +struct determinant_f +{ + typename MatT::value_type operator()(const MatT& M) const + { + /* Compute the LU factorization: */ + typename MatT::temporary_type LU = lu(M); + + /* The product of the diagonal entries is the determinant: */ + typename MatT::value_type det = LU(0,0); + for(size_t i = 1; i < LU.rows(); ++ i) + det *= LU(i,i); + return det; + } + +}; + +/* Generator for the determinant functional for fixed-size matrices: */ +template typename MatT::value_type +determinant(const MatT& M, fixed_size_tag) +{ + /* Require a square matrix: */ + cml::et::CheckedSquare(M, fixed_size_tag()); + return determinant_f()(M); +} + +/* Generator for the determinant functional for dynamic-size matrices: */ +template typename MatT::value_type +determinant(const MatT& M, dynamic_size_tag) +{ + /* Require a square matrix: */ + cml::et::CheckedSquare(M, dynamic_size_tag()); + + /* Dispatch based upon the matrix dimension: */ + switch(M.rows()) { + case 2: return determinant_f()(M); + case 3: return determinant_f()(M); + case 4: return determinant_f()(M); + default: return determinant_f()(M); // > 4x4. + } +} + +} // namespace detail + +/** Determinant of a matrix. */ +template inline E +determinant(const matrix& M) +{ + typedef typename matrix::size_tag size_tag; + return detail::determinant(M,size_tag()); +} + +/** Determinant of a matrix expression. */ +template inline typename XprT::value_type +determinant(const et::MatrixXpr& e) +{ + typedef typename et::MatrixXpr::size_tag size_tag; + return detail::determinant(e,size_tag()); +} + +} // namespace cml + +#endif + +// ------------------------------------------------------------------------- +// vim:ft=cpp