diff --git a/Eigen/Core b/Eigen/Core index 7cf431320..1ec749452 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -391,6 +391,7 @@ using std::ptrdiff_t; #include "src/Core/GeneralProduct.h" #include "src/Core/Solve.h" #include "src/Core/Inverse.h" +#include "src/Core/SolverBase.h" #include "src/Core/PermutationMatrix.h" #include "src/Core/Transpositions.h" #include "src/Core/TriangularMatrix.h" diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index f33aa9bf1..06421d5ed 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -170,6 +170,10 @@ class CholmodBase : public SparseSolverBase typedef typename MatrixType::RealScalar RealScalar; typedef MatrixType CholMatrixType; typedef typename MatrixType::StorageIndex StorageIndex; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; public: diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index a8b359085..42ad452f7 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -29,6 +29,7 @@ struct storage_kind_to_evaluator_kind { template struct storage_kind_to_shape; template<> struct storage_kind_to_shape { typedef DenseShape Shape; }; +template<> struct storage_kind_to_shape { typedef SolverShape Shape; }; template<> struct storage_kind_to_shape { typedef PermutationShape Shape; }; template<> struct storage_kind_to_shape { typedef TranspositionsShape Shape; }; diff --git a/Eigen/src/Core/Inverse.h b/Eigen/src/Core/Inverse.h index 8ba1a12d9..f3ec84990 100644 --- a/Eigen/src/Core/Inverse.h +++ b/Eigen/src/Core/Inverse.h @@ -48,6 +48,7 @@ public: typedef typename internal::ref_selector::type XprTypeNested; typedef typename internal::remove_all::type XprTypeNestedCleaned; typedef typename internal::ref_selector::type Nested; + typedef typename internal::remove_all::type NestedExpression; explicit Inverse(const XprType &xpr) : m_xpr(xpr) @@ -62,25 +63,16 @@ protected: XprTypeNested m_xpr; }; -/** \internal - * Specialization of the Inverse expression for dense expressions. - * Direct access to the coefficients are discared. - * FIXME this intermediate class is probably not needed anymore. - */ -template -class InverseImpl - : public MatrixBase > +// Generic API dispatcher +template +class InverseImpl + : public internal::generic_xpr_base >::type { - typedef Inverse Derived; - public: - - typedef MatrixBase Base; - EIGEN_DENSE_PUBLIC_INTERFACE(Derived) - typedef typename internal::remove_all::type NestedExpression; - + typedef typename internal::generic_xpr_base >::type Base; + typedef typename XprType::Scalar Scalar; private: - + Scalar coeff(Index row, Index col) const; Scalar coeff(Index i) const; }; diff --git a/Eigen/src/Core/Solve.h b/Eigen/src/Core/Solve.h index 2d163fe2a..ba2ee53b8 100644 --- a/Eigen/src/Core/Solve.h +++ b/Eigen/src/Core/Solve.h @@ -34,12 +34,11 @@ template struct s template struct solve_traits { - typedef typename Decomposition::MatrixType MatrixType; typedef Matrix PlainObject; }; @@ -145,6 +144,28 @@ struct Assignment, internal::assign_op +struct Assignment,RhsType>, internal::assign_op, Dense2Dense, Scalar> +{ + typedef Solve,RhsType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + { + src.dec().nestedExpression().template _solve_impl_transposed(src.rhs(), dst); + } +}; + +// Specialization for "dst = dec.adjoint().solve(rhs)" +template +struct Assignment, const Transpose >,RhsType>, internal::assign_op, Dense2Dense, Scalar> +{ + typedef Solve, const Transpose >,RhsType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + { + src.dec().nestedExpression().nestedExpression().template _solve_impl_transposed(src.rhs(), dst); + } +}; + } // end namepsace internal } // end namespace Eigen diff --git a/Eigen/src/Core/SolverBase.h b/Eigen/src/Core/SolverBase.h new file mode 100644 index 000000000..8a4adc229 --- /dev/null +++ b/Eigen/src/Core/SolverBase.h @@ -0,0 +1,130 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2015 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_SOLVERBASE_H +#define EIGEN_SOLVERBASE_H + +namespace Eigen { + +namespace internal { + + + +} // end namespace internal + +/** \class SolverBase + * \brief A base class for matrix decomposition and solvers + * + * \tparam Derived the actual type of the decomposition/solver. + * + * Any matrix decomposition inheriting this base class provide the following API: + * + * \code + * MatrixType A, b, x; + * DecompositionType dec(A); + * x = dec.solve(b); // solve A * x = b + * x = dec.transpose().solve(b); // solve A^T * x = b + * x = dec.adjoint().solve(b); // solve A' * x = b + * \endcode + * + * \warning Currently, any other usage of transpose() and adjoint() are not supported and will produce compilation errors. + * + * \sa class PartialPivLU, class FullPivLU + */ +template +class SolverBase : public EigenBase +{ + public: + + typedef EigenBase Base; + typedef typename internal::traits::Scalar Scalar; + typedef Scalar CoeffReturnType; + + enum { + RowsAtCompileTime = internal::traits::RowsAtCompileTime, + ColsAtCompileTime = internal::traits::ColsAtCompileTime, + SizeAtCompileTime = (internal::size_at_compile_time::RowsAtCompileTime, + internal::traits::ColsAtCompileTime>::ret), + MaxRowsAtCompileTime = internal::traits::MaxRowsAtCompileTime, + MaxColsAtCompileTime = internal::traits::MaxColsAtCompileTime, + MaxSizeAtCompileTime = (internal::size_at_compile_time::MaxRowsAtCompileTime, + internal::traits::MaxColsAtCompileTime>::ret), + IsVectorAtCompileTime = internal::traits::MaxRowsAtCompileTime == 1 + || internal::traits::MaxColsAtCompileTime == 1 + }; + + /** Default constructor */ + SolverBase() + {} + + ~SolverBase() + {} + + using Base::derived; + + /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A. + */ + template + inline const Solve + solve(const MatrixBase& b) const + { + eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); + return Solve(derived(), b.derived()); + } + + /** \internal the return type of transpose() */ + typedef typename internal::add_const >::type ConstTransposeReturnType; + /** \returns an expression of the transposed of the factored matrix. + * + * A typical usage is to solve for the transposed problem A^T x = b: + * \code x = dec.transpose().solve(b); \endcode + * + * \sa adjoint(), solve() + */ + inline ConstTransposeReturnType transpose() const + { + return ConstTransposeReturnType(derived()); + } + + /** \internal the return type of adjoint() */ + typedef typename internal::conditional::IsComplex, + CwiseUnaryOp, ConstTransposeReturnType>, + ConstTransposeReturnType + >::type AdjointReturnType; + /** \returns an expression of the adjoint of the factored matrix + * + * A typical usage is to solve for the adjoint problem A' x = b: + * \code x = dec.adjoint().solve(b); \endcode + * + * For real scalar types, this function is equivalent to transpose(). + * + * \sa transpose(), solve() + */ + inline AdjointReturnType adjoint() const + { + return AdjointReturnType(derived().transpose()); + } + + protected: +}; + +namespace internal { + +template +struct generic_xpr_base +{ + typedef SolverBase type; + +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SOLVERBASE_H diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 2152405d5..5b66eb5e1 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -39,7 +39,7 @@ struct traits > : public traits MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime, FlagsLvalueBit = is_lvalue::value ? LvalueBit : 0, - Flags0 = MatrixTypeNestedPlain::Flags & ~(LvalueBit | NestByRefBit), + Flags0 = traits::Flags & ~(LvalueBit | NestByRefBit), Flags1 = Flags0 | FlagsLvalueBit, Flags = Flags1 ^ RowMajorBit, InnerStrideAtCompileTime = inner_stride_at_compile_time::ret, diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 28852c8c3..a364f48d1 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -492,6 +492,9 @@ struct Dense {}; /** The type used to identify a general sparse storage. */ struct Sparse {}; +/** The type used to identify a general solver (foctored) storage. */ +struct SolverStorage {}; + /** The type used to identify a permutation storage. */ struct PermutationStorage {}; @@ -506,6 +509,7 @@ struct ArrayXpr {}; // An evaluator must define its shape. By default, it can be one of the following: struct DenseShape { static std::string debugName() { return "DenseShape"; } }; +struct SolverShape { static std::string debugName() { return "SolverShape"; } }; struct HomogeneousShape { static std::string debugName() { return "HomogeneousShape"; } }; struct DiagonalShape { static std::string debugName() { return "DiagonalShape"; } }; struct BandShape { static std::string debugName() { return "BandShape"; } }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 1aa81abf8..483af876f 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -132,6 +132,7 @@ template struct CommaInitializer; template class ReturnByValue; template class ArrayWrapper; template class MatrixWrapper; +template class SolverBase; template class InnerIterator; namespace internal { diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h index b850630a3..358444aff 100644 --- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h +++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h @@ -39,8 +39,10 @@ class DiagonalPreconditioner typedef Matrix Vector; public: typedef typename Vector::StorageIndex StorageIndex; - // this typedef is only to export the scalar type and compile-time dimensions to solve_retval - typedef Matrix MatrixType; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; DiagonalPreconditioner() : m_isInitialized(false) {} diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h index 8f549af82..284e37f13 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h @@ -57,12 +57,15 @@ class IncompleteCholesky : public SparseSolverBase FactorType; - typedef FactorType MatrixType; typedef Matrix VectorSx; typedef Matrix VectorRx; typedef Matrix VectorIx; typedef std::vector > VectorList; enum { UpLo = _UpLo }; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; public: /** Default constructor leaving the object in a partly non-initialized stage. diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 519472377..338e6f10a 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -109,11 +109,13 @@ class IncompleteLUT : public SparseSolverBase VectorI; typedef SparseMatrix FactorType; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; + public: - // this typedef is only to export the scalar type and compile-time dimensions to solve_retval - typedef Matrix MatrixType; - IncompleteLUT() : m_droptol(NumTraits::dummy_precision()), m_fillfactor(10), m_analysisIsOk(false), m_factorizationIsOk(false) diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 5f4bcea11..e51ff7280 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -31,6 +31,11 @@ public: typedef typename MatrixType::StorageIndex StorageIndex; typedef typename MatrixType::RealScalar RealScalar; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + public: using Base::derived; diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 4691efd2f..0c4d63923 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -16,6 +16,8 @@ namespace internal { template struct traits > : traits<_MatrixType> { + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; enum { Flags = 0 }; }; @@ -53,21 +55,18 @@ template struct traits > * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse() */ template class FullPivLU + : public SolverBase > { public: typedef _MatrixType MatrixType; + typedef SolverBase Base; + + EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU) + // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef typename internal::traits::StorageKind StorageKind; - // FIXME should be int - typedef typename MatrixType::StorageIndex StorageIndex; typedef typename internal::plain_row_type::type IntRowVectorType; typedef typename internal::plain_col_type::type IntColVectorType; typedef PermutationMatrix PermutationQType; @@ -223,6 +222,7 @@ template class FullPivLU * * \sa TriangularView::solve(), kernel(), inverse() */ + // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion. template inline const Solve solve(const MatrixBase& b) const diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 91abbc341..50e920609 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -17,6 +17,8 @@ namespace internal { template struct traits > : traits<_MatrixType> { + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; typedef traits<_MatrixType> BaseTraits; enum { Flags = BaseTraits::Flags & RowMajorBit, @@ -58,33 +60,29 @@ template struct traits > * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU */ template class PartialPivLU + : public SolverBase > { public: typedef _MatrixType MatrixType; + typedef SolverBase Base; + EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU) + // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef typename internal::traits::StorageKind StorageKind; - // FIXME should be int - typedef typename MatrixType::StorageIndex StorageIndex; typedef PermutationMatrix PermutationType; typedef Transpositions TranspositionType; typedef typename MatrixType::PlainObject PlainObject; /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via PartialPivLU::compute(const MatrixType&). - */ + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via PartialPivLU::compute(const MatrixType&). + */ PartialPivLU(); /** \brief Default Constructor with memory preallocation @@ -145,6 +143,7 @@ template class PartialPivLU * * \sa TriangularView::solve(), inverse(), computeInverse() */ + // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion. template inline const Solve solve(const MatrixBase& b) const @@ -508,7 +507,7 @@ MatrixType PartialPivLU::reconstructedMatrix() const return res; } -/***** Implementation of solve() *****************************************************/ +/***** Implementation details *****************************************************/ namespace internal { diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index cec4149e7..1999fd289 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -141,6 +141,10 @@ class PastixBase : public SparseSolverBase typedef typename MatrixType::StorageIndex StorageIndex; typedef Matrix Vector; typedef SparseMatrix ColSpMatrix; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; public: diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h index d2f053fa5..d9c3113e7 100644 --- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h @@ -68,6 +68,10 @@ class SPQR : public SparseSolverBase > typedef SuiteSparse_long StorageIndex ; typedef SparseMatrix MatrixType; typedef Map > PermutationType; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; public: SPQR() : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits::epsilon()), m_useDefaultThreshold(true) diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index ef612cf45..1343eb15c 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -71,6 +71,11 @@ class SimplicialCholeskyBase : public SparseSolverBase typedef Matrix VectorType; typedef Matrix VectorI; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + public: using Base::derived; diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 73368cba4..acd3ad100 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -90,6 +90,11 @@ class SparseLU : public SparseSolverBase >, typedef Matrix IndexVector; typedef PermutationMatrix PermutationType; typedef internal::SparseLUImpl Base; + + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; public: SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index bbd337c40..4f26c19ca 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -84,6 +84,12 @@ class SparseQR : public SparseSolverBase > typedef Matrix IndexVector; typedef Matrix ScalarVector; typedef PermutationMatrix PermutationType; + + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + public: SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { } diff --git a/Eigen/src/SuperLUSupport/SuperLUSupport.h b/Eigen/src/SuperLUSupport/SuperLUSupport.h index c145e25bd..b20da37f7 100644 --- a/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -304,6 +304,10 @@ class SuperLUBase : public SparseSolverBase typedef Matrix IntColVectorType; typedef Map > PermutationMap; typedef SparseMatrix LUMatrixType; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; public: diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index caac082f3..aaec8c6f1 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -146,6 +146,10 @@ class UmfPackLU : public SparseSolverBase > typedef SparseMatrix LUMatrixType; typedef SparseMatrix UmfpackMatrixType; typedef Ref UmfpackMatrixRef; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; public: diff --git a/test/lu.cpp b/test/lu.cpp index f753fc74a..f14435114 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -99,6 +99,9 @@ template void lu_non_invertible() m3 = MatrixType::Random(rows,cols2); lu.template _solve_impl_transposed(m2, m3); VERIFY_IS_APPROX(m2, m1.transpose()*m3); + m3 = MatrixType::Random(rows,cols2); + m3 = lu.transpose().solve(m2); + VERIFY_IS_APPROX(m2, m1.transpose()*m3); // test solve with conjugate transposed m3 = MatrixType::Random(rows,cols2); @@ -106,6 +109,9 @@ template void lu_non_invertible() m3 = MatrixType::Random(rows,cols2); lu.template _solve_impl_transposed(m2, m3); VERIFY_IS_APPROX(m2, m1.adjoint()*m3); + m3 = MatrixType::Random(rows,cols2); + m3 = lu.adjoint().solve(m2); + VERIFY_IS_APPROX(m2, m1.adjoint()*m3); } template void lu_invertible() @@ -138,12 +144,20 @@ template void lu_invertible() m2 = lu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); VERIFY_IS_APPROX(m2, lu.inverse()*m3); + // test solve with transposed lu.template _solve_impl_transposed(m3, m2); VERIFY_IS_APPROX(m3, m1.transpose()*m2); + m3 = MatrixType::Random(size,size); + m3 = lu.transpose().solve(m2); + VERIFY_IS_APPROX(m2, m1.transpose()*m3); + // test solve with conjugate transposed lu.template _solve_impl_transposed(m3, m2); VERIFY_IS_APPROX(m3, m1.adjoint()*m2); + m3 = MatrixType::Random(size,size); + m3 = lu.adjoint().solve(m2); + VERIFY_IS_APPROX(m2, m1.adjoint()*m3); // Regression test for Bug 302 MatrixType m4 = MatrixType::Random(size,size); @@ -168,12 +182,20 @@ template void lu_partial_piv() m2 = plu.solve(m3); VERIFY_IS_APPROX(m3, m1*m2); VERIFY_IS_APPROX(m2, plu.inverse()*m3); + // test solve with transposed plu.template _solve_impl_transposed(m3, m2); VERIFY_IS_APPROX(m3, m1.transpose()*m2); + m3 = MatrixType::Random(size,size); + m3 = plu.transpose().solve(m2); + VERIFY_IS_APPROX(m2, m1.transpose()*m3); + // test solve with conjugate transposed plu.template _solve_impl_transposed(m3, m2); VERIFY_IS_APPROX(m3, m1.adjoint()*m2); + m3 = MatrixType::Random(size,size); + m3 = plu.adjoint().solve(m2); + VERIFY_IS_APPROX(m2, m1.adjoint()*m3); } template void lu_verify_assert()