GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-11 19:29:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

* Bye bye MultiplierBase, extend a bit AnyMatrixBase to allow =, +=, and -=

Browse Source 
* This probably makes ReturnByValue needless
This commit is contained in:
Gael Guennebaud 2009-08-03 16:05:15 +02:00
parent ce1dc1ab16
commit 3cf5bb31f6
15 changed files with 126 additions and 113 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -57,7 +57,7 @@ struct ei_traits<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
};
template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
class BandMatrix : public MultiplierBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
{
public:

17
Eigen/src/Core/DiagonalMatrix.h
View File

@ -27,7 +27,7 @@
#define EIGEN_DIAGONALMATRIX_H
template<typename Derived>
class DiagonalBase : public MultiplierBase<Derived>
class DiagonalBase : public AnyMatrixBase<Derived>
{
public:
typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType;
@ -52,6 +52,12 @@ class DiagonalBase : public MultiplierBase<Derived>
DenseMatrixType toDenseMatrix() const { return derived(); }
template<typename DenseDerived>
void evalToDense(MatrixBase<DenseDerived> &other) const;
template<typename DenseDerived>
void addToDense(MatrixBase<DenseDerived> &other) const
{ other.diagonal() += diagonal(); }
template<typename DenseDerived>
void subToDense(MatrixBase<DenseDerived> &other) const
{ other.diagonal() -= diagonal(); }
inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
@ -84,6 +90,7 @@ void DiagonalBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const
*/
template<typename _Scalar, int _Size>
struct ei_traits<DiagonalMatrix<_Scalar,_Size> >
: ei_traits<Matrix<_Scalar,_Size,_Size> >
{
typedef Matrix<_Scalar,_Size,1> DiagonalVectorType;
};
@ -170,6 +177,14 @@ template<typename _DiagonalVectorType>
struct ei_traits<DiagonalWrapper<_DiagonalVectorType> >
{
typedef _DiagonalVectorType DiagonalVectorType;
typedef typename DiagonalVectorType::Scalar Scalar;
enum {
RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
MaxRowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
MaxColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
Flags = 0
};
};
template<typename _DiagonalVectorType>

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

@ -462,6 +462,7 @@ class Matrix
: m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols())
{
_check_template_params();
resize(other.rows(), other.cols());
*this = other;
}

52
Eigen/src/Core/MatrixBase.h
View File

@ -37,19 +37,31 @@
*/
template<typename Derived> struct AnyMatrixBase
{
typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
Derived& derived() { return *static_cast<Derived*>(this); }
const Derived& derived() const { return *static_cast<const Derived*>(this); }
};
/** Common base class for all classes T such that there are overloaded operator* allowing to
* multiply a MatrixBase by a T on both sides.
*
* In other words, an AnyMatrixBase object is an object that can be multiplied a MatrixBase, the result being again a MatrixBase.
*
* Besides MatrixBase-derived classes, this also includes certain special matrix classes, such as diagonal matrices.
*/
template<typename Derived> struct MultiplierBase : public AnyMatrixBase<Derived>
{
using AnyMatrixBase<Derived>::derived;
/** \returns the number of rows. \sa cols(), RowsAtCompileTime */
inline int rows() const { return derived().rows(); }
/** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
inline int cols() const { return derived().cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{ derived().evalTo(dst); }
template<typename Dest> inline void addToDense(Dest& dst) const
{
typename Dest::PlainMatrixType res(rows(),cols());
evalToDense(res);
dst += res;
}
template<typename Dest> inline void subToDense(Dest& dst) const
{
typename Dest::PlainMatrixType res(rows(),cols());
evalToDense(res);
dst -= res;
}
};
/** \class MatrixBase
@ -79,7 +91,7 @@ template<typename Derived> struct MultiplierBase : public AnyMatrixBase<Derived>
*/
template<typename Derived> class MatrixBase
#ifndef EIGEN_PARSED_BY_DOXYGEN
: public MultiplierBase<Derived>,
: public AnyMatrixBase<Derived>,
public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
#endif // not EIGEN_PARSED_BY_DOXYGEN
@ -298,12 +310,16 @@ template<typename Derived> class MatrixBase
Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
{ other.derived().evalToDense(derived()); return derived(); }
template<typename OtherDerived>
Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
{ other.derived().addToDense(derived()); return derived(); }
template<typename OtherDerived>
Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
{ other.derived().subToDense(derived()); return derived(); }
template<typename OtherDerived,typename OtherEvalType>
Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
template<typename OtherDerived,typename OtherEvalType>
Derived& operator+=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
template<typename OtherDerived,typename OtherEvalType>
Derived& operator-=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** Copies \a other into *this without evaluating other. \returns a reference to *this. */
@ -410,7 +426,7 @@ template<typename Derived> class MatrixBase
operator*(const MatrixBase<OtherDerived> &other) const;
template<typename OtherDerived>
Derived& operator*=(const MultiplierBase<OtherDerived>& other);
Derived& operator*=(const AnyMatrixBase<OtherDerived>& other);
template<typename DiagonalDerived>
const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
@ -645,7 +661,7 @@ template<typename Derived> class MatrixBase
void visit(Visitor& func) const;
#ifndef EIGEN_PARSED_BY_DOXYGEN
using MultiplierBase<Derived>::derived;
using AnyMatrixBase<Derived>::derived;
inline Derived& const_cast_derived() const
{ return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
#endif // not EIGEN_PARSED_BY_DOXYGEN

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

@ -294,7 +294,7 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
template<typename Derived>
template<typename OtherDerived>
inline Derived &
MatrixBase<Derived>::operator*=(const MultiplierBase<OtherDerived> &other)
MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
{
return derived() = derived() * other.derived();
}

32
Eigen/src/Core/ReturnByValue.h
View File

@ -48,14 +48,6 @@ template<typename Functor, typename EvalType> class ReturnByValue
public:
template<typename Dest> inline void evalTo(Dest& dst) const
{ static_cast<const Functor*>(this)->evalTo(dst); }
template<typename Dest> inline void addTo(Dest& dst) const
{ static_cast<const Functor*>(this)->_addTo(dst); }
template<typename Dest> inline void subTo(Dest& dst) const
{ static_cast<const Functor*>(this)->_subTo(dst); }
template<typename Dest> inline void _addTo(Dest& dst) const
{ EvalType res; evalTo(res); dst += res; }
template<typename Dest> inline void _subTo(Dest& dst) const
{ EvalType res; evalTo(res); dst -= res; }
};
template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols>
@ -68,14 +60,6 @@ template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int
template<typename Dest>
inline void evalTo(Dest& dst) const
{ static_cast<const Functor* const>(this)->evalTo(dst); }
template<typename Dest> inline void addTo(Dest& dst) const
{ static_cast<const Functor*>(this)->_addTo(dst); }
template<typename Dest> inline void subTo(Dest& dst) const
{ static_cast<const Functor*>(this)->_subTo(dst); }
template<typename Dest> inline void _addTo(Dest& dst) const
{ EvalType res; evalTo(res); dst += res; }
template<typename Dest> inline void _subTo(Dest& dst) const
{ EvalType res; evalTo(res); dst -= res; }
inline int rows() const { return static_cast<const Functor* const>(this)->rows(); }
inline int cols() const { return static_cast<const Functor* const>(this)->cols(); }
};
@ -88,20 +72,4 @@ Derived& MatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived,OtherEv
return derived();
}
template<typename Derived>
template<typename OtherDerived,typename OtherEvalType>
Derived& MatrixBase<Derived>::operator+=(const ReturnByValue<OtherDerived,OtherEvalType>& other)
{
other.addTo(derived());
return derived();
}
template<typename Derived>
template<typename OtherDerived,typename OtherEvalType>
Derived& MatrixBase<Derived>::operator-=(const ReturnByValue<OtherDerived,OtherEvalType>& other)
{
other.subTo(derived());
return derived();
}
#endif // EIGEN_RETURNBYVALUE_H

48
Eigen/src/Core/SelfAdjointView.h
View File

@ -55,7 +55,7 @@ struct ei_traits<SelfAdjointView<MatrixType, TriangularPart> > : ei_traits<Matri
template <typename Lhs, int LhsMode, bool LhsIsVector,
typename Rhs, int RhsMode, bool RhsIsVector>
struct ei_selfadjoint_product_returntype;
struct SelfadjointProductMatrix;
// FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
@ -100,20 +100,20 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
/** Efficient self-adjoint matrix times vector/matrix product */
template<typename OtherDerived>
ei_selfadjoint_product_returntype<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
operator*(const MatrixBase<OtherDerived>& rhs) const
{
return ei_selfadjoint_product_returntype
return SelfadjointProductMatrix
<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
(m_matrix, rhs.derived());
}
/** Efficient vector/matrix times self-adjoint matrix product */
template<typename OtherDerived> friend
ei_selfadjoint_product_returntype<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
{
return ei_selfadjoint_product_returntype
return SelfadjointProductMatrix
<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
(lhs.derived(),rhs.m_matrix);
}
@ -201,10 +201,13 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Dynami
***************************************************************************/
template<typename Lhs, int LhsMode, typename Rhs>
struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>
: public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>,
Matrix<typename ei_traits<Rhs>::Scalar,
Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> >
: ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
{};
template<typename Lhs, int LhsMode, typename Rhs>
struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
: public AnyMatrixBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> >
{
typedef typename Lhs::Scalar Scalar;
@ -224,19 +227,19 @@ struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>
LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit)
};
ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
{}
inline int rows() const { return m_lhs.rows(); }
inline int cols() const { return m_rhs.cols(); }
template<typename Dest> inline void _addTo(Dest& dst) const
template<typename Dest> inline void addToDense(Dest& dst) const
{ evalTo(dst,1); }
template<typename Dest> inline void _subTo(Dest& dst) const
template<typename Dest> inline void subToDense(Dest& dst) const
{ evalTo(dst,-1); }
template<typename Dest> void evalTo(Dest& dst) const
template<typename Dest> void evalToDense(Dest& dst) const
{
dst.setZero();
evalTo(dst,1);
@ -272,12 +275,15 @@ struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>
***************************************************************************/
template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>
: public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>,
Matrix<typename ei_traits<Rhs>::Scalar,
Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
: ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
{};
template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>
: public AnyMatrixBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
{
ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
{}
@ -305,12 +311,12 @@ struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>
RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit
};
template<typename Dest> inline void _addTo(Dest& dst) const
template<typename Dest> inline void addToDense(Dest& dst) const
{ evalTo(dst,1); }
template<typename Dest> inline void _subTo(Dest& dst) const
template<typename Dest> inline void subToDense(Dest& dst) const
{ evalTo(dst,-1); }
template<typename Dest> void evalTo(Dest& dst) const
template<typename Dest> void evalToDense(Dest& dst) const
{
dst.setZero();
evalTo(dst,1);

12
Eigen/src/Core/TriangularMatrix.h
View File

@ -43,7 +43,7 @@
*
* \sa MatrixBase::part()
*/
template<typename Derived> class TriangularBase : public MultiplierBase<Derived>
template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
{
public:
@ -145,7 +145,7 @@ struct ei_traits<TriangularView<MatrixType, _Mode> > : ei_traits<MatrixType>
template<int Mode, bool LhsIsTriangular,
typename Lhs, bool LhsIsVector,
typename Rhs, bool RhsIsVector>
struct ei_triangular_product_returntype;
struct TriangularProduct;
template<typename _MatrixType, unsigned int _Mode> class TriangularView
: public TriangularBase<TriangularView<_MatrixType, _Mode> >
@ -253,20 +253,20 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
/** Efficient triangular matrix times vector/matrix product */
template<typename OtherDerived>
ei_triangular_product_returntype<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
TriangularProduct<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
operator*(const MatrixBase<OtherDerived>& rhs) const
{
return ei_triangular_product_returntype
return TriangularProduct
<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
(m_matrix, rhs.derived());
}
/** Efficient vector/matrix times triangular matrix product */
template<typename OtherDerived> friend
ei_triangular_product_returntype<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
TriangularProduct<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
{
return ei_triangular_product_returntype
return TriangularProduct
<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
(lhs.derived(),rhs.m_matrix);
}

19
Eigen/src/Core/products/TriangularMatrixMatrix.h
View File

@ -321,12 +321,15 @@ struct ei_product_triangular_matrix_matrix<Scalar,Mode,false,
***************************************************************************/
template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
struct ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
: public ReturnByValue<ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false>,
Matrix<typename ei_traits<Rhs>::Scalar,
Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
struct ei_traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
: ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
{};
template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
struct TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
: public AnyMatrixBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
{
ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs)
TriangularProduct(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
{}
@ -347,12 +350,12 @@ struct ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false
typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
template<typename Dest> inline void _addTo(Dest& dst) const
template<typename Dest> inline void addToDense(Dest& dst) const
{ evalTo(dst,1); }
template<typename Dest> inline void _subTo(Dest& dst) const
template<typename Dest> inline void subToDense(Dest& dst) const
{ evalTo(dst,-1); }
template<typename Dest> void evalTo(Dest& dst) const
template<typename Dest> void evalToDense(Dest& dst) const
{
dst.resize(m_lhs.rows(), m_rhs.cols());
dst.setZero();

19
Eigen/src/Core/products/TriangularMatrixVector.h
View File

@ -118,10 +118,13 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs
***************************************************************************/
template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs>
struct ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>
: public ReturnByValue<ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>,
Matrix<typename ei_traits<Rhs>::Scalar,
Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
struct ei_traits<TriangularProduct<Mode,true,Lhs,false,Rhs,true> >
: ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
{};
template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs>
struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
: public AnyMatrixBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true> >
{
typedef typename Lhs::Scalar Scalar;
@ -137,19 +140,19 @@ struct ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>
typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs)
TriangularProduct(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
{}
inline int rows() const { return m_lhs.rows(); }
inline int cols() const { return m_rhs.cols(); }
template<typename Dest> inline void _addTo(Dest& dst) const
template<typename Dest> inline void addToDense(Dest& dst) const
{ evalTo(dst,1); }
template<typename Dest> inline void _subTo(Dest& dst) const
template<typename Dest> inline void subToDense(Dest& dst) const
{ evalTo(dst,-1); }
template<typename Dest> void evalTo(Dest& dst) const
template<typename Dest> void evalToDense(Dest& dst) const
{
dst.setZero();
evalTo(dst,1);

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

@ -29,7 +29,6 @@ template<typename T> struct ei_traits;
template<typename T> struct NumTraits;
template<typename Derived> struct AnyMatrixBase;
template<typename Derived> struct MultiplierBase;
template<typename _Scalar, int _Rows, int _Cols,
int _Options = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION | AutoAlign,

18
bench/btl/libs/C_BLAS/C_BLAS_interface.hh
View File

@ -237,6 +237,15 @@ public :
spotrf_(&uplo, &N, C, &N, &info);
}
static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
int N2 = N*N;
scopy_(&N2, X, &intone, C, &intone);
char uplo = 'L';
int info = 0;
int * ipiv = (int*)alloca(sizeof(int)*N);
sgetrf_(&N, &N, C, &N, ipiv, &info);
}
#ifdef HAS_LAPACK
static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
@ -249,14 +258,7 @@ public :
sgetc2_(&N, C, &N, ipiv, jpiv, &info);
}
static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
int N2 = N*N;
scopy_(&N2, X, &intone, C, &intone);
char uplo = 'L';
int info = 0;
int * ipiv = (int*)alloca(sizeof(int)*N);
sgetrf_(&N, &N, C, &N, ipiv, &info);
}
static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
#ifdef PUREBLAS

4
bench/btl/libs/C_BLAS/main.cpp
View File

@ -24,7 +24,7 @@
#include "action_cholesky.hh"
#include "action_lu_decomp.hh"
// #include "action_partial_lu_decomp.hh"
#include "action_partial_lu.hh"
#include "action_trisolve_matrix.hh"
#ifdef HAS_LAPACK
@ -52,10 +52,10 @@ int main()
bench<Action_trisolve_matrix<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_cholesky<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_partial_lu<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
#ifdef HAS_LAPACK
bench<Action_lu_decomp<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
// bench<Action_partial_lu_decomp<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_hessenberg<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_tridiagonalization<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
#endif

8
bench/btl/libs/eigen2/eigen2_interface.hh
View File

@ -190,10 +190,10 @@ public :
}
static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
// X = L.template triangularView<LowerTriangular>().solve(B);
X = B;
ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular>
::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride());
X = L.template triangularView<LowerTriangular>().solve(B);
//
// ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular>
// ::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride());
}
static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){

4
bench/btl/libs/eigen2/main_adv.cpp
View File

@ -23,7 +23,7 @@
#include "action_cholesky.hh"
#include "action_hessenberg.hh"
#include "action_lu_decomp.hh"
// #include "action_partial_lu_decomp.hh"
#include "action_partial_lu.hh"
BTL_MAIN;
@ -33,7 +33,7 @@ int main()
bench<Action_trisolve_matrix<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_cholesky<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_lu_decomp<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
// bench<Action_partial_lu_decomp<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_partial_lu<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_hessenberg<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_tridiagonalization<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);