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 --- src/armadillo/include/armadillo_bits/fn_trace.hpp | 663 ++++++++++++++++++++++ 1 file changed, 663 insertions(+) create mode 100644 src/armadillo/include/armadillo_bits/fn_trace.hpp (limited to 'src/armadillo/include/armadillo_bits/fn_trace.hpp') diff --git a/src/armadillo/include/armadillo_bits/fn_trace.hpp b/src/armadillo/include/armadillo_bits/fn_trace.hpp new file mode 100644 index 0000000..8a15bac --- /dev/null +++ b/src/armadillo/include/armadillo_bits/fn_trace.hpp @@ -0,0 +1,663 @@ +// 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 fn_trace +//! @{ + + +template +arma_warn_unused +inline +typename T1::elem_type +trace(const Base& X) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const Proxy P(X.get_ref()); + + const uword N = (std::min)(P.get_n_rows(), P.get_n_cols()); + + eT val1 = eT(0); + eT val2 = eT(0); + + uword i,j; + for(i=0, j=1; j +arma_warn_unused +inline +typename T1::elem_type +trace(const Op& X) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const diagmat_proxy A(X.m); + + const uword N = (std::min)(A.n_rows, A.n_cols); + + eT val = eT(0); + + for(uword i=0; i +arma_warn_unused +inline +typename enable_if2< is_cx::no, typename T1::elem_type>::result +trace(const Glue& X) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const partial_unwrap tmp1(X.A); + const partial_unwrap tmp2(X.B); + + const typename partial_unwrap::stored_type& A = tmp1.M; + const typename partial_unwrap::stored_type& B = tmp2.M; + + const bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times; + const eT alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0); + + arma_debug_assert_trans_mul_size< partial_unwrap::do_trans, partial_unwrap::do_trans >(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "matrix multiplication"); + + if( (A.n_elem == 0) || (B.n_elem == 0) ) { return eT(0); } + + const uword A_n_rows = A.n_rows; + const uword A_n_cols = A.n_cols; + + const uword B_n_rows = B.n_rows; + const uword B_n_cols = B.n_cols; + + eT acc = eT(0); + + if( (partial_unwrap::do_trans == false) && (partial_unwrap::do_trans == false) ) + { + const uword N = (std::min)(A_n_rows, B_n_cols); + + eT acc1 = eT(0); + eT acc2 = eT(0); + + for(uword k=0; k < N; ++k) + { + const eT* B_colptr = B.colptr(k); + + // condition: A_n_cols = B_n_rows + + uword j; + + for(j=1; j < A_n_cols; j+=2) + { + const uword i = (j-1); + + const eT tmp_i = B_colptr[i]; + const eT tmp_j = B_colptr[j]; + + acc1 += A.at(k, i) * tmp_i; + acc2 += A.at(k, j) * tmp_j; + } + + const uword i = (j-1); + + if(i < A_n_cols) + { + acc1 += A.at(k, i) * B_colptr[i]; + } + } + + acc = (acc1 + acc2); + } + else + if( (partial_unwrap::do_trans == true ) && (partial_unwrap::do_trans == false) ) + { + const uword N = (std::min)(A_n_cols, B_n_cols); + + for(uword k=0; k < N; ++k) + { + const eT* A_colptr = A.colptr(k); + const eT* B_colptr = B.colptr(k); + + // condition: A_n_rows = B_n_rows + acc += op_dot::direct_dot(A_n_rows, A_colptr, B_colptr); + } + } + else + if( (partial_unwrap::do_trans == false) && (partial_unwrap::do_trans == true ) ) + { + const uword N = (std::min)(A_n_rows, B_n_rows); + + for(uword k=0; k < N; ++k) + { + // condition: A_n_cols = B_n_cols + for(uword i=0; i < A_n_cols; ++i) + { + acc += A.at(k,i) * B.at(k,i); + } + } + } + else + if( (partial_unwrap::do_trans == true ) && (partial_unwrap::do_trans == true ) ) + { + const uword N = (std::min)(A_n_cols, B_n_rows); + + for(uword k=0; k < N; ++k) + { + const eT* A_colptr = A.colptr(k); + + // condition: A_n_rows = B_n_cols + for(uword i=0; i < A_n_rows; ++i) + { + acc += A_colptr[i] * B.at(k,i); + } + } + } + + return (use_alpha) ? (alpha * acc) : acc; + } + + + +//! speedup for trace(A*B); complex elements +template +arma_warn_unused +inline +typename enable_if2< is_cx::yes, typename T1::elem_type>::result +trace(const Glue& X) + { + arma_extra_debug_sigprint(); + + typedef typename T1::pod_type T; + typedef typename T1::elem_type eT; + + const partial_unwrap tmp1(X.A); + const partial_unwrap tmp2(X.B); + + const typename partial_unwrap::stored_type& A = tmp1.M; + const typename partial_unwrap::stored_type& B = tmp2.M; + + const bool use_alpha = partial_unwrap::do_times || partial_unwrap::do_times; + const eT alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0); + + arma_debug_assert_trans_mul_size< partial_unwrap::do_trans, partial_unwrap::do_trans >(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "matrix multiplication"); + + if( (A.n_elem == 0) || (B.n_elem == 0) ) { return eT(0); } + + const uword A_n_rows = A.n_rows; + const uword A_n_cols = A.n_cols; + + const uword B_n_rows = B.n_rows; + const uword B_n_cols = B.n_cols; + + eT acc = eT(0); + + if( (partial_unwrap::do_trans == false) && (partial_unwrap::do_trans == false) ) + { + const uword N = (std::min)(A_n_rows, B_n_cols); + + T acc_real = T(0); + T acc_imag = T(0); + + for(uword k=0; k < N; ++k) + { + const eT* B_colptr = B.colptr(k); + + // condition: A_n_cols = B_n_rows + + for(uword i=0; i < A_n_cols; ++i) + { + // acc += A.at(k, i) * B_colptr[i]; + + const std::complex& xx = A.at(k, i); + const std::complex& yy = B_colptr[i]; + + const T a = xx.real(); + const T b = xx.imag(); + + const T c = yy.real(); + const T d = yy.imag(); + + acc_real += (a*c) - (b*d); + acc_imag += (a*d) + (b*c); + } + } + + acc = std::complex(acc_real, acc_imag); + } + else + if( (partial_unwrap::do_trans == true) && (partial_unwrap::do_trans == false) ) + { + const uword N = (std::min)(A_n_cols, B_n_cols); + + T acc_real = T(0); + T acc_imag = T(0); + + for(uword k=0; k < N; ++k) + { + const eT* A_colptr = A.colptr(k); + const eT* B_colptr = B.colptr(k); + + // condition: A_n_rows = B_n_rows + + for(uword i=0; i < A_n_rows; ++i) + { + // acc += std::conj(A_colptr[i]) * B_colptr[i]; + + const std::complex& xx = A_colptr[i]; + const std::complex& yy = B_colptr[i]; + + const T a = xx.real(); + const T b = xx.imag(); + + const T c = yy.real(); + const T d = yy.imag(); + + // take into account the complex conjugate of xx + + acc_real += (a*c) + (b*d); + acc_imag += (a*d) - (b*c); + } + } + + acc = std::complex(acc_real, acc_imag); + } + else + if( (partial_unwrap::do_trans == false) && (partial_unwrap::do_trans == true) ) + { + const uword N = (std::min)(A_n_rows, B_n_rows); + + T acc_real = T(0); + T acc_imag = T(0); + + for(uword k=0; k < N; ++k) + { + // condition: A_n_cols = B_n_cols + for(uword i=0; i < A_n_cols; ++i) + { + // acc += A.at(k,i) * std::conj(B.at(k,i)); + + const std::complex& xx = A.at(k, i); + const std::complex& yy = B.at(k, i); + + const T a = xx.real(); + const T b = xx.imag(); + + const T c = yy.real(); + const T d = -yy.imag(); // take the conjugate + + acc_real += (a*c) - (b*d); + acc_imag += (a*d) + (b*c); + } + } + + acc = std::complex(acc_real, acc_imag); + } + else + if( (partial_unwrap::do_trans == true) && (partial_unwrap::do_trans == true) ) + { + const uword N = (std::min)(A_n_cols, B_n_rows); + + T acc_real = T(0); + T acc_imag = T(0); + + for(uword k=0; k < N; ++k) + { + const eT* A_colptr = A.colptr(k); + + // condition: A_n_rows = B_n_cols + for(uword i=0; i < A_n_rows; ++i) + { + // acc += std::conj(A_colptr[i]) * std::conj(B.at(k,i)); + + const std::complex& xx = A_colptr[i]; + const std::complex& yy = B.at(k, i); + + const T a = xx.real(); + const T b = -xx.imag(); // take the conjugate + + const T c = yy.real(); + const T d = -yy.imag(); // take the conjugate + + acc_real += (a*c) - (b*d); + acc_imag += (a*d) + (b*c); + } + } + + acc = std::complex(acc_real, acc_imag); + } + + return (use_alpha) ? eT(alpha * acc) : eT(acc); + } + + + +//! trace of sparse object; generic version +template +arma_warn_unused +inline +typename T1::elem_type +trace(const SpBase& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const SpProxy P(expr.get_ref()); + + const uword N = (std::min)(P.get_n_rows(), P.get_n_cols()); + + eT acc = eT(0); + + if( (is_SpMat::stored_type>::value) && (P.get_n_nonzero() >= 5*N) ) + { + const unwrap_spmat::stored_type> U(P.Q); + + const SpMat& X = U.M; + + for(uword i=0; i < N; ++i) + { + acc += X.at(i,i); // use binary search + } + } + else + { + typename SpProxy::const_iterator_type it = P.begin(); + + const uword P_n_nz = P.get_n_nonzero(); + + for(uword i=0; i < P_n_nz; ++i) + { + if(it.row() == it.col()) { acc += (*it); } + + ++it; + } + } + + return acc; + } + + + +//! trace of sparse object; speedup for trace(A + B) +template +arma_warn_unused +inline +typename T1::elem_type +trace(const SpGlue& expr) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat UA(expr.A); + const unwrap_spmat UB(expr.B); + + arma_debug_assert_same_size(UA.M.n_rows, UA.M.n_cols, UB.M.n_rows, UB.M.n_cols, "addition"); + + return (trace(UA.M) + trace(UB.M)); + } + + + +//! trace of sparse object; speedup for trace(A - B) +template +arma_warn_unused +inline +typename T1::elem_type +trace(const SpGlue& expr) + { + arma_extra_debug_sigprint(); + + const unwrap_spmat UA(expr.A); + const unwrap_spmat UB(expr.B); + + arma_debug_assert_same_size(UA.M.n_rows, UA.M.n_cols, UB.M.n_rows, UB.M.n_cols, "subtraction"); + + return (trace(UA.M) - trace(UB.M)); + } + + + +//! trace of sparse object; speedup for trace(A % B) +template +arma_warn_unused +inline +typename T1::elem_type +trace(const SpGlue& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const unwrap_spmat UA(expr.A); + const unwrap_spmat UB(expr.B); + + const SpMat& A = UA.M; + const SpMat& B = UB.M; + + arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "element-wise multiplication"); + + const uword N = (std::min)(A.n_rows, A.n_cols); + + eT acc = eT(0); + + for(uword i=0; i +arma_warn_unused +inline +typename T1::elem_type +trace(const SpGlue& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + // better-than-nothing implementation + + const unwrap_spmat UA(expr.A); + const unwrap_spmat UB(expr.B); + + const SpMat& A = UA.M; + const SpMat& B = UB.M; + + arma_debug_assert_mul_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "matrix multiplication"); + + if( (A.n_nonzero == 0) || (B.n_nonzero == 0) ) { return eT(0); } + + const uword N = (std::min)(A.n_rows, B.n_cols); + + eT acc = eT(0); + + // TODO: the threshold may need tuning for complex matrices + if( (A.n_nonzero >= 5*N) || (B.n_nonzero >= 5*N) ) + { + for(uword k=0; k < N; ++k) + { + typename SpMat::const_col_iterator B_it = B.begin_col_no_sync(k); + typename SpMat::const_col_iterator B_it_end = B.end_col_no_sync(k); + + while(B_it != B_it_end) + { + const eT B_val = (*B_it); + const uword i = B_it.row(); + + acc += A.at(k,i) * B_val; + + ++B_it; + } + } + } + else + { + const SpMat AB = A * B; + + acc = trace(AB); + } + + return acc; + } + + + +//! trace of sparse object; speedup for trace(A.t()*B); non-complex elements +template +arma_warn_unused +inline +typename enable_if2< is_cx::no, typename T1::elem_type>::result +trace(const SpGlue, T2, spglue_times>& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const unwrap_spmat UA(expr.A.m); + const unwrap_spmat UB(expr.B); + + const SpMat& A = UA.M; + const SpMat& B = UB.M; + + // NOTE: deliberately swapped A.n_rows and A.n_cols to take into account the requested transpose operation + arma_debug_assert_mul_size(A.n_cols, A.n_rows, B.n_rows, B.n_cols, "matrix multiplication"); + + if( (A.n_nonzero == 0) || (B.n_nonzero == 0) ) { return eT(0); } + + const uword N = (std::min)(A.n_cols, B.n_cols); + + eT acc = eT(0); + + if( (A.n_nonzero >= 5*N) || (B.n_nonzero >= 5*N) ) + { + for(uword k=0; k < N; ++k) + { + typename SpMat::const_col_iterator B_it = B.begin_col_no_sync(k); + typename SpMat::const_col_iterator B_it_end = B.end_col_no_sync(k); + + while(B_it != B_it_end) + { + const eT B_val = (*B_it); + const uword i = B_it.row(); + + acc += A.at(i,k) * B_val; + + ++B_it; + } + } + } + else + { + const SpMat AtB = A.t() * B; + + acc = trace(AtB); + } + + return acc; + } + + + +//! trace of sparse object; speedup for trace(A.t()*B); complex elements +template +arma_warn_unused +inline +typename enable_if2< is_cx::yes, typename T1::elem_type>::result +trace(const SpGlue, T2, spglue_times>& expr) + { + arma_extra_debug_sigprint(); + + typedef typename T1::elem_type eT; + + const unwrap_spmat UA(expr.A.m); + const unwrap_spmat UB(expr.B); + + const SpMat& A = UA.M; + const SpMat& B = UB.M; + + // NOTE: deliberately swapped A.n_rows and A.n_cols to take into account the requested transpose operation + arma_debug_assert_mul_size(A.n_cols, A.n_rows, B.n_rows, B.n_cols, "matrix multiplication"); + + if( (A.n_nonzero == 0) || (B.n_nonzero == 0) ) { return eT(0); } + + const uword N = (std::min)(A.n_cols, B.n_cols); + + eT acc = eT(0); + + // TODO: the threshold may need tuning for complex matrices + if( (A.n_nonzero >= 5*N) || (B.n_nonzero >= 5*N) ) + { + for(uword k=0; k < N; ++k) + { + typename SpMat::const_col_iterator B_it = B.begin_col_no_sync(k); + typename SpMat::const_col_iterator B_it_end = B.end_col_no_sync(k); + + while(B_it != B_it_end) + { + const eT B_val = (*B_it); + const uword i = B_it.row(); + + acc += std::conj(A.at(i,k)) * B_val; + + ++B_it; + } + } + } + else + { + const SpMat AtB = A.t() * B; + + acc = trace(AtB); + } + + return acc; + } + + + +//! @} -- cgit v1.2.1