aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/Core1
-rw-r--r--Eigen/src/Core/SolveTriangular.h4
-rw-r--r--Eigen/src/Core/TriangularMatrix.h15
-rw-r--r--Eigen/src/Core/products/SelfadjointMatrixVector.h8
-rw-r--r--Eigen/src/Core/products/TriangularMatrixVector.h161
-rw-r--r--Eigen/src/Core/util/Macros.h8
-rw-r--r--test/CMakeLists.txt1
-rw-r--r--test/product_triangular.cpp92
-rw-r--r--test/triangular.cpp4
9 files changed, 281 insertions, 13 deletions
diff --git a/Eigen/Core b/Eigen/Core
index 18e6a6045..8627ada53 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -181,6 +181,7 @@ namespace Eigen {
#include "src/Core/SelfAdjointView.h"
#include "src/Core/SolveTriangular.h"
#include "src/Core/products/SelfadjointRank2Update.h"
+#include "src/Core/products/TriangularMatrixVector.h"
} // namespace Eigen
diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h
index 200b4a325..76e4289de 100644
--- a/Eigen/src/Core/SolveTriangular.h
+++ b/Eigen/src/Core/SolveTriangular.h
@@ -45,7 +45,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,RowMajor>
};
static void run(const Lhs& lhs, Rhs& other)
{
- static const int PanelWidth = EIGEN_TUNE_TRSV_PANEL_WIDTH;
+ static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
const ActualLhsType& actualLhs = LhsProductTraits::extract(lhs);
const int size = lhs.cols();
@@ -104,7 +104,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,ColMajor>
static void run(const Lhs& lhs, Rhs& other)
{
- static const int PanelWidth = EIGEN_TUNE_TRSV_PANEL_WIDTH;
+ static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
const ActualLhsType& actualLhs = LhsProductTraits::extract(lhs);
const int size = lhs.cols();
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 81621dbcc..506be7965 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -142,13 +142,17 @@ struct ei_traits<TriangularView<MatrixType, _Mode> > : ei_traits<MatrixType>
};
};
-template<typename MatrixType, unsigned int _Mode> class TriangularView
- : public TriangularBase<TriangularView<MatrixType, _Mode> >
+template<typename Lhs,typename Rhs>
+struct ei_triangular_vector_product_returntype;
+
+template<typename _MatrixType, unsigned int _Mode> class TriangularView
+ : public TriangularBase<TriangularView<_MatrixType, _Mode> >
{
public:
typedef TriangularBase<TriangularView> Base;
typedef typename ei_traits<TriangularView>::Scalar Scalar;
+ typedef _MatrixType MatrixType;
typedef typename MatrixType::PlainMatrixType PlainMatrixType;
enum {
@@ -244,6 +248,13 @@ template<typename MatrixType, unsigned int _Mode> class TriangularView
}
template<typename OtherDerived>
+ ei_triangular_vector_product_returntype<TriangularView,OtherDerived>
+ operator*(const MatrixBase<OtherDerived>& rhs) const
+ {
+ return ei_triangular_vector_product_returntype<TriangularView,OtherDerived>(*this, rhs.derived(), 1);
+ }
+
+ template<typename OtherDerived>
typename ei_plain_matrix_type_column_major<OtherDerived>::type
solve(const MatrixBase<OtherDerived>& other) const;
diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h
index aa3187a07..8ac83772c 100644
--- a/Eigen/src/Core/products/SelfadjointMatrixVector.h
+++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h
@@ -25,7 +25,7 @@
#ifndef EIGEN_SELFADJOINT_MATRIX_VECTOR_H
#define EIGEN_SELFADJOINT_MATRIX_VECTOR_H
-/* Optimized col-major selfadjoint matrix * vector product:
+/* Optimized selfadjoint matrix * vector product:
* This algorithm processes 2 columns at onces that allows to both reduce
* the number of load/stores of the result by a factor 2 and to reduce
* the instruction dependency.
@@ -78,11 +78,11 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector(
size_t alignedStart = (starti) + ei_alignmentOffset(&res[starti], endi-starti);
alignedEnd = alignedStart + ((endi-alignedStart)/(PacketSize))*(PacketSize);
- res[j] += t0 * conj0(A0[j]);
+ res[j] += t0 * conj0(A0[j]);
if(FirstTriangular)
{
- res[j+1] += t1 * conj0(A1[j+1]);
- res[j] += t1 * conj0(A1[j]);
+ res[j+1] += t1 * conj0(A1[j+1]);
+ res[j] += t1 * conj0(A1[j]);
t3 += conj1(A1[j]) * rhs[j];
}
else
diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h
new file mode 100644
index 000000000..533aad170
--- /dev/null
+++ b/Eigen/src/Core/products/TriangularMatrixVector.h
@@ -0,0 +1,161 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRIANGULARMATRIXVECTOR_H
+#define EIGEN_TRIANGULARMATRIXVECTOR_H
+
+template<typename MatrixType, typename Rhs, typename Result,
+ int Mode, bool ConjLhs, bool ConjRhs, int StorageOrder>
+struct ei_product_triangular_vector_selector;
+
+template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs>
+struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,ColMajor>
+{
+ typedef typename Rhs::Scalar Scalar;
+ enum {
+ IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit),
+ HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit
+ };
+ static void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha)
+ {
+ static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
+ typename ei_conj_expr_if<ConjLhs,Lhs>::ret cjLhs(lhs);
+ typename ei_conj_expr_if<ConjRhs,Rhs>::ret cjRhs(rhs);
+
+ int size = lhs.cols();
+ for (int pi=0; pi<size; pi+=PanelWidth)
+ {
+ int actualPanelWidth = std::min(PanelWidth, size-pi);
+ for (int k=0; k<actualPanelWidth; ++k)
+ {
+ int i = pi + k;
+ int s = IsLowerTriangular ? (HasUnitDiag ? i+1 : i ) : pi;
+ int r = IsLowerTriangular ? actualPanelWidth-k : k+1;
+ if ((!HasUnitDiag) || (--r)>0)
+ res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r);
+ if (HasUnitDiag)
+ res.coeffRef(i) += alpha * cjRhs.coeff(i);
+ }
+ int r = IsLowerTriangular ? size - pi - actualPanelWidth : pi;
+ if (r>0)
+ {
+ int s = IsLowerTriangular ? pi+actualPanelWidth : 0;
+ ei_cache_friendly_product_colmajor_times_vector<ConjLhs,ConjRhs>(
+ r,
+ &(lhs.const_cast_derived().coeffRef(s,pi)), lhs.stride(),
+ rhs.segment(pi, actualPanelWidth),
+ &(res.coeffRef(s)),
+ alpha);
+ }
+ }
+ }
+};
+
+template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs>
+struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,RowMajor>
+{
+ typedef typename Rhs::Scalar Scalar;
+ enum {
+ IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit),
+ HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit
+ };
+ static void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha)
+ {
+ static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
+ typename ei_conj_expr_if<ConjLhs,Lhs>::ret cjLhs(lhs);
+ typename ei_conj_expr_if<ConjRhs,Rhs>::ret cjRhs(rhs);
+ int size = lhs.cols();
+ for (int pi=0; pi<size; pi+=PanelWidth)
+ {
+ int actualPanelWidth = std::min(PanelWidth, size-pi);
+ for (int k=0; k<actualPanelWidth; ++k)
+ {
+ int i = pi + k;
+ int s = IsLowerTriangular ? pi : (HasUnitDiag ? i+1 : i);
+ int r = IsLowerTriangular ? k+1 : actualPanelWidth-k;
+ if ((!HasUnitDiag) || (--r)>0)
+ res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwise() * cjRhs.segment(s,r).transpose()).sum();
+ if (HasUnitDiag)
+ res.coeffRef(i) += alpha * cjRhs.coeff(i);
+ }
+ int r = IsLowerTriangular ? pi : size - pi - actualPanelWidth;
+ if (r>0)
+ {
+ int s = IsLowerTriangular ? 0 : pi + actualPanelWidth;
+ Block<Result,Dynamic,1> target(res,pi,0,actualPanelWidth,1);
+ ei_cache_friendly_product_rowmajor_times_vector<ConjLhs,ConjRhs>(
+ &(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(),
+ &(rhs.const_cast_derived().coeffRef(s)), r,
+ target, alpha);
+ }
+ }
+ }
+};
+
+template<typename Lhs,typename Rhs>
+struct ei_triangular_vector_product_returntype
+ : public ReturnByValue<ei_triangular_vector_product_returntype<Lhs,Rhs>,
+ Matrix<typename ei_traits<Rhs>::Scalar,
+ Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{
+ typedef typename Lhs::Scalar Scalar;
+ typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
+ ei_triangular_vector_product_returntype(const Lhs& lhs, const Rhs& rhs, Scalar alpha)
+ : m_lhs(lhs), m_rhs(rhs), m_alpha(alpha)
+ {}
+
+ template<typename Dest> void evalTo(Dest& dst) const
+ {
+ typedef typename Lhs::MatrixType MatrixType;
+
+ typedef ei_blas_traits<MatrixType> LhsBlasTraits;
+ typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
+ typedef typename ei_cleantype<ActualLhsType>::type _ActualLhsType;
+ const ActualLhsType actualLhs = LhsBlasTraits::extract(m_lhs._expression());
+
+ typedef ei_blas_traits<Rhs> RhsBlasTraits;
+ typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
+ typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
+ const ActualRhsType actualRhs = RhsBlasTraits::extract(m_rhs);
+
+ Scalar actualAlpha = m_alpha * LhsBlasTraits::extractScalarFactor(m_lhs._expression())
+ * RhsBlasTraits::extractScalarFactor(m_rhs);
+
+ dst.resize(m_rhs.rows(), m_rhs.cols());
+ dst.setZero();
+ ei_product_triangular_vector_selector
+ <_ActualLhsType,_ActualRhsType,Dest,
+ ei_traits<Lhs>::Mode,
+ LhsBlasTraits::NeedToConjugate,
+ RhsBlasTraits::NeedToConjugate,
+ ei_traits<Lhs>::Flags&RowMajorBit>
+ ::run(actualLhs,actualRhs,dst,actualAlpha);
+ }
+
+ const Lhs m_lhs;
+ const typename Rhs::Nested m_rhs;
+ const Scalar m_alpha;
+};
+
+#endif // EIGEN_TRIANGULARMATRIXVECTOR_H
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index bb73b03cc..0cc5ae9aa 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -94,11 +94,11 @@
#define EIGEN_TUNE_FOR_CPU_CACHE_SIZE (sizeof(float)*256*256)
#endif
-/** Defines the maximal width of the blocks used in the triangular solver
- * for vectors (level 2 blas xTRSV). The default is 8.
+/** Defines the maximal width of the blocks used in the triangular product and solver
+ * for vectors (level 2 blas xTRMV and xTRSV). The default is 8.
*/
-#ifndef EIGEN_TUNE_TRSV_PANEL_WIDTH
-#define EIGEN_TUNE_TRSV_PANEL_WIDTH 8
+#ifndef EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH
+#define EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH 8
#endif
/** Allows to disable some optimizations which might affect the accuracy of the result.
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 229c7123b..886731d45 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -110,6 +110,7 @@ ei_add_test(array)
ei_add_test(array_replicate)
ei_add_test(array_reverse)
ei_add_test(triangular)
+ei_add_test(product_triangular)
ei_add_test(cholesky " " "${GSL_LIBRARIES}")
ei_add_test(lu ${EI_OFLAG})
ei_add_test(determinant)
diff --git a/test/product_triangular.cpp b/test/product_triangular.cpp
new file mode 100644
index 000000000..876fb4388
--- /dev/null
+++ b/test/product_triangular.cpp
@@ -0,0 +1,92 @@
+// This file is triangularView of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+template<typename MatrixType> void product_triangular(const MatrixType& m)
+{
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+
+ RealScalar largerEps = 10*test_precision<RealScalar>();
+
+ int rows = m.rows();
+ int cols = m.cols();
+
+ MatrixType m1 = MatrixType::Random(rows, cols),
+ m3(rows, cols);
+ VectorType v1 = VectorType::Random(rows);
+
+ Scalar s1 = ei_random<Scalar>();
+
+ m1 = MatrixType::Random(rows, cols);
+
+ // check with a column-major matrix
+ m3 = m1.template triangularView<Eigen::LowerTriangular>();
+ VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::UpperTriangular>();
+ VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UpperTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::UnitLowerTriangular>();
+ VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
+ VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps));
+
+ // check conjugated and scalar multiple expressions (col-major)
+ m3 = m1.template triangularView<Eigen::LowerTriangular>();
+ VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::UpperTriangular>();
+ VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::UpperTriangular>() * v1.conjugate(), largerEps));
+
+ // check with a row-major matrix
+ m3 = m1.template triangularView<Eigen::UpperTriangular>();
+ VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::LowerTriangular>();
+ VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UpperTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
+ VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::UnitLowerTriangular>();
+ VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps));
+
+ // check conjugated and scalar multiple expressions (row-major)
+ m3 = m1.template triangularView<Eigen::UpperTriangular>();
+ VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
+ m3 = m1.template triangularView<Eigen::LowerTriangular>();
+ VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::UpperTriangular>() * (s1*v1.conjugate()), largerEps));
+ m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
+
+ // TODO check with sub-matrices
+}
+
+void test_product_triangular()
+{
+ for(int i = 0; i < g_repeat ; i++) {
+ CALL_SUBTEST( product_triangular(Matrix<float, 1, 1>()) );
+ CALL_SUBTEST( product_triangular(Matrix<float, 2, 2>()) );
+ CALL_SUBTEST( product_triangular(Matrix3d()) );
+ CALL_SUBTEST( product_triangular(Matrix<std::complex<float>,23, 23>()) );
+ CALL_SUBTEST( product_triangular(MatrixXcd(17,17)) );
+ CALL_SUBTEST( product_triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(19, 19)) );
+ }
+}
diff --git a/test/triangular.cpp b/test/triangular.cpp
index 7c680a8ed..6385bffd1 100644
--- a/test/triangular.cpp
+++ b/test/triangular.cpp
@@ -51,6 +51,8 @@ template<typename MatrixType> void triangular(const MatrixType& m)
v2 = VectorType::Random(rows),
vzero = VectorType::Zero(rows);
+ Scalar s1 = ei_random<Scalar>();
+
MatrixType m1up = m1.template triangularView<Eigen::UpperTriangular>();
MatrixType m2up = m2.template triangularView<Eigen::UpperTriangular>();
@@ -142,7 +144,7 @@ void test_triangular()
CALL_SUBTEST( triangular(Matrix3d()) );
CALL_SUBTEST( triangular(MatrixXcf(4, 4)) );
CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) );
- CALL_SUBTEST( triangular(MatrixXd(17,17)) );
+ CALL_SUBTEST( triangular(MatrixXcd(17,17)) );
CALL_SUBTEST( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) );
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-05 02:47:29 +0000