From a556b45abf18f1bd509daaf63b66b7d55e9fd291 Mon Sep 17 00:00:00 2001 From: jjesswan Date: Mon, 22 Apr 2024 21:56:26 -0400 Subject: add engine version --- engine-ocean/Eigen/src/Core/SolveTriangular.h | 235 ++++++++++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100644 engine-ocean/Eigen/src/Core/SolveTriangular.h (limited to 'engine-ocean/Eigen/src/Core/SolveTriangular.h') diff --git a/engine-ocean/Eigen/src/Core/SolveTriangular.h b/engine-ocean/Eigen/src/Core/SolveTriangular.h new file mode 100644 index 0000000..dfbf995 --- /dev/null +++ b/engine-ocean/Eigen/src/Core/SolveTriangular.h @@ -0,0 +1,235 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SOLVETRIANGULAR_H +#define EIGEN_SOLVETRIANGULAR_H + +namespace Eigen { + +namespace internal { + +// Forward declarations: +// The following two routines are implemented in the products/TriangularSolver*.h files +template +struct triangular_solve_vector; + +template +struct triangular_solve_matrix; + +// small helper struct extracting some traits on the underlying solver operation +template +class trsolve_traits +{ + private: + enum { + RhsIsVectorAtCompileTime = (Side==OnTheLeft ? Rhs::ColsAtCompileTime : Rhs::RowsAtCompileTime)==1 + }; + public: + enum { + Unrolling = (RhsIsVectorAtCompileTime && Rhs::SizeAtCompileTime != Dynamic && Rhs::SizeAtCompileTime <= 8) + ? CompleteUnrolling : NoUnrolling, + RhsVectors = RhsIsVectorAtCompileTime ? 1 : Dynamic + }; +}; + +template::Unrolling, + int RhsVectors = trsolve_traits::RhsVectors + > +struct triangular_solver_selector; + +template +struct triangular_solver_selector +{ + typedef typename Lhs::Scalar LhsScalar; + typedef typename Rhs::Scalar RhsScalar; + typedef blas_traits LhsProductTraits; + typedef typename LhsProductTraits::ExtractType ActualLhsType; + typedef Map, Aligned> MappedRhs; + static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs) + { + ActualLhsType actualLhs = LhsProductTraits::extract(lhs); + + // FIXME find a way to allow an inner stride if packet_traits::size==1 + + bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1; + + ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhs,rhs.size(), + (useRhsDirectly ? rhs.data() : 0)); + + if(!useRhsDirectly) + MappedRhs(actualRhs,rhs.size()) = rhs; + + triangular_solve_vector + ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs); + + if(!useRhsDirectly) + rhs = MappedRhs(actualRhs, rhs.size()); + } +}; + +// the rhs is a matrix +template +struct triangular_solver_selector +{ + typedef typename Rhs::Scalar Scalar; + typedef blas_traits LhsProductTraits; + typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType; + + static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs) + { + typename internal::add_const_on_value_type::type actualLhs = LhsProductTraits::extract(lhs); + + const Index size = lhs.rows(); + const Index othersize = Side==OnTheLeft? rhs.cols() : rhs.rows(); + + typedef internal::gemm_blocking_space<(Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor,Scalar,Scalar, + Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType; + + BlockingType blocking(rhs.rows(), rhs.cols(), size, 1, false); + + triangular_solve_matrix + ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.innerStride(), rhs.outerStride(), blocking); + } +}; + +/*************************************************************************** +* meta-unrolling implementation +***************************************************************************/ + +template +struct triangular_solver_unroller; + +template +struct triangular_solver_unroller { + enum { + IsLower = ((Mode&Lower)==Lower), + DiagIndex = IsLower ? LoopIndex : Size - LoopIndex - 1, + StartIndex = IsLower ? 0 : DiagIndex+1 + }; + static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs) + { + if (LoopIndex>0) + rhs.coeffRef(DiagIndex) -= lhs.row(DiagIndex).template segment(StartIndex).transpose() + .cwiseProduct(rhs.template segment(StartIndex)).sum(); + + if(!(Mode & UnitDiag)) + rhs.coeffRef(DiagIndex) /= lhs.coeff(DiagIndex,DiagIndex); + + triangular_solver_unroller::run(lhs,rhs); + } +}; + +template +struct triangular_solver_unroller { + static EIGEN_DEVICE_FUNC void run(const Lhs&, Rhs&) {} +}; + +template +struct triangular_solver_selector { + static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs) + { triangular_solver_unroller::run(lhs,rhs); } +}; + +template +struct triangular_solver_selector { + static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs) + { + Transpose trLhs(lhs); + Transpose trRhs(rhs); + + triangular_solver_unroller,Transpose, + ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), + 0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs); + } +}; + +} // end namespace internal + +/*************************************************************************** +* TriangularView methods +***************************************************************************/ + +#ifndef EIGEN_PARSED_BY_DOXYGEN +template +template +EIGEN_DEVICE_FUNC void TriangularViewImpl::solveInPlace(const MatrixBase& _other) const +{ + OtherDerived& other = _other.const_cast_derived(); + eigen_assert( derived().cols() == derived().rows() && ((Side==OnTheLeft && derived().cols() == other.rows()) || (Side==OnTheRight && derived().cols() == other.cols())) ); + eigen_assert((!(int(Mode) & int(ZeroDiag))) && bool(int(Mode) & (int(Upper) | int(Lower)))); + // If solving for a 0x0 matrix, nothing to do, simply return. + if (derived().cols() == 0) + return; + + enum { copy = (internal::traits::Flags & RowMajorBit) && OtherDerived::IsVectorAtCompileTime && OtherDerived::SizeAtCompileTime!=1}; + typedef typename internal::conditional::type, OtherDerived&>::type OtherCopy; + OtherCopy otherCopy(other); + + internal::triangular_solver_selector::type, + Side, Mode>::run(derived().nestedExpression(), otherCopy); + + if (copy) + other = otherCopy; +} + +template +template +const internal::triangular_solve_retval,Other> +TriangularViewImpl::solve(const MatrixBase& other) const +{ + return internal::triangular_solve_retval(derived(), other.derived()); +} +#endif + +namespace internal { + + +template +struct traits > +{ + typedef typename internal::plain_matrix_type_column_major::type ReturnType; +}; + +template struct triangular_solve_retval + : public ReturnByValue > +{ + typedef typename remove_all::type RhsNestedCleaned; + typedef ReturnByValue Base; + + triangular_solve_retval(const TriangularType& tri, const Rhs& rhs) + : m_triangularMatrix(tri), m_rhs(rhs) + {} + + inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_rhs.rows(); } + inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_rhs.cols(); } + + template inline void evalTo(Dest& dst) const + { + if(!is_same_dense(dst,m_rhs)) + dst = m_rhs; + m_triangularMatrix.template solveInPlace(dst); + } + + protected: + const TriangularType& m_triangularMatrix; + typename Rhs::Nested m_rhs; +}; + +} // namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SOLVETRIANGULAR_H -- cgit v1.2.3-70-g09d2