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
00027 (
00028 typename T1::elem_type& out_val,
00029 typename T1::pod_type& out_sign,
00030 const Base<typename T1::elem_type,T1>& X,
00031 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
00032 )
00033 {
00034 arma_extra_debug_sigprint();
00035
00036 typedef typename T1::elem_type eT;
00037
00038 const unwrap<T1> tmp(X.get_ref());
00039 const Mat<eT>& A = tmp.M;
00040
00041 arma_debug_check( !A.is_square(), "log_det(): matrix must be square" );
00042
00043 auxlib::log_det(out_val, out_sign, A);
00044 }
00045
00046
00047
00048 template<typename T1>
00049 inline
00050 void
00051 log_det
00052 (
00053 typename T1::elem_type& out_val,
00054 typename T1::pod_type& out_sign,
00055 const Op<T1,op_diagmat>& X,
00056 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
00057 )
00058 {
00059 arma_extra_debug_sigprint();
00060
00061 typedef typename T1::elem_type eT;
00062 typedef typename T1::pod_type T;
00063
00064 const diagmat_proxy<T1> A(X.m);
00065
00066 const u32 N = A.n_elem;
00067
00068 arma_debug_check( (N == 0), "log_det(): given matrix has no elements" );
00069
00070 const eT x = A[0];
00071
00072 T sign = (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
00073 eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
00074
00075 for(u32 i=1; i<N; ++i)
00076 {
00077 const eT x = A[i];
00078
00079 sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
00080 val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
00081 }
00082
00083 out_val = val;
00084 out_sign = sign;
00085 }
00086
00087
00088
00089