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_powmat_meat.hpp | 261 +++++++++++++++++++++ 1 file changed, 261 insertions(+) create mode 100644 src/armadillo/include/armadillo_bits/op_powmat_meat.hpp (limited to 'src/armadillo/include/armadillo_bits/op_powmat_meat.hpp') diff --git a/src/armadillo/include/armadillo_bits/op_powmat_meat.hpp b/src/armadillo/include/armadillo_bits/op_powmat_meat.hpp new file mode 100644 index 0000000..323db32 --- /dev/null +++ b/src/armadillo/include/armadillo_bits/op_powmat_meat.hpp @@ -0,0 +1,261 @@ +// 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_powmat +//! @{ + + +template +inline +void +op_powmat::apply(Mat& out, const Op& expr) + { + arma_extra_debug_sigprint(); + + const uword y = expr.aux_uword_a; + const bool y_neg = (expr.aux_uword_b == uword(1)); + + const bool status = op_powmat::apply_direct(out, expr.m, y, y_neg); + + if(status == false) + { + out.soft_reset(); + arma_stop_runtime_error("powmat(): transformation failed"); + } + } + + + +template +inline +bool +op_powmat::apply_direct(Mat& out, const Base& X, const uword y, const bool y_neg) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + if(y_neg) + { + if(y == uword(1)) + { + return op_inv_gen_default::apply_direct(out, X.get_ref(), "powmat()"); + } + else + { + Mat X_inv; + + const bool inv_status = op_inv_gen_default::apply_direct(X_inv, X.get_ref(), "powmat()"); + + if(inv_status == false) { return false; } + + op_powmat::apply_direct_positive(out, X_inv, y); + } + } + else + { + const quasi_unwrap U(X.get_ref()); + + arma_debug_check( (U.M.is_square() == false), "powmat(): given matrix must be square sized" ); + + op_powmat::apply_direct_positive(out, U.M, y); + } + + return true; + } + + + +template +inline +void +op_powmat::apply_direct_positive(Mat& out, const Mat& X, const uword y) + { + arma_extra_debug_sigprint(); + + const uword N = X.n_rows; + + if(y == uword(0)) { out.eye(N,N); return; } + if(y == uword(1)) { out = X; return; } + + if(X.is_diagmat()) + { + arma_extra_debug_print("op_powmat: detected diagonal matrix"); + + podarray tmp(N); // use temporary array in case we have aliasing + + for(uword i=0; i tmp = X*X; out = X*tmp; } + else if(y == uword(4)) { const Mat tmp = X*X; out = tmp*tmp; } + else if(y == uword(5)) { const Mat tmp = X*X; out = X*tmp*tmp; } + else + { + Mat tmp = X; + + out = X; + + uword z = y-1; + + while(z > 0) + { + if(z & 1) { out = tmp * out; } + + z /= uword(2); + + if(z > 0) { tmp = tmp * tmp; } + } + } + } + } + + + +template +inline +void +op_powmat_cx::apply(Mat< std::complex >& out, const mtOp,T1,op_powmat_cx>& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::pod_type in_T; + + const in_T y = std::real(expr.aux_out_eT); + + const bool status = op_powmat_cx::apply_direct(out, expr.m, y); + + if(status == false) + { + out.soft_reset(); + arma_stop_runtime_error("powmat(): transformation failed"); + } + } + + + +template +inline +bool +op_powmat_cx::apply_direct(Mat< std::complex >& out, const Base& X, const typename T1::pod_type y) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type in_eT; + typedef typename T1::pod_type in_T; + typedef std::complex out_eT; + + if( y == in_T(int(y)) ) + { + arma_extra_debug_print("op_powmat_cx::apply_direct(): integer exponent detected; redirecting to op_powmat"); + + const uword y_val = (y < int(0)) ? uword(-y) : uword(y); + const bool y_neg = (y < int(0)); + + Mat tmp; + + const bool status = op_powmat::apply_direct(tmp, X.get_ref(), y_val, y_neg); + + if(status == false) { return false; } + + out = conv_to< Mat >::from(tmp); + + return true; + } + + const quasi_unwrap U(X.get_ref()); + const Mat& A = U.M; + + arma_debug_check( (A.is_square() == false), "powmat(): given matrix must be square sized" ); + + const uword N = A.n_rows; + + if(A.is_diagmat()) + { + arma_extra_debug_print("op_powmat_cx: detected diagonal matrix"); + + podarray tmp(N); // use temporary array in case we have aliasing + + for(uword i=0; i(A.at(i,i)), y) ; } + + out.zeros(N,N); + + for(uword i=0; i eigval; + Mat eigvec; + + const bool eig_status = eig_sym(eigval, eigvec, A); + + if(eig_status) + { + eigval = pow(eigval, y); + + const Mat tmp = diagmat(eigval) * eigvec.t(); + + out = conv_to< Mat >::from(eigvec * tmp); + + return true; + } + + arma_extra_debug_print("op_powmat_cx: sympd optimisation failed"); + + // fallthrough if optimisation failed + } + + bool powmat_status = false; + + Col eigval; + Mat eigvec; + + const bool eig_status = eig_gen(eigval, eigvec, A); + + if(eig_status) + { + eigval = pow(eigval, y); + + Mat eigvec_t = trans(eigvec); + Mat tmp = diagmat(conj(eigval)) * eigvec_t; + + const bool solve_status = auxlib::solve_square_fast(out, eigvec_t, tmp); + + if(solve_status) { out = trans(out); powmat_status = true; } + } + + return powmat_status; + } + + + +//! @} -- cgit v1.2.1