Chris@49: // Copyright (C) 2010-2011 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2010-2011 Conrad Sanderson
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup fn_log_det
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! log determinant of mat
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: bool
Chris@49: log_det
Chris@49:   (
Chris@49:         typename T1::elem_type&          out_val,
Chris@49:         typename T1::pod_type&           out_sign,
Chris@49:   const Base<typename T1::elem_type,T1>& X,
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   return auxlib::log_det(out_val, out_sign, X);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: void
Chris@49: log_det
Chris@49:   (
Chris@49:         typename T1::elem_type& out_val,
Chris@49:         typename T1::pod_type&  out_sign,
Chris@49:   const Op<T1,op_diagmat>&      X,
Chris@49:   const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   typedef typename T1::pod_type   T;
Chris@49:   
Chris@49:   const diagmat_proxy<T1> A(X.m);
Chris@49:   
Chris@49:   const uword N = A.n_elem;
Chris@49:   
Chris@49:   if(N == 0)
Chris@49:     {
Chris@49:     out_val  = eT(0);
Chris@49:     out_sign =  T(1);
Chris@49:     
Chris@49:     return;
Chris@49:     }
Chris@49:   
Chris@49:   eT x = A[0];
Chris@49:   
Chris@49:   T  sign = (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
Chris@49:   eT val  = (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
Chris@49:   
Chris@49:   for(uword i=1; i<N; ++i)
Chris@49:     {
Chris@49:     x = A[i];
Chris@49:     
Chris@49:     sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
Chris@49:     val  += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
Chris@49:     }
Chris@49:   
Chris@49:   out_val  = val;
Chris@49:   out_sign = sign;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}