"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

---