fn_log_det.hpp
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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