devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merged local branch with main trunk

Browse Source
This commit is contained in:
Desire NUENTSA
2013-01-25 19:05:33 +01:00
parent e4ec63aee7 81d4bfa8d9
commit d58056bde4
8 changed files with 50 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -157,27 +157,6 @@ class PastixBase : internal::noncopyable
template<typename Rhs,typename Dest>
bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
/** \internal */
template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
{
eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
eigen_assert(rows()==b.rows());
// we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
static const int NbColsAtOnce = 1;
int rhsCols = b.cols();
int size = b.rows();
Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
for(int k=0; k<rhsCols; k+=NbColsAtOnce)
{
int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
}
}
Derived& derived()
{
return *static_cast<Derived*>(this);
@@ -731,7 +710,7 @@ struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec()._solve_sparse(rhs(),dst);
this->defaultEvalTo(dst);
}
};

25
Eigen/src/PardisoSupport/PardisoSupport.h
View File

@@ -206,29 +206,6 @@ class PardisoImpl
template<typename BDerived, typename XDerived>
bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
/** \internal */
template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
{
eigen_assert(m_size==b.rows());
// we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
static const int NbColsAtOnce = 4;
int rhsCols = b.cols();
int size = b.rows();
// Pardiso cannot solve in-place,
// so we need two temporaries
Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_rhs(size,rhsCols);
Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_res(size,rhsCols);
for(int k=0; k<rhsCols; k+=NbColsAtOnce)
{
int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
}
}
protected:
void pardisoRelease()
{
@@ -604,7 +581,7 @@ struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec().derived()._solve_sparse(rhs(),dst);
this->defaultEvalTo(dst);
}
};

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

@@ -215,27 +215,6 @@ class SimplicialCholeskyBase : internal::noncopyable
dest = m_Pinv * dest;
}
/** \internal */
template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
{
eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
eigen_assert(m_matrix.rows()==b.rows());
// we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
static const int NbColsAtOnce = 4;
int rhsCols = b.cols();
int size = b.rows();
Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
for(int k=0; k<rhsCols; k+=NbColsAtOnce)
{
int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
}
}
#endif // EIGEN_PARSED_BY_DOXYGEN
protected:
@@ -864,7 +843,7 @@ struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec().derived()._solve_sparse(rhs(),dst);
this->defaultEvalTo(dst);
}
};

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

@@ -353,14 +353,14 @@ class SuperLUBase : internal::noncopyable
*
* \sa compute()
*/
// template<typename Rhs>
// inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
// {
// eigen_assert(m_isInitialized && "SuperLU is not initialized.");
// eigen_assert(rows()==b.rows()
// && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
// return internal::sparse_solve_retval<SuperLU, Rhs>(*this, b.derived());
// }
template<typename Rhs>
inline const internal::sparse_solve_retval<SuperLUBase, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "SuperLU is not initialized.");
eigen_assert(rows()==b.rows()
&& "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
return internal::sparse_solve_retval<SuperLUBase, Rhs>(*this, b.derived());
}
/** Performs a symbolic decomposition on the sparcity of \a matrix.
*
@@ -1015,7 +1015,7 @@ struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec().derived()._solve(rhs(),dst);
this->defaultEvalTo(dst);
}
};

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

@@ -215,14 +215,14 @@ class UmfPackLU : internal::noncopyable
*
* \sa compute()
*/
// template<typename Rhs>
// inline const internal::sparse_solve_retval<UmfPAckLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
// {
// eigen_assert(m_isInitialized && "UmfPAckLU is not initialized.");
// eigen_assert(rows()==b.rows()
// && "UmfPAckLU::solve(): invalid number of rows of the right hand side matrix b");
// return internal::sparse_solve_retval<UmfPAckLU, Rhs>(*this, b.derived());
// }
template<typename Rhs>
inline const internal::sparse_solve_retval<UmfPackLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
eigen_assert(rows()==b.rows()
&& "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
return internal::sparse_solve_retval<UmfPackLU, Rhs>(*this, b.derived());
}
/** Performs a symbolic decomposition on the sparcity of \a matrix.
*
@@ -381,7 +381,8 @@ bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDe
const int rhsCols = b.cols();
eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
int errorCode;
for (int j=0; j<rhsCols; ++j)
{
@@ -420,7 +421,7 @@ struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec()._solve(rhs(),dst);
this->defaultEvalTo(dst);
}
};

17
Eigen/src/misc/SparseSolve.h
View File

@@ -47,6 +47,23 @@ template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_b
}
protected:
template<typename DestScalar, int DestOptions, typename DestIndex>
inline void defaultEvalTo(SparseMatrix<DestScalar,DestOptions,DestIndex>& dst) const
{
// we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
static const int NbColsAtOnce = 4;
int rhsCols = m_rhs.cols();
int size = m_rhs.rows();
Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,rhsCols);
for(int k=0; k<rhsCols; k+=NbColsAtOnce)
{
int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
tmp.leftCols(actualCols) = m_rhs.middleCols(k,actualCols);
tmpX.leftCols(actualCols) = m_dec.solve(tmp.leftCols(actualCols));
dst.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView();
}
}
const DecompositionType& m_dec;
typename Rhs::Nested m_rhs;
};