Chris@49
|
1 // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2008-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 op_inv
|
Chris@49
|
10 //! @{
|
Chris@49
|
11
|
Chris@49
|
12
|
Chris@49
|
13 //! immediate inverse of a matrix, storing the result in a dense matrix
|
Chris@49
|
14 template<typename eT>
|
Chris@49
|
15 inline
|
Chris@49
|
16 void
|
Chris@49
|
17 op_inv::apply(Mat<eT>& out, const Mat<eT>& A, const bool slow)
|
Chris@49
|
18 {
|
Chris@49
|
19 arma_extra_debug_sigprint();
|
Chris@49
|
20
|
Chris@49
|
21 // no need to check for aliasing, due to:
|
Chris@49
|
22 // - auxlib::inv() copies A to out before inversion
|
Chris@49
|
23 // - for 2x2 and 3x3 matrices the code is alias safe
|
Chris@49
|
24
|
Chris@49
|
25 bool status = auxlib::inv(out, A, slow);
|
Chris@49
|
26
|
Chris@49
|
27 if(status == false)
|
Chris@49
|
28 {
|
Chris@49
|
29 out.reset();
|
Chris@49
|
30 arma_bad("inv(): matrix appears to be singular");
|
Chris@49
|
31 }
|
Chris@49
|
32 }
|
Chris@49
|
33
|
Chris@49
|
34
|
Chris@49
|
35
|
Chris@49
|
36 //! immediate inverse of T1, storing the result in a dense matrix
|
Chris@49
|
37 template<typename T1>
|
Chris@49
|
38 inline
|
Chris@49
|
39 void
|
Chris@49
|
40 op_inv::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv>& X)
|
Chris@49
|
41 {
|
Chris@49
|
42 arma_extra_debug_sigprint();
|
Chris@49
|
43
|
Chris@49
|
44 const strip_diagmat<T1> strip(X.m);
|
Chris@49
|
45
|
Chris@49
|
46 if(strip.do_diagmat == true)
|
Chris@49
|
47 {
|
Chris@49
|
48 op_inv::apply_diag(out, strip.M);
|
Chris@49
|
49 }
|
Chris@49
|
50 else
|
Chris@49
|
51 {
|
Chris@49
|
52 const uword mode = X.aux_uword_a;
|
Chris@49
|
53
|
Chris@49
|
54 const bool status = (mode == 0) ? auxlib::inv(out, X.m) : auxlib::inv(out, X.m, true);
|
Chris@49
|
55
|
Chris@49
|
56 if(status == false)
|
Chris@49
|
57 {
|
Chris@49
|
58 out.reset();
|
Chris@49
|
59 arma_bad("inv(): matrix appears to be singular");
|
Chris@49
|
60 }
|
Chris@49
|
61 }
|
Chris@49
|
62 }
|
Chris@49
|
63
|
Chris@49
|
64
|
Chris@49
|
65
|
Chris@49
|
66 template<typename T1>
|
Chris@49
|
67 inline
|
Chris@49
|
68 void
|
Chris@49
|
69 op_inv::apply_diag(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type, T1>& X)
|
Chris@49
|
70 {
|
Chris@49
|
71 arma_extra_debug_sigprint();
|
Chris@49
|
72
|
Chris@49
|
73 typedef typename T1::elem_type eT;
|
Chris@49
|
74
|
Chris@49
|
75 const diagmat_proxy_check<T1> A(X.get_ref(), out);
|
Chris@49
|
76
|
Chris@49
|
77 const uword N = A.n_elem;
|
Chris@49
|
78
|
Chris@49
|
79 out.set_size(N,N);
|
Chris@49
|
80
|
Chris@49
|
81 for(uword col=0; col<N; ++col)
|
Chris@49
|
82 {
|
Chris@49
|
83 for(uword row=0; row<col; ++row) { out.at(row,col) = eT(0); }
|
Chris@49
|
84
|
Chris@49
|
85 out.at(col,col) = eT(1) / A[col];
|
Chris@49
|
86
|
Chris@49
|
87 for(uword row=col+1; row<N; ++row) { out.at(row,col) = eT(0); }
|
Chris@49
|
88 }
|
Chris@49
|
89
|
Chris@49
|
90 }
|
Chris@49
|
91
|
Chris@49
|
92
|
Chris@49
|
93
|
Chris@49
|
94 //! inverse of T1 (triangular matrices)
|
Chris@49
|
95 template<typename T1>
|
Chris@49
|
96 inline
|
Chris@49
|
97 void
|
Chris@49
|
98 op_inv_tr::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_tr>& X)
|
Chris@49
|
99 {
|
Chris@49
|
100 arma_extra_debug_sigprint();
|
Chris@49
|
101
|
Chris@49
|
102 const bool status = auxlib::inv_tr(out, X.m, X.aux_uword_a);
|
Chris@49
|
103
|
Chris@49
|
104 if(status == false)
|
Chris@49
|
105 {
|
Chris@49
|
106 out.reset();
|
Chris@49
|
107 arma_bad("inv(): matrix appears to be singular");
|
Chris@49
|
108 }
|
Chris@49
|
109 }
|
Chris@49
|
110
|
Chris@49
|
111
|
Chris@49
|
112
|
Chris@49
|
113 //! inverse of T1 (symmetric positive definite matrices)
|
Chris@49
|
114 template<typename T1>
|
Chris@49
|
115 inline
|
Chris@49
|
116 void
|
Chris@49
|
117 op_inv_sympd::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_sympd>& X)
|
Chris@49
|
118 {
|
Chris@49
|
119 arma_extra_debug_sigprint();
|
Chris@49
|
120
|
Chris@49
|
121 const bool status = auxlib::inv_sympd(out, X.m, X.aux_uword_a);
|
Chris@49
|
122
|
Chris@49
|
123 if(status == false)
|
Chris@49
|
124 {
|
Chris@49
|
125 out.reset();
|
Chris@49
|
126 arma_bad("inv(): matrix appears to be singular");
|
Chris@49
|
127 }
|
Chris@49
|
128 }
|
Chris@49
|
129
|
Chris@49
|
130
|
Chris@49
|
131
|
Chris@49
|
132 //! @}
|