|
max@0
|
1 // Copyright (C) 2010-2011 NICTA (www.nicta.com.au)
|
|
max@0
|
2 // Copyright (C) 2010-2011 Conrad Sanderson
|
|
max@0
|
3 //
|
|
max@0
|
4 // This file is part of the Armadillo C++ library.
|
|
max@0
|
5 // It is provided without any warranty of fitness
|
|
max@0
|
6 // for any purpose. You can redistribute this file
|
|
max@0
|
7 // and/or modify it under the terms of the GNU
|
|
max@0
|
8 // Lesser General Public License (LGPL) as published
|
|
max@0
|
9 // by the Free Software Foundation, either version 3
|
|
max@0
|
10 // of the License or (at your option) any later version.
|
|
max@0
|
11 // (see http://www.opensource.org/licenses for more info)
|
|
max@0
|
12
|
|
max@0
|
13
|
|
max@0
|
14 //! \addtogroup fn_log_det
|
|
max@0
|
15 //! @{
|
|
max@0
|
16
|
|
max@0
|
17
|
|
max@0
|
18
|
|
max@0
|
19 //! log determinant of mat
|
|
max@0
|
20 template<typename T1>
|
|
max@0
|
21 inline
|
|
max@0
|
22 bool
|
|
max@0
|
23 log_det
|
|
max@0
|
24 (
|
|
max@0
|
25 typename T1::elem_type& out_val,
|
|
max@0
|
26 typename T1::pod_type& out_sign,
|
|
max@0
|
27 const Base<typename T1::elem_type,T1>& X,
|
|
max@0
|
28 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
|
|
max@0
|
29 )
|
|
max@0
|
30 {
|
|
max@0
|
31 arma_extra_debug_sigprint();
|
|
max@0
|
32 arma_ignore(junk);
|
|
max@0
|
33
|
|
max@0
|
34 return auxlib::log_det(out_val, out_sign, X);
|
|
max@0
|
35 }
|
|
max@0
|
36
|
|
max@0
|
37
|
|
max@0
|
38
|
|
max@0
|
39 template<typename T1>
|
|
max@0
|
40 inline
|
|
max@0
|
41 void
|
|
max@0
|
42 log_det
|
|
max@0
|
43 (
|
|
max@0
|
44 typename T1::elem_type& out_val,
|
|
max@0
|
45 typename T1::pod_type& out_sign,
|
|
max@0
|
46 const Op<T1,op_diagmat>& X,
|
|
max@0
|
47 const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
|
|
max@0
|
48 )
|
|
max@0
|
49 {
|
|
max@0
|
50 arma_extra_debug_sigprint();
|
|
max@0
|
51 arma_ignore(junk);
|
|
max@0
|
52
|
|
max@0
|
53 typedef typename T1::elem_type eT;
|
|
max@0
|
54 typedef typename T1::pod_type T;
|
|
max@0
|
55
|
|
max@0
|
56 const diagmat_proxy<T1> A(X.m);
|
|
max@0
|
57
|
|
max@0
|
58 const uword N = A.n_elem;
|
|
max@0
|
59
|
|
max@0
|
60 if(N == 0)
|
|
max@0
|
61 {
|
|
max@0
|
62 out_val = eT(0);
|
|
max@0
|
63 out_sign = T(1);
|
|
max@0
|
64
|
|
max@0
|
65 return;
|
|
max@0
|
66 }
|
|
max@0
|
67
|
|
max@0
|
68 const eT x = A[0];
|
|
max@0
|
69
|
|
max@0
|
70 T sign = (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
|
|
max@0
|
71 eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
|
|
max@0
|
72
|
|
max@0
|
73 for(uword i=1; i<N; ++i)
|
|
max@0
|
74 {
|
|
max@0
|
75 const eT x = A[i];
|
|
max@0
|
76
|
|
max@0
|
77 sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1;
|
|
max@0
|
78 val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x);
|
|
max@0
|
79 }
|
|
max@0
|
80
|
|
max@0
|
81 out_val = val;
|
|
max@0
|
82 out_sign = sign;
|
|
max@0
|
83 }
|
|
max@0
|
84
|
|
max@0
|
85
|
|
max@0
|
86
|
|
max@0
|
87 //! @}
|
