// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008-2011 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/.
// for compatibility with super old version of umfpack, // not sure this is really needed, but this is harmless. #ifndef SuiteSparse_long #ifdef UF_long #define SuiteSparse_long UF_long #else #error neither SuiteSparse_long nor UF_long are defined #endif #endif
namespace Eigen {
/* TODO extract L, extract U, compute det, etc... */
// Get Lunz inlineint umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
{ return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
}
inlineint umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
{ return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
}
// Get Numeric inlineint umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
{ return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
}
/** \ingroup UmfPackSupport_Module * \brief A sparse LU factorization and solver based on UmfPack * * This class allows to solve for A.X = B sparse linear problems via a LU factorization * using the UmfPack library. The sparse matrix A must be squared and full rank. * The vectors or matrices X and B can be either dense or sparse. * * \warning The input matrix A should be in a \b compressed and \b column-major form. * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix. * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * * \implsparsesolverconcept * * \sa \ref TutorialSparseSolverConcept, class SparseLU
*/ template<typename _MatrixType> class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
{ protected: typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base; using Base::m_isInitialized; public: using Base::_solve_impl; typedef _MatrixType MatrixType; typedeftypename MatrixType::Scalar Scalar; typedeftypename MatrixType::RealScalar RealScalar; typedeftypename MatrixType::StorageIndex StorageIndex; typedef Matrix<Scalar,Dynamic,1> Vector; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; typedef SparseMatrix<Scalar> LUMatrixType; typedef SparseMatrix<Scalar,ColMajor,StorageIndex> UmfpackMatrixType; typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef; enum {
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
inline Index rows() const { return mp_matrix.rows(); } inline Index cols() const { return mp_matrix.cols(); }
/** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was successful, * \c NumericalIssue if the matrix.appears to be negative.
*/
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info;
}
/** Computes the sparse Cholesky decomposition of \a matrix * Note that the matrix should be column-major, and in compressed format for best performance. * \sa SparseMatrix::makeCompressed().
*/ template<typename InputMatrixType> void compute(const InputMatrixType& matrix)
{ if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex()); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
grab(matrix.derived());
analyzePattern_impl();
factorize_impl();
}
/** Performs a symbolic decomposition on the sparcity of \a matrix. * * This function is particularly useful when solving for several problems having the same structure. * * \sa factorize(), compute()
*/ template<typename InputMatrixType> void analyzePattern(const InputMatrixType& matrix)
{ if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex()); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
grab(matrix.derived());
analyzePattern_impl();
}
/** Provides the return status code returned by UmfPack during the numeric * factorization. * * \sa factorize(), compute()
*/ inlineint umfpackFactorizeReturncode() const
{
eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()"); return m_fact_errorCode;
}
/** Provides access to the control settings array used by UmfPack. * * If this array contains NaN's, the default values are used. * * See UMFPACK documentation for details.
*/ inlineconst UmfpackControl& umfpackControl() const
{ return m_control;
}
/** Provides access to the control settings array used by UmfPack. * * If this array contains NaN's, the default values are used. * * See UMFPACK documentation for details.
*/ inline UmfpackControl& umfpackControl()
{ return m_control;
}
/** Performs a numeric decomposition of \a matrix * * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed. * * \sa analyzePattern(), compute()
*/ template<typename InputMatrixType> void factorize(const InputMatrixType& matrix)
{
eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); if(m_numeric)
umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
grab(matrix.derived());
factorize_impl();
}
/** Prints the current UmfPack control settings. * * \sa umfpackControl()
*/ void printUmfpackControl()
{
umfpack_report_control(m_control.data(), Scalar(),StorageIndex());
}
/** Prints statistics collected by UmfPack. * * \sa analyzePattern(), compute()
*/ void printUmfpackInfo()
{
eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(),StorageIndex());
}
/** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization). * * \sa analyzePattern(), compute()
*/ void printUmfpackStatus() {
eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(),StorageIndex());
}
template<typename MatrixType> template<typename BDerived,typename XDerived> bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
{
Index 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");
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
