fn_log_det.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // 
00007 // This file is part of the Armadillo C++ library.
00008 // It is provided without any warranty of fitness
00009 // for any purpose. You can redistribute this file
00010 // and/or modify it under the terms of the GNU
00011 // Lesser General Public License (LGPL) as published
00012 // by the Free Software Foundation, either version 3
00013 // of the License or (at your option) any later version.
00014 // (see http://www.opensource.org/licenses for more info)
00015 
00016 
00017 //! \addtogroup fn_log_det
00018 //! @{
00019 
00020 
00021 
00022 //! log determinant of mat
00023 template<typename T1>
00024 inline
00025 void
00026 log_det(typename T1::elem_type& out_val, typename T1::pod_type& out_sign, const Base<typename T1::elem_type,T1>& X)
00027   {
00028   arma_extra_debug_sigprint();
00029   
00030   typedef typename T1::elem_type eT;
00031   
00032   const unwrap<T1>   tmp(X.get_ref());
00033   const Mat<eT>& A = tmp.M;
00034   
00035   arma_debug_check( !A.is_square(), "log_det(): matrix must be square" );
00036   
00037   auxlib::log_det(out_val, out_sign, A);
00038   }
00039 
00040 
00041 
00042 template<typename T1>
00043 inline
00044 void
00045 log_det(typename T1::elem_type& out_val, typename T1::pod_type& out_sign, const Op<T1,op_diagmat>& X)
00046   {
00047   arma_extra_debug_sigprint();
00048   
00049   typedef typename T1::elem_type eT;
00050   typedef typename T1::pod_type   T;
00051   
00052   const diagmat_proxy<T1> A(X.m);
00053   
00054   const u32 N = A.n_elem;
00055   
00056   arma_debug_check( (N == 0), "log_det(): given matrix has no elements" );
00057   
00058   const eT x = A[0];
00059   
00060   T  sign = (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
00061   eT val  = (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
00062 
00063   for(u32 i=1; i<N; ++i)
00064     {
00065     const eT x = A[i];
00066     
00067     sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
00068     val  += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
00069     }
00070   
00071   out_val  = val;
00072   out_sign = sign;
00073   }
00074 
00075 
00076 
00077 //! @}