+++ /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 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 <cml/matrix/lu.h>
-
-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<typename MatT, int N> struct determinant_f;
-
-/* 2x2 determinant. Despite being marked for fixed_size matrices, this can
- * be used for dynamic-sized ones also:
- */
-template<typename MatT>
-struct determinant_f<MatT,2>
-{
- 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<typename MatT>
-struct determinant_f<MatT,3>
-{
- /* [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<typename MatT>
-struct determinant_f<MatT,4>
-{
- /* [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<typename MatT, int N>
-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> 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<MatT,MatT::array_rows>()(M);
-}
-
-/* Generator for the determinant functional for dynamic-size matrices: */
-template<typename MatT> 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<MatT,2>()(M);
- case 3: return determinant_f<MatT,3>()(M);
- case 4: return determinant_f<MatT,4>()(M);
- default: return determinant_f<MatT,0>()(M); // > 4x4.
- }
-}
-
-} // namespace detail
-
-/** Determinant of a matrix. */
-template<typename E, class AT, class BO, class L> inline E
-determinant(const matrix<E,AT,BO,L>& M)
-{
- typedef typename matrix<E,AT,BO,L>::size_tag size_tag;
- return detail::determinant(M,size_tag());
-}
-
-/** Determinant of a matrix expression. */
-template<typename XprT> inline typename XprT::value_type
-determinant(const et::MatrixXpr<XprT>& e)
-{
- typedef typename et::MatrixXpr<XprT>::size_tag size_tag;
- return detail::determinant(e,size_tag());
-}
-
-} // namespace cml
-
-#endif
-
-// -------------------------------------------------------------------------
-// vim:ft=cpp
projects / chaz / yoink / blobdiff