"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "include/armadillo_bits/op_expmat_meat.hpp" between
armadillo-10.8.2.tar.xz and armadillo-11.0.0.tar.xz

About: Armadillo is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use.

op_expmat_meat.hpp  (armadillo-10.8.2.tar.xz):op_expmat_meat.hpp  (armadillo-11.0.0.tar.xz)
skipping to change at line 62 skipping to change at line 62
if(is_op_diagmat<T1>::value) if(is_op_diagmat<T1>::value)
{ {
out = expr.get_ref(); // force the evaluation of diagmat() out = expr.get_ref(); // force the evaluation of diagmat()
arma_debug_check( (out.is_square() == false), "expmat(): given matrix must b e square sized" ); arma_debug_check( (out.is_square() == false), "expmat(): given matrix must b e square sized" );
const uword N = (std::min)(out.n_rows, out.n_cols); const uword N = (std::min)(out.n_rows, out.n_cols);
for(uword i=0; i<N; ++i) { out.at(i,i) = std::exp( out.at(i,i) ); } for(uword i=0; i<N; ++i) { out.at(i,i) = std::exp( out.at(i,i) ); }
return true;
} }
else
{
Mat<eT> A = expr.get_ref();
arma_debug_check( (A.is_square() == false), "expmat(): given matrix must be square sized" ); Mat<eT> A = expr.get_ref();
if(A.is_diagmat()) arma_debug_check( (A.is_square() == false), "expmat(): given matrix must be sq
{ uare sized" );
arma_extra_debug_print("op_expmat: detected diagonal matrix");
const uword N = (std::min)(A.n_rows, A.n_cols); if(A.is_diagmat())
{
arma_extra_debug_print("op_expmat: detected diagonal matrix");
out.zeros(N,N); const uword N = (std::min)(A.n_rows, A.n_cols);
for(uword i=0; i<N; ++i) { out.at(i,i) = std::exp( A.at(i,i) ); } out.zeros(N,N);
return true; for(uword i=0; i<N; ++i) { out.at(i,i) = std::exp( A.at(i,i) ); }
}
#if defined(ARMA_OPTIMISE_SYMPD) return true;
const bool try_sympd = sympd_helper::guess_sympd(A); }
#else
const bool try_sympd = false;
#endif
if(try_sympd)
{
arma_extra_debug_print("op_expmat: attempting sympd optimisation");
// if matrix A is sympd, all its eigenvalues are positive const bool try_sympd = arma_config::optimise_sympd && sympd_helper::guess_symp d(A);
Col< T> eigval; if(try_sympd)
Mat<eT> eigvec; {
arma_extra_debug_print("op_expmat: attempting sympd optimisation");
const bool eig_status = eig_sym_helper(eigval, eigvec, A, 'd', "expmat()") ; // if matrix A is sympd, all its eigenvalues are positive
if(eig_status) Col< T> eigval;
{ Mat<eT> eigvec;
eigval = exp(eigval);
out = eigvec * diagmat(eigval) * eigvec.t(); const bool eig_status = eig_sym_helper(eigval, eigvec, A, 'd', "expmat()");
return true; if(eig_status)
} {
eigval = exp(eigval);
arma_extra_debug_print("op_expmat: sympd optimisation failed"); out = eigvec * diagmat(eigval) * eigvec.t();
// fallthrough if eigen decomposition failed return true;
} }
const T norm_val = arma::norm(A, "inf"); arma_extra_debug_print("op_expmat: sympd optimisation failed");
const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) // fallthrough if eigen decomposition failed
: double(0); }
int exponent = int(0); std::frexp(log2_val, &exponent); const T norm_val = arma::norm(A, "inf");
const uword s = uword( (std::max)(int(0), exponent + int(1)) ); const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) : double(0);
A /= eT(eop_aux::pow(double(2), double(s))); int exponent = int(0); std::frexp(log2_val, &exponent);
T c = T(0.5); const uword s = uword( (std::max)(int(0), exponent + int(1)) );
Mat<eT> E(A.n_rows, A.n_rows, fill::eye); E += c * A; A /= eT(eop_aux::pow(double(2), double(s)));
Mat<eT> D(A.n_rows, A.n_rows, fill::eye); D -= c * A;
Mat<eT> X = A; T c = T(0.5);
bool positive = true; Mat<eT> E(A.n_rows, A.n_rows, fill::eye); E += c * A;
Mat<eT> D(A.n_rows, A.n_rows, fill::eye); D -= c * A;
const uword N = 6; Mat<eT> X = A;
for(uword i = 2; i <= N; ++i) bool positive = true;
{
c = c * T(N - i + 1) / T(i * (2*N - i + 1));
X = A * X; const uword N = 6;
E += c * X; for(uword i = 2; i <= N; ++i)
{
c = c * T(N - i + 1) / T(i * (2*N - i + 1));
if(positive) { D += c * X; } else { D -= c * X; } X = A * X;
positive = (positive) ? false : true; E += c * X;
}
if( (D.is_finite() == false) || (E.is_finite() == false) ) { return false; } if(positive) { D += c * X; } else { D -= c * X; }
const bool status = solve(out, D, E, solve_opts::no_approx); positive = (positive) ? false : true;
}
if(status == false) { return false; } if( (D.is_finite() == false) || (E.is_finite() == false) ) { return false; }
for(uword i=0; i < s; ++i) { out = out * out; } const bool status = solve(out, D, E, solve_opts::no_approx);
}
if(status == false) { return false; }
for(uword i=0; i < s; ++i) { out = out * out; }
return true; return true;
} }
template<typename T1> template<typename T1>
inline inline
void void
op_expmat_sym::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_expmat_sym >& in) op_expmat_sym::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_expmat_sym >& in)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
skipping to change at line 184 skipping to change at line 180
template<typename T1> template<typename T1>
inline inline
bool bool
op_expmat_sym::apply_direct(Mat<typename T1::elem_type>& out, const Base<typenam e T1::elem_type,T1>& expr) op_expmat_sym::apply_direct(Mat<typename T1::elem_type>& out, const Base<typenam e T1::elem_type,T1>& expr)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
#if defined(ARMA_USE_LAPACK) #if defined(ARMA_USE_LAPACK)
{ {
typedef typename T1::pod_type T;
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
typedef typename T1::pod_type T;
const unwrap<T1> U(expr.get_ref()); const unwrap<T1> U(expr.get_ref());
const Mat<eT>& X = U.M; const Mat<eT>& X = U.M;
arma_debug_check( (X.is_square() == false), "expmat_sym(): given matrix must be square sized" ); arma_debug_check( (X.is_square() == false), "expmat_sym(): given matrix must be square sized" );
if((arma_config::debug) && (arma_config::warn_level > 0) && (is_cx<eT>::yes)
&& (sympd_helper::check_diag_imag(X) == false))
{
arma_debug_warn_level(1, "inv_sympd(): imaginary components on diagonal ar
e non-zero");
}
if(is_op_diagmat<T1>::value || X.is_diagmat())
{
arma_extra_debug_print("op_expmat_sym: detected diagonal matrix");
out = X;
eT* colmem = out.memptr();
const uword N = X.n_rows;
for(uword i=0; i<N; ++i)
{
eT& out_ii = colmem[i];
T out_ii_real = access::tmp_real(out_ii);
out_ii = eT( std::exp(out_ii_real) );
colmem += N;
}
return true;
}
Col< T> eigval; Col< T> eigval;
Mat<eT> eigvec; Mat<eT> eigvec;
const bool status = eig_sym_helper(eigval, eigvec, X, 'd', "expmat_sym()"); const bool status = eig_sym_helper(eigval, eigvec, X, 'd', "expmat_sym()");
if(status == false) { return false; } if(status == false) { return false; }
eigval = exp(eigval); eigval = exp(eigval);
out = eigvec * diagmat(eigval) * eigvec.t(); out = eigvec * diagmat(eigval) * eigvec.t();
 End of changes. 39 change blocks. 
59 lines changed or deleted 85 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)