aboutsummaryrefslogtreecommitdiffhomepage
path: root/third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h
diff options
context:
space:
mode:
authorGravatar Benoit Steiner <benoit.steiner.goog@gmail.com>2016-01-12 11:11:40 -0800
committerGravatar Benoit Steiner <benoit.steiner.goog@gmail.com>2016-01-12 11:11:40 -0800
commit42673c3a588ce4cc20b02ab20e5e9d38b64a3cb4 (patch)
tree1545b8e2411a774728685a4da519058897d49ee5 /third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h
parentccf01f9d77b28b649777f5a937a295f6dee2a130 (diff)
Deleted the remainder of the local copy of eigen that is shipped with
TensorFlow.
Diffstat (limited to 'third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h')
-rw-r--r--third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h1026
1 files changed, 0 insertions, 1026 deletions
diff --git a/third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h b/third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h
deleted file mode 100644
index bcb355760c..0000000000
--- a/third_party/eigen3/Eigen/src/SuperLUSupport/SuperLUSupport.h
+++ /dev/null
@@ -1,1026 +0,0 @@
-// 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/.
-
-#ifndef EIGEN_SUPERLUSUPPORT_H
-#define EIGEN_SUPERLUSUPPORT_H
-
-namespace Eigen {
-
-#define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \
- extern "C" { \
- typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t; \
- extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
- char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
- void *, int, SuperMatrix *, SuperMatrix *, \
- FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \
- PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \
- } \
- inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \
- int *perm_c, int *perm_r, int *etree, char *equed, \
- FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
- SuperMatrix *U, void *work, int lwork, \
- SuperMatrix *B, SuperMatrix *X, \
- FLOATTYPE *recip_pivot_growth, \
- FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
- SuperLUStat_t *stats, int *info, KEYTYPE) { \
- PREFIX##mem_usage_t mem_usage; \
- PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
- U, work, lwork, B, X, recip_pivot_growth, rcond, \
- ferr, berr, &mem_usage, stats, info); \
- return mem_usage.for_lu; /* bytes used by the factor storage */ \
- }
-
-DECL_GSSVX(s,float,float)
-DECL_GSSVX(c,float,std::complex<float>)
-DECL_GSSVX(d,double,double)
-DECL_GSSVX(z,double,std::complex<double>)
-
-#ifdef MILU_ALPHA
-#define EIGEN_SUPERLU_HAS_ILU
-#endif
-
-#ifdef EIGEN_SUPERLU_HAS_ILU
-
-// similarly for the incomplete factorization using gsisx
-#define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \
- extern "C" { \
- extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \
- char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \
- void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \
- PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \
- } \
- inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
- int *perm_c, int *perm_r, int *etree, char *equed, \
- FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
- SuperMatrix *U, void *work, int lwork, \
- SuperMatrix *B, SuperMatrix *X, \
- FLOATTYPE *recip_pivot_growth, \
- FLOATTYPE *rcond, \
- SuperLUStat_t *stats, int *info, KEYTYPE) { \
- PREFIX##mem_usage_t mem_usage; \
- PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \
- U, work, lwork, B, X, recip_pivot_growth, rcond, \
- &mem_usage, stats, info); \
- return mem_usage.for_lu; /* bytes used by the factor storage */ \
- }
-
-DECL_GSISX(s,float,float)
-DECL_GSISX(c,float,std::complex<float>)
-DECL_GSISX(d,double,double)
-DECL_GSISX(z,double,std::complex<double>)
-
-#endif
-
-template<typename MatrixType>
-struct SluMatrixMapHelper;
-
-/** \internal
- *
- * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
- * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
- *
- * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
- */
-struct SluMatrix : SuperMatrix
-{
- SluMatrix()
- {
- Store = &storage;
- }
-
- SluMatrix(const SluMatrix& other)
- : SuperMatrix(other)
- {
- Store = &storage;
- storage = other.storage;
- }
-
- SluMatrix& operator=(const SluMatrix& other)
- {
- SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
- Store = &storage;
- storage = other.storage;
- return *this;
- }
-
- struct
- {
- union {int nnz;int lda;};
- void *values;
- int *innerInd;
- int *outerInd;
- } storage;
-
- void setStorageType(Stype_t t)
- {
- Stype = t;
- if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
- Store = &storage;
- else
- {
- eigen_assert(false && "storage type not supported");
- Store = 0;
- }
- }
-
- template<typename Scalar>
- void setScalarType()
- {
- if (internal::is_same<Scalar,float>::value)
- Dtype = SLU_S;
- else if (internal::is_same<Scalar,double>::value)
- Dtype = SLU_D;
- else if (internal::is_same<Scalar,std::complex<float> >::value)
- Dtype = SLU_C;
- else if (internal::is_same<Scalar,std::complex<double> >::value)
- Dtype = SLU_Z;
- else
- {
- eigen_assert(false && "Scalar type not supported by SuperLU");
- }
- }
-
- template<typename MatrixType>
- static SluMatrix Map(MatrixBase<MatrixType>& _mat)
- {
- MatrixType& mat(_mat.derived());
- eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
- SluMatrix res;
- res.setStorageType(SLU_DN);
- res.setScalarType<typename MatrixType::Scalar>();
- res.Mtype = SLU_GE;
-
- res.nrow = mat.rows();
- res.ncol = mat.cols();
-
- res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
- res.storage.values = (void*)(mat.data());
- return res;
- }
-
- template<typename MatrixType>
- static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
- {
- SluMatrix res;
- if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
- {
- res.setStorageType(SLU_NR);
- res.nrow = mat.cols();
- res.ncol = mat.rows();
- }
- else
- {
- res.setStorageType(SLU_NC);
- res.nrow = mat.rows();
- res.ncol = mat.cols();
- }
-
- res.Mtype = SLU_GE;
-
- res.storage.nnz = mat.nonZeros();
- res.storage.values = mat.derived().valuePtr();
- res.storage.innerInd = mat.derived().innerIndexPtr();
- res.storage.outerInd = mat.derived().outerIndexPtr();
-
- res.setScalarType<typename MatrixType::Scalar>();
-
- // FIXME the following is not very accurate
- if (MatrixType::Flags & Upper)
- res.Mtype = SLU_TRU;
- if (MatrixType::Flags & Lower)
- res.Mtype = SLU_TRL;
-
- eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
-
- return res;
- }
-};
-
-template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
-struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
-{
- typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
- static void run(MatrixType& mat, SluMatrix& res)
- {
- eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
- res.setStorageType(SLU_DN);
- res.setScalarType<Scalar>();
- res.Mtype = SLU_GE;
-
- res.nrow = mat.rows();
- res.ncol = mat.cols();
-
- res.storage.lda = mat.outerStride();
- res.storage.values = mat.data();
- }
-};
-
-template<typename Derived>
-struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
-{
- typedef Derived MatrixType;
- static void run(MatrixType& mat, SluMatrix& res)
- {
- if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
- {
- res.setStorageType(SLU_NR);
- res.nrow = mat.cols();
- res.ncol = mat.rows();
- }
- else
- {
- res.setStorageType(SLU_NC);
- res.nrow = mat.rows();
- res.ncol = mat.cols();
- }
-
- res.Mtype = SLU_GE;
-
- res.storage.nnz = mat.nonZeros();
- res.storage.values = mat.valuePtr();
- res.storage.innerInd = mat.innerIndexPtr();
- res.storage.outerInd = mat.outerIndexPtr();
-
- res.setScalarType<typename MatrixType::Scalar>();
-
- // FIXME the following is not very accurate
- if (MatrixType::Flags & Upper)
- res.Mtype = SLU_TRU;
- if (MatrixType::Flags & Lower)
- res.Mtype = SLU_TRL;
-
- eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
- }
-};
-
-namespace internal {
-
-template<typename MatrixType>
-SluMatrix asSluMatrix(MatrixType& mat)
-{
- return SluMatrix::Map(mat);
-}
-
-/** View a Super LU matrix as an Eigen expression */
-template<typename Scalar, int Flags, typename Index>
-MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
-{
- eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
- || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
-
- Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
-
- return MappedSparseMatrix<Scalar,Flags,Index>(
- sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
- sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
-}
-
-} // end namespace internal
-
-/** \ingroup SuperLUSupport_Module
- * \class SuperLUBase
- * \brief The base class for the direct and incomplete LU factorization of SuperLU
- */
-template<typename _MatrixType, typename Derived>
-class SuperLUBase : internal::noncopyable
-{
- public:
- typedef _MatrixType MatrixType;
- typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef Matrix<Scalar,Dynamic,1> Vector;
- typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
- typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
- typedef SparseMatrix<Scalar> LUMatrixType;
-
- public:
-
- SuperLUBase() {}
-
- ~SuperLUBase()
- {
- clearFactors();
- }
-
- Derived& derived() { return *static_cast<Derived*>(this); }
- const Derived& derived() const { return *static_cast<const Derived*>(this); }
-
- inline Index rows() const { return m_matrix.rows(); }
- inline Index cols() const { return m_matrix.cols(); }
-
- /** \returns a reference to the Super LU option object to configure the Super LU algorithms. */
- inline superlu_options_t& options() { return m_sluOptions; }
-
- /** \brief Reports whether previous computation was successful.
- *
- * \returns \c Success if computation was succesful,
- * \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 */
- void compute(const MatrixType& matrix)
- {
- derived().analyzePattern(matrix);
- derived().factorize(matrix);
- }
-
- /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
- *
- * \sa compute()
- */
- template<typename Rhs>
- inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<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::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
- }
-
- /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
- *
- * \sa compute()
- */
- 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.
- *
- * This function is particularly useful when solving for several problems having the same structure.
- *
- * \sa factorize()
- */
- void analyzePattern(const MatrixType& /*matrix*/)
- {
- m_isInitialized = true;
- m_info = Success;
- m_analysisIsOk = true;
- m_factorizationIsOk = false;
- }
-
- template<typename Stream>
- void dumpMemory(Stream& /*s*/)
- {}
-
- protected:
-
- void initFactorization(const MatrixType& a)
- {
- set_default_options(&this->m_sluOptions);
-
- const int size = a.rows();
- m_matrix = a;
-
- m_sluA = internal::asSluMatrix(m_matrix);
- clearFactors();
-
- m_p.resize(size);
- m_q.resize(size);
- m_sluRscale.resize(size);
- m_sluCscale.resize(size);
- m_sluEtree.resize(size);
-
- // set empty B and X
- m_sluB.setStorageType(SLU_DN);
- m_sluB.setScalarType<Scalar>();
- m_sluB.Mtype = SLU_GE;
- m_sluB.storage.values = 0;
- m_sluB.nrow = 0;
- m_sluB.ncol = 0;
- m_sluB.storage.lda = size;
- m_sluX = m_sluB;
-
- m_extractedDataAreDirty = true;
- }
-
- void init()
- {
- m_info = InvalidInput;
- m_isInitialized = false;
- m_sluL.Store = 0;
- m_sluU.Store = 0;
- }
-
- void extractData() const;
-
- void clearFactors()
- {
- if(m_sluL.Store)
- Destroy_SuperNode_Matrix(&m_sluL);
- if(m_sluU.Store)
- Destroy_CompCol_Matrix(&m_sluU);
-
- m_sluL.Store = 0;
- m_sluU.Store = 0;
-
- memset(&m_sluL,0,sizeof m_sluL);
- memset(&m_sluU,0,sizeof m_sluU);
- }
-
- // cached data to reduce reallocation, etc.
- mutable LUMatrixType m_l;
- mutable LUMatrixType m_u;
- mutable IntColVectorType m_p;
- mutable IntRowVectorType m_q;
-
- mutable LUMatrixType m_matrix; // copy of the factorized matrix
- mutable SluMatrix m_sluA;
- mutable SuperMatrix m_sluL, m_sluU;
- mutable SluMatrix m_sluB, m_sluX;
- mutable SuperLUStat_t m_sluStat;
- mutable superlu_options_t m_sluOptions;
- mutable std::vector<int> m_sluEtree;
- mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
- mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
- mutable char m_sluEqued;
-
- mutable ComputationInfo m_info;
- bool m_isInitialized;
- int m_factorizationIsOk;
- int m_analysisIsOk;
- mutable bool m_extractedDataAreDirty;
-
- private:
- SuperLUBase(SuperLUBase& ) { }
-};
-
-
-/** \ingroup SuperLUSupport_Module
- * \class SuperLU
- * \brief A sparse direct LU factorization and solver based on the SuperLU library
- *
- * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
- * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
- * X and B can be either dense or sparse.
- *
- * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
- *
- * \sa \ref TutorialSparseDirectSolvers
- */
-template<typename _MatrixType>
-class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
-{
- public:
- typedef SuperLUBase<_MatrixType,SuperLU> Base;
- typedef _MatrixType MatrixType;
- typedef typename Base::Scalar Scalar;
- typedef typename Base::RealScalar RealScalar;
- typedef typename Base::Index Index;
- typedef typename Base::IntRowVectorType IntRowVectorType;
- typedef typename Base::IntColVectorType IntColVectorType;
- typedef typename Base::LUMatrixType LUMatrixType;
- typedef TriangularView<LUMatrixType, Lower|UnitDiag> LMatrixType;
- typedef TriangularView<LUMatrixType, Upper> UMatrixType;
-
- public:
-
- SuperLU() : Base() { init(); }
-
- SuperLU(const MatrixType& matrix) : Base()
- {
- init();
- Base::compute(matrix);
- }
-
- ~SuperLU()
- {
- }
-
- /** 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()
- */
- void analyzePattern(const MatrixType& matrix)
- {
- m_info = InvalidInput;
- m_isInitialized = false;
- Base::analyzePattern(matrix);
- }
-
- /** Performs a numeric decomposition of \a matrix
- *
- * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
- *
- * \sa analyzePattern()
- */
- void factorize(const MatrixType& matrix);
-
- #ifndef EIGEN_PARSED_BY_DOXYGEN
- /** \internal */
- template<typename Rhs,typename Dest>
- void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
- #endif // EIGEN_PARSED_BY_DOXYGEN
-
- inline const LMatrixType& matrixL() const
- {
- if (m_extractedDataAreDirty) this->extractData();
- return m_l;
- }
-
- inline const UMatrixType& matrixU() const
- {
- if (m_extractedDataAreDirty) this->extractData();
- return m_u;
- }
-
- inline const IntColVectorType& permutationP() const
- {
- if (m_extractedDataAreDirty) this->extractData();
- return m_p;
- }
-
- inline const IntRowVectorType& permutationQ() const
- {
- if (m_extractedDataAreDirty) this->extractData();
- return m_q;
- }
-
- Scalar determinant() const;
-
- protected:
-
- using Base::m_matrix;
- using Base::m_sluOptions;
- using Base::m_sluA;
- using Base::m_sluB;
- using Base::m_sluX;
- using Base::m_p;
- using Base::m_q;
- using Base::m_sluEtree;
- using Base::m_sluEqued;
- using Base::m_sluRscale;
- using Base::m_sluCscale;
- using Base::m_sluL;
- using Base::m_sluU;
- using Base::m_sluStat;
- using Base::m_sluFerr;
- using Base::m_sluBerr;
- using Base::m_l;
- using Base::m_u;
-
- using Base::m_analysisIsOk;
- using Base::m_factorizationIsOk;
- using Base::m_extractedDataAreDirty;
- using Base::m_isInitialized;
- using Base::m_info;
-
- void init()
- {
- Base::init();
-
- set_default_options(&this->m_sluOptions);
- m_sluOptions.PrintStat = NO;
- m_sluOptions.ConditionNumber = NO;
- m_sluOptions.Trans = NOTRANS;
- m_sluOptions.ColPerm = COLAMD;
- }
-
-
- private:
- SuperLU(SuperLU& ) { }
-};
-
-template<typename MatrixType>
-void SuperLU<MatrixType>::factorize(const MatrixType& a)
-{
- eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
- if(!m_analysisIsOk)
- {
- m_info = InvalidInput;
- return;
- }
-
- this->initFactorization(a);
-
- m_sluOptions.ColPerm = COLAMD;
- int info = 0;
- RealScalar recip_pivot_growth, rcond;
- RealScalar ferr, berr;
-
- StatInit(&m_sluStat);
- SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
- &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
- &m_sluL, &m_sluU,
- NULL, 0,
- &m_sluB, &m_sluX,
- &recip_pivot_growth, &rcond,
- &ferr, &berr,
- &m_sluStat, &info, Scalar());
- StatFree(&m_sluStat);
-
- m_extractedDataAreDirty = true;
-
- // FIXME how to better check for errors ???
- m_info = info == 0 ? Success : NumericalIssue;
- m_factorizationIsOk = true;
-}
-
-template<typename MatrixType>
-template<typename Rhs,typename Dest>
-void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
-{
- eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
-
- const int size = m_matrix.rows();
- const int rhsCols = b.cols();
- eigen_assert(size==b.rows());
-
- m_sluOptions.Trans = NOTRANS;
- m_sluOptions.Fact = FACTORED;
- m_sluOptions.IterRefine = NOREFINE;
-
-
- m_sluFerr.resize(rhsCols);
- m_sluBerr.resize(rhsCols);
- m_sluB = SluMatrix::Map(b.const_cast_derived());
- m_sluX = SluMatrix::Map(x.derived());
-
- typename Rhs::PlainObject b_cpy;
- if(m_sluEqued!='N')
- {
- b_cpy = b;
- m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
- }
-
- StatInit(&m_sluStat);
- int info = 0;
- RealScalar recip_pivot_growth, rcond;
- SuperLU_gssvx(&m_sluOptions, &m_sluA,
- m_q.data(), m_p.data(),
- &m_sluEtree[0], &m_sluEqued,
- &m_sluRscale[0], &m_sluCscale[0],
- &m_sluL, &m_sluU,
- NULL, 0,
- &m_sluB, &m_sluX,
- &recip_pivot_growth, &rcond,
- &m_sluFerr[0], &m_sluBerr[0],
- &m_sluStat, &info, Scalar());
- StatFree(&m_sluStat);
- m_info = info==0 ? Success : NumericalIssue;
-}
-
-// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
-//
-// Copyright (c) 1994 by Xerox Corporation. All rights reserved.
-//
-// THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
-// EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
-//
-template<typename MatrixType, typename Derived>
-void SuperLUBase<MatrixType,Derived>::extractData() const
-{
- eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
- if (m_extractedDataAreDirty)
- {
- int upper;
- int fsupc, istart, nsupr;
- int lastl = 0, lastu = 0;
- SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
- NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
- Scalar *SNptr;
-
- const int size = m_matrix.rows();
- m_l.resize(size,size);
- m_l.resizeNonZeros(Lstore->nnz);
- m_u.resize(size,size);
- m_u.resizeNonZeros(Ustore->nnz);
-
- int* Lcol = m_l.outerIndexPtr();
- int* Lrow = m_l.innerIndexPtr();
- Scalar* Lval = m_l.valuePtr();
-
- int* Ucol = m_u.outerIndexPtr();
- int* Urow = m_u.innerIndexPtr();
- Scalar* Uval = m_u.valuePtr();
-
- Ucol[0] = 0;
- Ucol[0] = 0;
-
- /* for each supernode */
- for (int k = 0; k <= Lstore->nsuper; ++k)
- {
- fsupc = L_FST_SUPC(k);
- istart = L_SUB_START(fsupc);
- nsupr = L_SUB_START(fsupc+1) - istart;
- upper = 1;
-
- /* for each column in the supernode */
- for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
- {
- SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
-
- /* Extract U */
- for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
- {
- Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
- /* Matlab doesn't like explicit zero. */
- if (Uval[lastu] != 0.0)
- Urow[lastu++] = U_SUB(i);
- }
- for (int i = 0; i < upper; ++i)
- {
- /* upper triangle in the supernode */
- Uval[lastu] = SNptr[i];
- /* Matlab doesn't like explicit zero. */
- if (Uval[lastu] != 0.0)
- Urow[lastu++] = L_SUB(istart+i);
- }
- Ucol[j+1] = lastu;
-
- /* Extract L */
- Lval[lastl] = 1.0; /* unit diagonal */
- Lrow[lastl++] = L_SUB(istart + upper - 1);
- for (int i = upper; i < nsupr; ++i)
- {
- Lval[lastl] = SNptr[i];
- /* Matlab doesn't like explicit zero. */
- if (Lval[lastl] != 0.0)
- Lrow[lastl++] = L_SUB(istart+i);
- }
- Lcol[j+1] = lastl;
-
- ++upper;
- } /* for j ... */
-
- } /* for k ... */
-
- // squeeze the matrices :
- m_l.resizeNonZeros(lastl);
- m_u.resizeNonZeros(lastu);
-
- m_extractedDataAreDirty = false;
- }
-}
-
-template<typename MatrixType>
-typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
-{
- eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
-
- if (m_extractedDataAreDirty)
- this->extractData();
-
- Scalar det = Scalar(1);
- for (int j=0; j<m_u.cols(); ++j)
- {
- if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
- {
- int lastId = m_u.outerIndexPtr()[j+1]-1;
- eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
- if (m_u.innerIndexPtr()[lastId]==j)
- det *= m_u.valuePtr()[lastId];
- }
- }
- if(m_sluEqued!='N')
- return det/m_sluRscale.prod()/m_sluCscale.prod();
- else
- return det;
-}
-
-#ifdef EIGEN_PARSED_BY_DOXYGEN
-#define EIGEN_SUPERLU_HAS_ILU
-#endif
-
-#ifdef EIGEN_SUPERLU_HAS_ILU
-
-/** \ingroup SuperLUSupport_Module
- * \class SuperILU
- * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
- *
- * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
- * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
- *
- * \warning This class requires SuperLU 4 or later.
- *
- * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
- *
- * \sa \ref TutorialSparseDirectSolvers, class ConjugateGradient, class BiCGSTAB
- */
-
-template<typename _MatrixType>
-class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
-{
- public:
- typedef SuperLUBase<_MatrixType,SuperILU> Base;
- typedef _MatrixType MatrixType;
- typedef typename Base::Scalar Scalar;
- typedef typename Base::RealScalar RealScalar;
- typedef typename Base::Index Index;
-
- public:
-
- SuperILU() : Base() { init(); }
-
- SuperILU(const MatrixType& matrix) : Base()
- {
- init();
- Base::compute(matrix);
- }
-
- ~SuperILU()
- {
- }
-
- /** 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()
- */
- void analyzePattern(const MatrixType& matrix)
- {
- Base::analyzePattern(matrix);
- }
-
- /** Performs a numeric decomposition of \a matrix
- *
- * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
- *
- * \sa analyzePattern()
- */
- void factorize(const MatrixType& matrix);
-
- #ifndef EIGEN_PARSED_BY_DOXYGEN
- /** \internal */
- template<typename Rhs,typename Dest>
- void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
- #endif // EIGEN_PARSED_BY_DOXYGEN
-
- protected:
-
- using Base::m_matrix;
- using Base::m_sluOptions;
- using Base::m_sluA;
- using Base::m_sluB;
- using Base::m_sluX;
- using Base::m_p;
- using Base::m_q;
- using Base::m_sluEtree;
- using Base::m_sluEqued;
- using Base::m_sluRscale;
- using Base::m_sluCscale;
- using Base::m_sluL;
- using Base::m_sluU;
- using Base::m_sluStat;
- using Base::m_sluFerr;
- using Base::m_sluBerr;
- using Base::m_l;
- using Base::m_u;
-
- using Base::m_analysisIsOk;
- using Base::m_factorizationIsOk;
- using Base::m_extractedDataAreDirty;
- using Base::m_isInitialized;
- using Base::m_info;
-
- void init()
- {
- Base::init();
-
- ilu_set_default_options(&m_sluOptions);
- m_sluOptions.PrintStat = NO;
- m_sluOptions.ConditionNumber = NO;
- m_sluOptions.Trans = NOTRANS;
- m_sluOptions.ColPerm = MMD_AT_PLUS_A;
-
- // no attempt to preserve column sum
- m_sluOptions.ILU_MILU = SILU;
- // only basic ILU(k) support -- no direct control over memory consumption
- // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
- // and set ILU_FillFactor to max memory growth
- m_sluOptions.ILU_DropRule = DROP_BASIC;
- m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
- }
-
- private:
- SuperILU(SuperILU& ) { }
-};
-
-template<typename MatrixType>
-void SuperILU<MatrixType>::factorize(const MatrixType& a)
-{
- eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
- if(!m_analysisIsOk)
- {
- m_info = InvalidInput;
- return;
- }
-
- this->initFactorization(a);
-
- int info = 0;
- RealScalar recip_pivot_growth, rcond;
-
- StatInit(&m_sluStat);
- SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
- &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
- &m_sluL, &m_sluU,
- NULL, 0,
- &m_sluB, &m_sluX,
- &recip_pivot_growth, &rcond,
- &m_sluStat, &info, Scalar());
- StatFree(&m_sluStat);
-
- // FIXME how to better check for errors ???
- m_info = info == 0 ? Success : NumericalIssue;
- m_factorizationIsOk = true;
-}
-
-template<typename MatrixType>
-template<typename Rhs,typename Dest>
-void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
-{
- eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
-
- const int size = m_matrix.rows();
- const int rhsCols = b.cols();
- eigen_assert(size==b.rows());
-
- m_sluOptions.Trans = NOTRANS;
- m_sluOptions.Fact = FACTORED;
- m_sluOptions.IterRefine = NOREFINE;
-
- m_sluFerr.resize(rhsCols);
- m_sluBerr.resize(rhsCols);
- m_sluB = SluMatrix::Map(b.const_cast_derived());
- m_sluX = SluMatrix::Map(x.derived());
-
- typename Rhs::PlainObject b_cpy;
- if(m_sluEqued!='N')
- {
- b_cpy = b;
- m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
- }
-
- int info = 0;
- RealScalar recip_pivot_growth, rcond;
-
- StatInit(&m_sluStat);
- SuperLU_gsisx(&m_sluOptions, &m_sluA,
- m_q.data(), m_p.data(),
- &m_sluEtree[0], &m_sluEqued,
- &m_sluRscale[0], &m_sluCscale[0],
- &m_sluL, &m_sluU,
- NULL, 0,
- &m_sluB, &m_sluX,
- &recip_pivot_growth, &rcond,
- &m_sluStat, &info, Scalar());
- StatFree(&m_sluStat);
-
- m_info = info==0 ? Success : NumericalIssue;
-}
-#endif
-
-namespace internal {
-
-template<typename _MatrixType, typename Derived, typename Rhs>
-struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
- : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
-{
- typedef SuperLUBase<_MatrixType,Derived> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec().derived()._solve(rhs(),dst);
- }
-};
-
-template<typename _MatrixType, typename Derived, typename Rhs>
-struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
- : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
-{
- typedef SuperLUBase<_MatrixType,Derived> Dec;
- EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- this->defaultEvalTo(dst);
- }
-};
-
-} // end namespace internal
-
-} // end namespace Eigen
-
-#endif // EIGEN_SUPERLUSUPPORT_H