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
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 //! @}