mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add LU::transpose().solve() and LU::adjoint().solve() API.

Browse Source
This commit is contained in:
Gael Guennebaud 2015-12-01 14:38:47 +01:00
parent 1663d15da7
commit 0bb12fa614
23 changed files with 267 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Eigen/Core
View File

@ -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"

4
Eigen/src/CholmodSupport/CholmodSupport.h
View File

@ -170,6 +170,10 @@ class CholmodBase : public SparseSolverBase<Derived>
typedef typename MatrixType::RealScalar RealScalar;
typedef MatrixType CholMatrixType;
typedef typename MatrixType::StorageIndex StorageIndex;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:

1
Eigen/src/Core/CoreEvaluators.h
View File

@ -29,6 +29,7 @@ struct storage_kind_to_evaluator_kind {
template<typename StorageKind> struct storage_kind_to_shape;
template<> struct storage_kind_to_shape<Dense> { typedef DenseShape Shape; };
template<> struct storage_kind_to_shape<SolverStorage> { typedef SolverShape Shape; };
template<> struct storage_kind_to_shape<PermutationStorage> { typedef PermutationShape Shape; };
template<> struct storage_kind_to_shape<TranspositionsStorage> { typedef TranspositionsShape Shape; };

24
Eigen/src/Core/Inverse.h
View File

@ -48,6 +48,7 @@ public:
typedef typename internal::ref_selector<XprType>::type XprTypeNested;
typedef typename internal::remove_all<XprTypeNested>::type XprTypeNestedCleaned;
typedef typename internal::ref_selector<Inverse>::type Nested;
typedef typename internal::remove_all<XprType>::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<typename XprType>
class InverseImpl<XprType,Dense>
: public MatrixBase<Inverse<XprType> >
// Generic API dispatcher
template<typename XprType, typename StorageKind>
class InverseImpl
: public internal::generic_xpr_base<Inverse<XprType> >::type
{
typedef Inverse<XprType> Derived;
public:
typedef MatrixBase<Derived> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
typedef typename internal::remove_all<XprType>::type NestedExpression;
typedef typename internal::generic_xpr_base<Inverse<XprType> >::type Base;
typedef typename XprType::Scalar Scalar;
private:
Scalar coeff(Index row, Index col) const;
Scalar coeff(Index i) const;
};

27
Eigen/src/Core/Solve.h
View File

@ -34,12 +34,11 @@ template<typename Decomposition, typename RhsType,typename StorageKind> struct s
template<typename Decomposition, typename RhsType>
struct solve_traits<Decomposition,RhsType,Dense>
{
typedef typename Decomposition::MatrixType MatrixType;
typedef Matrix<typename RhsType::Scalar,
MatrixType::ColsAtCompileTime,
Decomposition::ColsAtCompileTime,
RhsType::ColsAtCompileTime,
RhsType::PlainObject::Options,
MatrixType::MaxColsAtCompileTime,
Decomposition::MaxColsAtCompileTime,
RhsType::MaxColsAtCompileTime> PlainObject;
};
@ -145,6 +144,28 @@ struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar
}
};
// Specialization for "dst = dec.transpose().solve(rhs)"
template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
struct Assignment<DstXprType, Solve<Transpose<const DecType>,RhsType>, internal::assign_op<Scalar>, Dense2Dense, Scalar>
{
typedef Solve<Transpose<const DecType>,RhsType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
src.dec().nestedExpression().template _solve_impl_transposed<false>(src.rhs(), dst);
}
};
// Specialization for "dst = dec.adjoint().solve(rhs)"
template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
struct Assignment<DstXprType, Solve<CwiseUnaryOp<internal::scalar_conjugate_op<typename DecType::Scalar>, const Transpose<const DecType> >,RhsType>, internal::assign_op<Scalar>, Dense2Dense, Scalar>
{
typedef Solve<CwiseUnaryOp<internal::scalar_conjugate_op<typename DecType::Scalar>, const Transpose<const DecType> >,RhsType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
src.dec().nestedExpression().nestedExpression().template _solve_impl_transposed<true>(src.rhs(), dst);
}
};
} // end namepsace internal
} // end namespace Eigen

130
Eigen/src/Core/SolverBase.h Normal file
View File

@ -0,0 +1,130 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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<typename Derived>
class SolverBase : public EigenBase<Derived>
{
public:
typedef EigenBase<Derived> Base;
typedef typename internal::traits<Derived>::Scalar Scalar;
typedef Scalar CoeffReturnType;
enum {
RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
internal::traits<Derived>::ColsAtCompileTime>::ret),
MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
MaxSizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::MaxRowsAtCompileTime,
internal::traits<Derived>::MaxColsAtCompileTime>::ret),
IsVectorAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime == 1
|| internal::traits<Derived>::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<typename Rhs>
inline const Solve<Derived, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
return Solve<Derived, Rhs>(derived(), b.derived());
}
/** \internal the return type of transpose() */
typedef typename internal::add_const<Transpose<const Derived> >::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<NumTraits<Scalar>::IsComplex,
CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, 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<typename Derived>
struct generic_xpr_base<Derived, MatrixXpr, SolverStorage>
{
typedef SolverBase<Derived> type;
};
} // end namespace internal
} // end namespace Eigen
#endif // EIGEN_SOLVERBASE_H

2
Eigen/src/Core/Transpose.h
View File

@ -39,7 +39,7 @@ struct traits<Transpose<MatrixType> > : public traits<MatrixType>
MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
Flags0 = MatrixTypeNestedPlain::Flags & ~(LvalueBit | NestByRefBit),
Flags0 = traits<MatrixTypeNestedPlain>::Flags & ~(LvalueBit | NestByRefBit),
Flags1 = Flags0 | FlagsLvalueBit,
Flags = Flags1 ^ RowMajorBit,
InnerStrideAtCompileTime = inner_stride_at_compile_time<MatrixType>::ret,

4
Eigen/src/Core/util/Constants.h
View File

@ -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"; } };

1
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -132,6 +132,7 @@ template<typename MatrixType> struct CommaInitializer;
template<typename Derived> class ReturnByValue;
template<typename ExpressionType> class ArrayWrapper;
template<typename ExpressionType> class MatrixWrapper;
template<typename Derived> class SolverBase;
template<typename XprType> class InnerIterator;
namespace internal {

6
Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
View File

@ -39,8 +39,10 @@ class DiagonalPreconditioner
typedef Matrix<Scalar,Dynamic,1> 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<Scalar,Dynamic,Dynamic> MatrixType;
enum {
ColsAtCompileTime = Dynamic,
MaxColsAtCompileTime = Dynamic
};
DiagonalPreconditioner() : m_isInitialized(false) {}

5
Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h
View File

@ -57,12 +57,15 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
typedef typename OrderingType::PermutationType PermutationType;
typedef typename PermutationType::StorageIndex StorageIndex;
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
typedef FactorType MatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorSx;
typedef Matrix<RealScalar,Dynamic,1> VectorRx;
typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
typedef std::vector<std::list<StorageIndex> > VectorList;
enum { UpLo = _UpLo };
enum {
ColsAtCompileTime = Dynamic,
MaxColsAtCompileTime = Dynamic
};
public:
/** Default constructor leaving the object in a partly non-initialized stage.

8
Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
View File

@ -109,11 +109,13 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageInd
typedef Matrix<StorageIndex,Dynamic,1> VectorI;
typedef SparseMatrix<Scalar,RowMajor,StorageIndex> 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<Scalar,Dynamic,Dynamic> MatrixType;
IncompleteLUT()
: m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10),
m_analysisIsOk(false), m_factorizationIsOk(false)

5
Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
View File

@ -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;

16
Eigen/src/LU/FullPivLU.h
View File

@ -16,6 +16,8 @@ namespace internal {
template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> >
: traits<_MatrixType>
{
typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind;
enum { Flags = 0 };
};
@ -53,21 +55,18 @@ template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> >
* \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
*/
template<typename _MatrixType> class FullPivLU
: public SolverBase<FullPivLU<_MatrixType> >
{
public:
typedef _MatrixType MatrixType;
typedef SolverBase<FullPivLU> 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<typename MatrixType::Scalar>::Real RealScalar;
typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
// FIXME should be int
typedef typename MatrixType::StorageIndex StorageIndex;
typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
@ -223,6 +222,7 @@ template<typename _MatrixType> class FullPivLU
*
* \sa TriangularView::solve(), kernel(), inverse()
*/
// FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
template<typename Rhs>
inline const Solve<FullPivLU, Rhs>
solve(const MatrixBase<Rhs>& b) const

27
Eigen/src/LU/PartialPivLU.h
View File

@ -17,6 +17,8 @@ namespace internal {
template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> >
: traits<_MatrixType>
{
typedef MatrixXpr XprKind;
typedef SolverStorage StorageKind;
typedef traits<_MatrixType> BaseTraits;
enum {
Flags = BaseTraits::Flags & RowMajorBit,
@ -58,33 +60,29 @@ template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> >
* \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
*/
template<typename _MatrixType> class PartialPivLU
: public SolverBase<PartialPivLU<_MatrixType> >
{
public:
typedef _MatrixType MatrixType;
typedef SolverBase<PartialPivLU> 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<typename MatrixType::Scalar>::Real RealScalar;
typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
// FIXME should be int
typedef typename MatrixType::StorageIndex StorageIndex;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> 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<typename _MatrixType> class PartialPivLU
*
* \sa TriangularView::solve(), inverse(), computeInverse()
*/
// FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
template<typename Rhs>
inline const Solve<PartialPivLU, Rhs>
solve(const MatrixBase<Rhs>& b) const
@ -508,7 +507,7 @@ MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
return res;
}
/***** Implementation of solve() *****************************************************/
/***** Implementation details *****************************************************/
namespace internal {

4
Eigen/src/PaStiXSupport/PaStiXSupport.h
View File

@ -141,6 +141,10 @@ class PastixBase : public SparseSolverBase<Derived>
typedef typename MatrixType::StorageIndex StorageIndex;
typedef Matrix<Scalar,Dynamic,1> Vector;
typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:

4
Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
View File

@ -68,6 +68,10 @@ class SPQR : public SparseSolverBase<SPQR<_MatrixType> >
typedef SuiteSparse_long StorageIndex ;
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType;
typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType;
enum {
ColsAtCompileTime = Dynamic,
MaxColsAtCompileTime = Dynamic
};
public:
SPQR()
: m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)

5
Eigen/src/SparseCholesky/SimplicialCholesky.h
View File

@ -71,6 +71,11 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived>
typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef Matrix<StorageIndex,Dynamic,1> VectorI;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:
using Base::derived;

5
Eigen/src/SparseLU/SparseLU.h
View File

@ -90,6 +90,11 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
typedef internal::SparseLUImpl<Scalar, StorageIndex> 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)

6
Eigen/src/SparseQR/SparseQR.h
View File

@ -84,6 +84,12 @@ class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
typedef Matrix<StorageIndex, Dynamic, 1> IndexVector;
typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:
SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
{ }

4
Eigen/src/SuperLUSupport/SuperLUSupport.h
View File

@ -304,6 +304,10 @@ class SuperLUBase : public SparseSolverBase<Derived>
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Map<PermutationMatrix<Dynamic,Dynamic,int> > PermutationMap;
typedef SparseMatrix<Scalar> LUMatrixType;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:

4
Eigen/src/UmfPackSupport/UmfPackSupport.h
View File

@ -146,6 +146,10 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
typedef SparseMatrix<Scalar> LUMatrixType;
typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType;
typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
public:

22
test/lu.cpp
View File

@ -99,6 +99,9 @@ template<typename MatrixType> void lu_non_invertible()
m3 = MatrixType::Random(rows,cols2);
lu.template _solve_impl_transposed<false>(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<typename MatrixType> void lu_non_invertible()
m3 = MatrixType::Random(rows,cols2);
lu.template _solve_impl_transposed<true>(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<typename MatrixType> void lu_invertible()
@ -138,12 +144,20 @@ template<typename MatrixType> 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<false>(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<true>(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<typename MatrixType> 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<false>(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<true>(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<typename MatrixType> void lu_verify_assert()