Chris@49
|
1 // Copyright (C) 2010-2011 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2010-2011 Conrad Sanderson
|
Chris@49
|
3 //
|
Chris@49
|
4 // This Source Code Form is subject to the terms of the Mozilla Public
|
Chris@49
|
5 // License, v. 2.0. If a copy of the MPL was not distributed with this
|
Chris@49
|
6 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
Chris@49
|
7
|
Chris@49
|
8
|
Chris@49
|
9 //! \addtogroup fn_log_det
|
Chris@49
|
10 //! @{
|
Chris@49
|
11
|
Chris@49
|
12
|
Chris@49
|
13
|
Chris@49
|
14 //! log determinant of mat
|
Chris@49
|
15 template<typename T1>
|
Chris@49
|
16 inline
|
Chris@49
|
17 bool
|
Chris@49
|
18 log_det
|
Chris@49
|
19 (
|
Chris@49
|
20 typename T1::elem_type& out_val,
|
Chris@49
|
21 typename T1::pod_type& out_sign,
|
Chris@49
|
22 const Base<typename T1::elem_type,T1>& X,
|
Chris@49
|
23 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
|
Chris@49
|
24 )
|
Chris@49
|
25 {
|
Chris@49
|
26 arma_extra_debug_sigprint();
|
Chris@49
|
27 arma_ignore(junk);
|
Chris@49
|
28
|
Chris@49
|
29 return auxlib::log_det(out_val, out_sign, X);
|
Chris@49
|
30 }
|
Chris@49
|
31
|
Chris@49
|
32
|
Chris@49
|
33
|
Chris@49
|
34 template<typename T1>
|
Chris@49
|
35 inline
|
Chris@49
|
36 void
|
Chris@49
|
37 log_det
|
Chris@49
|
38 (
|
Chris@49
|
39 typename T1::elem_type& out_val,
|
Chris@49
|
40 typename T1::pod_type& out_sign,
|
Chris@49
|
41 const Op<T1,op_diagmat>& X,
|
Chris@49
|
42 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
|
Chris@49
|
43 )
|
Chris@49
|
44 {
|
Chris@49
|
45 arma_extra_debug_sigprint();
|
Chris@49
|
46 arma_ignore(junk);
|
Chris@49
|
47
|
Chris@49
|
48 typedef typename T1::elem_type eT;
|
Chris@49
|
49 typedef typename T1::pod_type T;
|
Chris@49
|
50
|
Chris@49
|
51 const diagmat_proxy<T1> A(X.m);
|
Chris@49
|
52
|
Chris@49
|
53 const uword N = A.n_elem;
|
Chris@49
|
54
|
Chris@49
|
55 if(N == 0)
|
Chris@49
|
56 {
|
Chris@49
|
57 out_val = eT(0);
|
Chris@49
|
58 out_sign = T(1);
|
Chris@49
|
59
|
Chris@49
|
60 return;
|
Chris@49
|
61 }
|
Chris@49
|
62
|
Chris@49
|
63 eT x = A[0];
|
Chris@49
|
64
|
Chris@49
|
65 T sign = (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
|
Chris@49
|
66 eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
|
Chris@49
|
67
|
Chris@49
|
68 for(uword i=1; i<N; ++i)
|
Chris@49
|
69 {
|
Chris@49
|
70 x = A[i];
|
Chris@49
|
71
|
Chris@49
|
72 sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
|
Chris@49
|
73 val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
|
Chris@49
|
74 }
|
Chris@49
|
75
|
Chris@49
|
76 out_val = val;
|
Chris@49
|
77 out_sign = sign;
|
Chris@49
|
78 }
|
Chris@49
|
79
|
Chris@49
|
80
|
Chris@49
|
81
|
Chris@49
|
82 //! @}
|