From eda5bc26f44ee9a6f83dcf8c91f17296d7fc509d Mon Sep 17 00:00:00 2001 From: Nao Pross Date: Mon, 12 Feb 2024 14:52:43 +0100 Subject: Move into version control --- .../include/armadillo_bits/op_expmat_meat.hpp | 256 +++++++++++++++++++++ 1 file changed, 256 insertions(+) create mode 100644 src/armadillo/include/armadillo_bits/op_expmat_meat.hpp (limited to 'src/armadillo/include/armadillo_bits/op_expmat_meat.hpp') diff --git a/src/armadillo/include/armadillo_bits/op_expmat_meat.hpp b/src/armadillo/include/armadillo_bits/op_expmat_meat.hpp new file mode 100644 index 0000000..d45fb36 --- /dev/null +++ b/src/armadillo/include/armadillo_bits/op_expmat_meat.hpp @@ -0,0 +1,256 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// ------------------------------------------------------------------------ + + + +//! \addtogroup op_expmat +//! @{ + + +//! implementation based on: +//! Cleve Moler, Charles Van Loan. +//! Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. +//! SIAM Review, Vol. 45, No. 1, 2003, pp. 3-49. +//! http://dx.doi.org/10.1137/S00361445024180 + + +template +inline +void +op_expmat::apply(Mat& out, const Op& expr) + { + arma_extra_debug_sigprint(); + + const bool status = op_expmat::apply_direct(out, expr.m); + + if(status == false) + { + out.soft_reset(); + arma_stop_runtime_error("expmat(): given matrix appears ill-conditioned"); + } + } + + + +template +inline +bool +op_expmat::apply_direct(Mat& out, const Base& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; + + if(is_op_diagmat::value) + { + out = expr.get_ref(); // force the evaluation of diagmat() + + arma_debug_check( (out.is_square() == false), "expmat(): given matrix must be square sized", [&](){ out.soft_reset(); } ); + + const uword N = (std::min)(out.n_rows, out.n_cols); + + for(uword i=0; i A = expr.get_ref(); + + arma_debug_check( (A.is_square() == false), "expmat(): given matrix must be square sized" ); + + if(A.is_diagmat()) + { + arma_extra_debug_print("op_expmat: detected diagonal matrix"); + + const uword N = (std::min)(A.n_rows, A.n_cols); + + out.zeros(N,N); + + for(uword i=0; i::no) ? (is_approx_sym) : (is_approx_sym && is_approx_sympd)); + } + + if(do_sym) + { + arma_extra_debug_print("op_expmat: symmetric/hermitian optimisation"); + + Col< T> eigval; + Mat eigvec; + + const bool eig_status = eig_sym_helper(eigval, eigvec, A, 'd', "expmat()"); + + if(eig_status == false) { return false; } + + eigval = exp(eigval); + + out = eigvec * diagmat(eigval) * eigvec.t(); + + return true; + } + + const T norm_val = arma::norm(A, "inf"); + + if(arma_isfinite(norm_val) == false) { return false; } + + const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) : double(0); + + int exponent = int(0); std::frexp(log2_val, &exponent); + + const uword s = uword( (std::max)(int(0), exponent + int(1)) ); + + A /= eT(eop_aux::pow(double(2), double(s))); + + T c = T(0.5); + + Mat E(A.n_rows, A.n_rows, fill::eye); E += c * A; + Mat D(A.n_rows, A.n_rows, fill::eye); D -= c * A; + + Mat X = A; + + bool positive = true; + + const uword N = 6; + + for(uword i = 2; i <= N; ++i) + { + c = c * T(N - i + 1) / T(i * (2*N - i + 1)); + + X = A * X; + + E += c * X; + + if(positive) { D += c * X; } else { D -= c * X; } + + positive = (positive) ? false : true; + } + + if( (D.internal_has_nonfinite()) || (E.internal_has_nonfinite()) ) { return false; } + + 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; + } + + + +template +inline +void +op_expmat_sym::apply(Mat& out, const Op& in) + { + arma_extra_debug_sigprint(); + + const bool status = op_expmat_sym::apply_direct(out, in.m); + + if(status == false) + { + out.soft_reset(); + arma_stop_runtime_error("expmat_sym(): transformation failed"); + } + } + + + +template +inline +bool +op_expmat_sym::apply_direct(Mat& out, const Base& expr) + { + arma_extra_debug_sigprint(); + + #if defined(ARMA_USE_LAPACK) + { + typedef typename T1::elem_type eT; + typedef typename T1::pod_type T; + + const unwrap U(expr.get_ref()); + const Mat& X = U.M; + + 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::yes) && (sym_helper::check_diag_imag(X) == false)) + { + arma_debug_warn_level(1, "inv_sympd(): imaginary components on diagonal are non-zero"); + } + + if(is_op_diagmat::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 eigval; + Mat eigvec; + + const bool status = eig_sym_helper(eigval, eigvec, X, 'd', "expmat_sym()"); + + if(status == false) { return false; } + + eigval = exp(eigval); + + out = eigvec * diagmat(eigval) * eigvec.t(); + + return true; + } + #else + { + arma_ignore(out); + arma_ignore(expr); + arma_stop_logic_error("expmat_sym(): use of LAPACK must be enabled"); + return false; + } + #endif + } + + + +//! @} -- cgit v1.2.1