1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19
20 //! \addtogroup op_expmat
21 //! @{
22
23
24 //! implementation based on:
25 //! Cleve Moler, Charles Van Loan.
26 //! Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later.
27 //! SIAM Review, Vol. 45, No. 1, 2003, pp. 3-49.
28 //! http://dx.doi.org/10.1137/S00361445024180
29
30
31 template<typename T1>
32 inline
33 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_expmat> & expr)34 op_expmat::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_expmat>& expr)
35 {
36 arma_extra_debug_sigprint();
37
38 const bool status = op_expmat::apply_direct(out, expr.m);
39
40 if(status == false)
41 {
42 out.soft_reset();
43 arma_stop_runtime_error("expmat(): given matrix appears ill-conditioned");
44 }
45 }
46
47
48
49 template<typename T1>
50 inline
51 bool
apply_direct(Mat<typename T1::elem_type> & out,const Base<typename T1::elem_type,T1> & expr)52 op_expmat::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type, T1>& expr)
53 {
54 arma_extra_debug_sigprint();
55
56 typedef typename T1::elem_type eT;
57 typedef typename T1::pod_type T;
58
59 if(is_op_diagmat<T1>::value)
60 {
61 out = expr.get_ref(); // force the evaluation of diagmat()
62
63 arma_debug_check( (out.is_square() == false), "expmat(): given matrix must be square sized" );
64
65 const uword N = (std::min)(out.n_rows, out.n_cols);
66
67 for(uword i=0; i<N; ++i) { out.at(i,i) = std::exp( out.at(i,i) ); }
68 }
69 else
70 {
71 Mat<eT> A = expr.get_ref();
72
73 arma_debug_check( (A.is_square() == false), "expmat(): given matrix must be square sized" );
74
75 if(A.is_diagmat())
76 {
77 const uword N = (std::min)(A.n_rows, A.n_cols);
78
79 out.zeros(N,N);
80
81 for(uword i=0; i<N; ++i) { out.at(i,i) = std::exp( A.at(i,i) ); }
82
83 return true;
84 }
85
86 #if defined(ARMA_OPTIMISE_SYMPD)
87 const bool try_sympd = sympd_helper::guess_sympd_anysize(A);
88 #else
89 const bool try_sympd = false;
90 #endif
91
92 if(try_sympd)
93 {
94 arma_extra_debug_print("op_expmat: attempting sympd optimisation");
95
96 // if matrix A is sympd, all its eigenvalues are positive
97
98 Col< T> eigval;
99 Mat<eT> eigvec;
100
101 const bool eig_status = eig_sym_helper(eigval, eigvec, A, 'd', "expmat()");
102
103 if(eig_status)
104 {
105 eigval = exp(eigval);
106
107 out = eigvec * diagmat(eigval) * eigvec.t();
108
109 return true;
110 }
111
112 arma_extra_debug_print("op_expmat: sympd optimisation failed");
113
114 // fallthrough if eigen decomposition failed
115 }
116
117 const T norm_val = arma::norm(A, "inf");
118
119 const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) : double(0);
120
121 int exponent = int(0); std::frexp(log2_val, &exponent);
122
123 const uword s = uword( (std::max)(int(0), exponent + int(1)) );
124
125 A /= eT(eop_aux::pow(double(2), double(s)));
126
127 T c = T(0.5);
128
129 Mat<eT> E(A.n_rows, A.n_rows, fill::eye); E += c * A;
130 Mat<eT> D(A.n_rows, A.n_rows, fill::eye); D -= c * A;
131
132 Mat<eT> X = A;
133
134 bool positive = true;
135
136 const uword N = 6;
137
138 for(uword i = 2; i <= N; ++i)
139 {
140 c = c * T(N - i + 1) / T(i * (2*N - i + 1));
141
142 X = A * X;
143
144 E += c * X;
145
146 if(positive) { D += c * X; } else { D -= c * X; }
147
148 positive = (positive) ? false : true;
149 }
150
151 if( (D.is_finite() == false) || (E.is_finite() == false) ) { return false; }
152
153 const bool status = solve(out, D, E, solve_opts::no_approx);
154
155 if(status == false) { return false; }
156
157 for(uword i=0; i < s; ++i) { out = out * out; }
158 }
159
160 return true;
161 }
162
163
164
165 template<typename T1>
166 inline
167 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_expmat_sym> & in)168 op_expmat_sym::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_expmat_sym>& in)
169 {
170 arma_extra_debug_sigprint();
171
172 const bool status = op_expmat_sym::apply_direct(out, in.m);
173
174 if(status == false)
175 {
176 out.soft_reset();
177 arma_stop_runtime_error("expmat_sym(): transformation failed");
178 }
179 }
180
181
182
183 template<typename T1>
184 inline
185 bool
apply_direct(Mat<typename T1::elem_type> & out,const Base<typename T1::elem_type,T1> & expr)186 op_expmat_sym::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1>& expr)
187 {
188 arma_extra_debug_sigprint();
189
190 #if defined(ARMA_USE_LAPACK)
191 {
192 typedef typename T1::pod_type T;
193 typedef typename T1::elem_type eT;
194
195 const unwrap<T1> U(expr.get_ref());
196 const Mat<eT>& X = U.M;
197
198 arma_debug_check( (X.is_square() == false), "expmat_sym(): given matrix must be square sized" );
199
200 Col< T> eigval;
201 Mat<eT> eigvec;
202
203 const bool status = eig_sym_helper(eigval, eigvec, X, 'd', "expmat_sym()");
204
205 if(status == false) { return false; }
206
207 eigval = exp(eigval);
208
209 out = eigvec * diagmat(eigval) * eigvec.t();
210
211 return true;
212 }
213 #else
214 {
215 arma_ignore(out);
216 arma_ignore(expr);
217 arma_stop_logic_error("expmat_sym(): use of LAPACK must be enabled");
218 return false;
219 }
220 #endif
221 }
222
223
224
225 //! @}
226