Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> |
| 5 | // |
| 6 | // This Source Code Form is subject to the terms of the Mozilla |
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed |
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 9 | |
| 10 | #ifndef EIGEN_SUPERLUSUPPORT_H |
| 11 | #define EIGEN_SUPERLUSUPPORT_H |
| 12 | |
| 13 | namespace Eigen { |
| 14 | |
| 15 | #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE) \ |
| 16 | extern "C" { \ |
| 17 | typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t; \ |
| 18 | extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \ |
| 19 | char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \ |
| 20 | void *, int, SuperMatrix *, SuperMatrix *, \ |
| 21 | FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, \ |
| 22 | PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \ |
| 23 | } \ |
| 24 | inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A, \ |
| 25 | int *perm_c, int *perm_r, int *etree, char *equed, \ |
| 26 | FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \ |
| 27 | SuperMatrix *U, void *work, int lwork, \ |
| 28 | SuperMatrix *B, SuperMatrix *X, \ |
| 29 | FLOATTYPE *recip_pivot_growth, \ |
| 30 | FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \ |
| 31 | SuperLUStat_t *stats, int *info, KEYTYPE) { \ |
| 32 | PREFIX##mem_usage_t mem_usage; \ |
| 33 | PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L, \ |
| 34 | U, work, lwork, B, X, recip_pivot_growth, rcond, \ |
| 35 | ferr, berr, &mem_usage, stats, info); \ |
| 36 | return mem_usage.for_lu; /* bytes used by the factor storage */ \ |
| 37 | } |
| 38 | |
| 39 | DECL_GSSVX(s,float,float) |
| 40 | DECL_GSSVX(c,float,std::complex<float>) |
| 41 | DECL_GSSVX(d,double,double) |
| 42 | DECL_GSSVX(z,double,std::complex<double>) |
| 43 | |
| 44 | #ifdef MILU_ALPHA |
| 45 | #define EIGEN_SUPERLU_HAS_ILU |
| 46 | #endif |
| 47 | |
| 48 | #ifdef EIGEN_SUPERLU_HAS_ILU |
| 49 | |
| 50 | // similarly for the incomplete factorization using gsisx |
| 51 | #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE) \ |
| 52 | extern "C" { \ |
| 53 | extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *, \ |
| 54 | char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *, \ |
| 55 | void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *, \ |
| 56 | PREFIX##mem_usage_t *, SuperLUStat_t *, int *); \ |
| 57 | } \ |
| 58 | inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \ |
| 59 | int *perm_c, int *perm_r, int *etree, char *equed, \ |
| 60 | FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \ |
| 61 | SuperMatrix *U, void *work, int lwork, \ |
| 62 | SuperMatrix *B, SuperMatrix *X, \ |
| 63 | FLOATTYPE *recip_pivot_growth, \ |
| 64 | FLOATTYPE *rcond, \ |
| 65 | SuperLUStat_t *stats, int *info, KEYTYPE) { \ |
| 66 | PREFIX##mem_usage_t mem_usage; \ |
| 67 | PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L, \ |
| 68 | U, work, lwork, B, X, recip_pivot_growth, rcond, \ |
| 69 | &mem_usage, stats, info); \ |
| 70 | return mem_usage.for_lu; /* bytes used by the factor storage */ \ |
| 71 | } |
| 72 | |
| 73 | DECL_GSISX(s,float,float) |
| 74 | DECL_GSISX(c,float,std::complex<float>) |
| 75 | DECL_GSISX(d,double,double) |
| 76 | DECL_GSISX(z,double,std::complex<double>) |
| 77 | |
| 78 | #endif |
| 79 | |
| 80 | template<typename MatrixType> |
| 81 | struct SluMatrixMapHelper; |
| 82 | |
| 83 | /** \internal |
| 84 | * |
| 85 | * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices |
| 86 | * and dense matrices. Supernodal and other fancy format are not supported by this wrapper. |
| 87 | * |
| 88 | * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure. |
| 89 | */ |
| 90 | struct SluMatrix : SuperMatrix |
| 91 | { |
| 92 | SluMatrix() |
| 93 | { |
| 94 | Store = &storage; |
| 95 | } |
| 96 | |
| 97 | SluMatrix(const SluMatrix& other) |
| 98 | : SuperMatrix(other) |
| 99 | { |
| 100 | Store = &storage; |
| 101 | storage = other.storage; |
| 102 | } |
| 103 | |
| 104 | SluMatrix& operator=(const SluMatrix& other) |
| 105 | { |
| 106 | SuperMatrix::operator=(static_cast<const SuperMatrix&>(other)); |
| 107 | Store = &storage; |
| 108 | storage = other.storage; |
| 109 | return *this; |
| 110 | } |
| 111 | |
| 112 | struct |
| 113 | { |
| 114 | union {int nnz;int lda;}; |
| 115 | void *values; |
| 116 | int *innerInd; |
| 117 | int *outerInd; |
| 118 | } storage; |
| 119 | |
| 120 | void setStorageType(Stype_t t) |
| 121 | { |
| 122 | Stype = t; |
| 123 | if (t==SLU_NC || t==SLU_NR || t==SLU_DN) |
| 124 | Store = &storage; |
| 125 | else |
| 126 | { |
| 127 | eigen_assert(false && "storage type not supported"); |
| 128 | Store = 0; |
| 129 | } |
| 130 | } |
| 131 | |
| 132 | template<typename Scalar> |
| 133 | void setScalarType() |
| 134 | { |
| 135 | if (internal::is_same<Scalar,float>::value) |
| 136 | Dtype = SLU_S; |
| 137 | else if (internal::is_same<Scalar,double>::value) |
| 138 | Dtype = SLU_D; |
| 139 | else if (internal::is_same<Scalar,std::complex<float> >::value) |
| 140 | Dtype = SLU_C; |
| 141 | else if (internal::is_same<Scalar,std::complex<double> >::value) |
| 142 | Dtype = SLU_Z; |
| 143 | else |
| 144 | { |
| 145 | eigen_assert(false && "Scalar type not supported by SuperLU"); |
| 146 | } |
| 147 | } |
| 148 | |
| 149 | template<typename MatrixType> |
| 150 | static SluMatrix Map(MatrixBase<MatrixType>& _mat) |
| 151 | { |
| 152 | MatrixType& mat(_mat.derived()); |
| 153 | eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU"); |
| 154 | SluMatrix res; |
| 155 | res.setStorageType(SLU_DN); |
| 156 | res.setScalarType<typename MatrixType::Scalar>(); |
| 157 | res.Mtype = SLU_GE; |
| 158 | |
| 159 | res.nrow = mat.rows(); |
| 160 | res.ncol = mat.cols(); |
| 161 | |
| 162 | res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 163 | res.storage.values = (void*)(mat.data()); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 164 | return res; |
| 165 | } |
| 166 | |
| 167 | template<typename MatrixType> |
| 168 | static SluMatrix Map(SparseMatrixBase<MatrixType>& mat) |
| 169 | { |
| 170 | SluMatrix res; |
| 171 | if ((MatrixType::Flags&RowMajorBit)==RowMajorBit) |
| 172 | { |
| 173 | res.setStorageType(SLU_NR); |
| 174 | res.nrow = mat.cols(); |
| 175 | res.ncol = mat.rows(); |
| 176 | } |
| 177 | else |
| 178 | { |
| 179 | res.setStorageType(SLU_NC); |
| 180 | res.nrow = mat.rows(); |
| 181 | res.ncol = mat.cols(); |
| 182 | } |
| 183 | |
| 184 | res.Mtype = SLU_GE; |
| 185 | |
| 186 | res.storage.nnz = mat.nonZeros(); |
| 187 | res.storage.values = mat.derived().valuePtr(); |
| 188 | res.storage.innerInd = mat.derived().innerIndexPtr(); |
| 189 | res.storage.outerInd = mat.derived().outerIndexPtr(); |
| 190 | |
| 191 | res.setScalarType<typename MatrixType::Scalar>(); |
| 192 | |
| 193 | // FIXME the following is not very accurate |
| 194 | if (MatrixType::Flags & Upper) |
| 195 | res.Mtype = SLU_TRU; |
| 196 | if (MatrixType::Flags & Lower) |
| 197 | res.Mtype = SLU_TRL; |
| 198 | |
| 199 | eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU"); |
| 200 | |
| 201 | return res; |
| 202 | } |
| 203 | }; |
| 204 | |
| 205 | template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols> |
| 206 | struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> > |
| 207 | { |
| 208 | typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType; |
| 209 | static void run(MatrixType& mat, SluMatrix& res) |
| 210 | { |
| 211 | eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU"); |
| 212 | res.setStorageType(SLU_DN); |
| 213 | res.setScalarType<Scalar>(); |
| 214 | res.Mtype = SLU_GE; |
| 215 | |
| 216 | res.nrow = mat.rows(); |
| 217 | res.ncol = mat.cols(); |
| 218 | |
| 219 | res.storage.lda = mat.outerStride(); |
| 220 | res.storage.values = mat.data(); |
| 221 | } |
| 222 | }; |
| 223 | |
| 224 | template<typename Derived> |
| 225 | struct SluMatrixMapHelper<SparseMatrixBase<Derived> > |
| 226 | { |
| 227 | typedef Derived MatrixType; |
| 228 | static void run(MatrixType& mat, SluMatrix& res) |
| 229 | { |
| 230 | if ((MatrixType::Flags&RowMajorBit)==RowMajorBit) |
| 231 | { |
| 232 | res.setStorageType(SLU_NR); |
| 233 | res.nrow = mat.cols(); |
| 234 | res.ncol = mat.rows(); |
| 235 | } |
| 236 | else |
| 237 | { |
| 238 | res.setStorageType(SLU_NC); |
| 239 | res.nrow = mat.rows(); |
| 240 | res.ncol = mat.cols(); |
| 241 | } |
| 242 | |
| 243 | res.Mtype = SLU_GE; |
| 244 | |
| 245 | res.storage.nnz = mat.nonZeros(); |
| 246 | res.storage.values = mat.valuePtr(); |
| 247 | res.storage.innerInd = mat.innerIndexPtr(); |
| 248 | res.storage.outerInd = mat.outerIndexPtr(); |
| 249 | |
| 250 | res.setScalarType<typename MatrixType::Scalar>(); |
| 251 | |
| 252 | // FIXME the following is not very accurate |
| 253 | if (MatrixType::Flags & Upper) |
| 254 | res.Mtype = SLU_TRU; |
| 255 | if (MatrixType::Flags & Lower) |
| 256 | res.Mtype = SLU_TRL; |
| 257 | |
| 258 | eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU"); |
| 259 | } |
| 260 | }; |
| 261 | |
| 262 | namespace internal { |
| 263 | |
| 264 | template<typename MatrixType> |
| 265 | SluMatrix asSluMatrix(MatrixType& mat) |
| 266 | { |
| 267 | return SluMatrix::Map(mat); |
| 268 | } |
| 269 | |
| 270 | /** View a Super LU matrix as an Eigen expression */ |
| 271 | template<typename Scalar, int Flags, typename Index> |
| 272 | MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat) |
| 273 | { |
| 274 | eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR |
| 275 | || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC); |
| 276 | |
| 277 | Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow; |
| 278 | |
| 279 | return MappedSparseMatrix<Scalar,Flags,Index>( |
| 280 | sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize], |
| 281 | sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) ); |
| 282 | } |
| 283 | |
| 284 | } // end namespace internal |
| 285 | |
| 286 | /** \ingroup SuperLUSupport_Module |
| 287 | * \class SuperLUBase |
| 288 | * \brief The base class for the direct and incomplete LU factorization of SuperLU |
| 289 | */ |
| 290 | template<typename _MatrixType, typename Derived> |
| 291 | class SuperLUBase : internal::noncopyable |
| 292 | { |
| 293 | public: |
| 294 | typedef _MatrixType MatrixType; |
| 295 | typedef typename MatrixType::Scalar Scalar; |
| 296 | typedef typename MatrixType::RealScalar RealScalar; |
| 297 | typedef typename MatrixType::Index Index; |
| 298 | typedef Matrix<Scalar,Dynamic,1> Vector; |
| 299 | typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; |
| 300 | typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; |
| 301 | typedef SparseMatrix<Scalar> LUMatrixType; |
| 302 | |
| 303 | public: |
| 304 | |
| 305 | SuperLUBase() {} |
| 306 | |
| 307 | ~SuperLUBase() |
| 308 | { |
| 309 | clearFactors(); |
| 310 | } |
| 311 | |
| 312 | Derived& derived() { return *static_cast<Derived*>(this); } |
| 313 | const Derived& derived() const { return *static_cast<const Derived*>(this); } |
| 314 | |
| 315 | inline Index rows() const { return m_matrix.rows(); } |
| 316 | inline Index cols() const { return m_matrix.cols(); } |
| 317 | |
| 318 | /** \returns a reference to the Super LU option object to configure the Super LU algorithms. */ |
| 319 | inline superlu_options_t& options() { return m_sluOptions; } |
| 320 | |
| 321 | /** \brief Reports whether previous computation was successful. |
| 322 | * |
| 323 | * \returns \c Success if computation was succesful, |
| 324 | * \c NumericalIssue if the matrix.appears to be negative. |
| 325 | */ |
| 326 | ComputationInfo info() const |
| 327 | { |
| 328 | eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| 329 | return m_info; |
| 330 | } |
| 331 | |
| 332 | /** Computes the sparse Cholesky decomposition of \a matrix */ |
| 333 | void compute(const MatrixType& matrix) |
| 334 | { |
| 335 | derived().analyzePattern(matrix); |
| 336 | derived().factorize(matrix); |
| 337 | } |
| 338 | |
| 339 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. |
| 340 | * |
| 341 | * \sa compute() |
| 342 | */ |
| 343 | template<typename Rhs> |
| 344 | inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const |
| 345 | { |
| 346 | eigen_assert(m_isInitialized && "SuperLU is not initialized."); |
| 347 | eigen_assert(rows()==b.rows() |
| 348 | && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); |
| 349 | return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived()); |
| 350 | } |
| 351 | |
| 352 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. |
| 353 | * |
| 354 | * \sa compute() |
| 355 | */ |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 356 | template<typename Rhs> |
| 357 | inline const internal::sparse_solve_retval<SuperLUBase, Rhs> solve(const SparseMatrixBase<Rhs>& b) const |
| 358 | { |
| 359 | eigen_assert(m_isInitialized && "SuperLU is not initialized."); |
| 360 | eigen_assert(rows()==b.rows() |
| 361 | && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); |
| 362 | return internal::sparse_solve_retval<SuperLUBase, Rhs>(*this, b.derived()); |
| 363 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 364 | |
| 365 | /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| 366 | * |
| 367 | * This function is particularly useful when solving for several problems having the same structure. |
| 368 | * |
| 369 | * \sa factorize() |
| 370 | */ |
| 371 | void analyzePattern(const MatrixType& /*matrix*/) |
| 372 | { |
| 373 | m_isInitialized = true; |
| 374 | m_info = Success; |
| 375 | m_analysisIsOk = true; |
| 376 | m_factorizationIsOk = false; |
| 377 | } |
| 378 | |
| 379 | template<typename Stream> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 380 | void dumpMemory(Stream& /*s*/) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 381 | {} |
| 382 | |
| 383 | protected: |
| 384 | |
| 385 | void initFactorization(const MatrixType& a) |
| 386 | { |
| 387 | set_default_options(&this->m_sluOptions); |
| 388 | |
| 389 | const int size = a.rows(); |
| 390 | m_matrix = a; |
| 391 | |
| 392 | m_sluA = internal::asSluMatrix(m_matrix); |
| 393 | clearFactors(); |
| 394 | |
| 395 | m_p.resize(size); |
| 396 | m_q.resize(size); |
| 397 | m_sluRscale.resize(size); |
| 398 | m_sluCscale.resize(size); |
| 399 | m_sluEtree.resize(size); |
| 400 | |
| 401 | // set empty B and X |
| 402 | m_sluB.setStorageType(SLU_DN); |
| 403 | m_sluB.setScalarType<Scalar>(); |
| 404 | m_sluB.Mtype = SLU_GE; |
| 405 | m_sluB.storage.values = 0; |
| 406 | m_sluB.nrow = 0; |
| 407 | m_sluB.ncol = 0; |
| 408 | m_sluB.storage.lda = size; |
| 409 | m_sluX = m_sluB; |
| 410 | |
| 411 | m_extractedDataAreDirty = true; |
| 412 | } |
| 413 | |
| 414 | void init() |
| 415 | { |
| 416 | m_info = InvalidInput; |
| 417 | m_isInitialized = false; |
| 418 | m_sluL.Store = 0; |
| 419 | m_sluU.Store = 0; |
| 420 | } |
| 421 | |
| 422 | void extractData() const; |
| 423 | |
| 424 | void clearFactors() |
| 425 | { |
| 426 | if(m_sluL.Store) |
| 427 | Destroy_SuperNode_Matrix(&m_sluL); |
| 428 | if(m_sluU.Store) |
| 429 | Destroy_CompCol_Matrix(&m_sluU); |
| 430 | |
| 431 | m_sluL.Store = 0; |
| 432 | m_sluU.Store = 0; |
| 433 | |
| 434 | memset(&m_sluL,0,sizeof m_sluL); |
| 435 | memset(&m_sluU,0,sizeof m_sluU); |
| 436 | } |
| 437 | |
| 438 | // cached data to reduce reallocation, etc. |
| 439 | mutable LUMatrixType m_l; |
| 440 | mutable LUMatrixType m_u; |
| 441 | mutable IntColVectorType m_p; |
| 442 | mutable IntRowVectorType m_q; |
| 443 | |
| 444 | mutable LUMatrixType m_matrix; // copy of the factorized matrix |
| 445 | mutable SluMatrix m_sluA; |
| 446 | mutable SuperMatrix m_sluL, m_sluU; |
| 447 | mutable SluMatrix m_sluB, m_sluX; |
| 448 | mutable SuperLUStat_t m_sluStat; |
| 449 | mutable superlu_options_t m_sluOptions; |
| 450 | mutable std::vector<int> m_sluEtree; |
| 451 | mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale; |
| 452 | mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr; |
| 453 | mutable char m_sluEqued; |
| 454 | |
| 455 | mutable ComputationInfo m_info; |
| 456 | bool m_isInitialized; |
| 457 | int m_factorizationIsOk; |
| 458 | int m_analysisIsOk; |
| 459 | mutable bool m_extractedDataAreDirty; |
| 460 | |
| 461 | private: |
| 462 | SuperLUBase(SuperLUBase& ) { } |
| 463 | }; |
| 464 | |
| 465 | |
| 466 | /** \ingroup SuperLUSupport_Module |
| 467 | * \class SuperLU |
| 468 | * \brief A sparse direct LU factorization and solver based on the SuperLU library |
| 469 | * |
| 470 | * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization |
| 471 | * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices |
| 472 | * X and B can be either dense or sparse. |
| 473 | * |
| 474 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| 475 | * |
| 476 | * \sa \ref TutorialSparseDirectSolvers |
| 477 | */ |
| 478 | template<typename _MatrixType> |
| 479 | class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > |
| 480 | { |
| 481 | public: |
| 482 | typedef SuperLUBase<_MatrixType,SuperLU> Base; |
| 483 | typedef _MatrixType MatrixType; |
| 484 | typedef typename Base::Scalar Scalar; |
| 485 | typedef typename Base::RealScalar RealScalar; |
| 486 | typedef typename Base::Index Index; |
| 487 | typedef typename Base::IntRowVectorType IntRowVectorType; |
| 488 | typedef typename Base::IntColVectorType IntColVectorType; |
| 489 | typedef typename Base::LUMatrixType LUMatrixType; |
| 490 | typedef TriangularView<LUMatrixType, Lower|UnitDiag> LMatrixType; |
| 491 | typedef TriangularView<LUMatrixType, Upper> UMatrixType; |
| 492 | |
| 493 | public: |
| 494 | |
| 495 | SuperLU() : Base() { init(); } |
| 496 | |
| 497 | SuperLU(const MatrixType& matrix) : Base() |
| 498 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 499 | init(); |
| 500 | Base::compute(matrix); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 501 | } |
| 502 | |
| 503 | ~SuperLU() |
| 504 | { |
| 505 | } |
| 506 | |
| 507 | /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| 508 | * |
| 509 | * This function is particularly useful when solving for several problems having the same structure. |
| 510 | * |
| 511 | * \sa factorize() |
| 512 | */ |
| 513 | void analyzePattern(const MatrixType& matrix) |
| 514 | { |
| 515 | m_info = InvalidInput; |
| 516 | m_isInitialized = false; |
| 517 | Base::analyzePattern(matrix); |
| 518 | } |
| 519 | |
| 520 | /** Performs a numeric decomposition of \a matrix |
| 521 | * |
| 522 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| 523 | * |
| 524 | * \sa analyzePattern() |
| 525 | */ |
| 526 | void factorize(const MatrixType& matrix); |
| 527 | |
| 528 | #ifndef EIGEN_PARSED_BY_DOXYGEN |
| 529 | /** \internal */ |
| 530 | template<typename Rhs,typename Dest> |
| 531 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; |
| 532 | #endif // EIGEN_PARSED_BY_DOXYGEN |
| 533 | |
| 534 | inline const LMatrixType& matrixL() const |
| 535 | { |
| 536 | if (m_extractedDataAreDirty) this->extractData(); |
| 537 | return m_l; |
| 538 | } |
| 539 | |
| 540 | inline const UMatrixType& matrixU() const |
| 541 | { |
| 542 | if (m_extractedDataAreDirty) this->extractData(); |
| 543 | return m_u; |
| 544 | } |
| 545 | |
| 546 | inline const IntColVectorType& permutationP() const |
| 547 | { |
| 548 | if (m_extractedDataAreDirty) this->extractData(); |
| 549 | return m_p; |
| 550 | } |
| 551 | |
| 552 | inline const IntRowVectorType& permutationQ() const |
| 553 | { |
| 554 | if (m_extractedDataAreDirty) this->extractData(); |
| 555 | return m_q; |
| 556 | } |
| 557 | |
| 558 | Scalar determinant() const; |
| 559 | |
| 560 | protected: |
| 561 | |
| 562 | using Base::m_matrix; |
| 563 | using Base::m_sluOptions; |
| 564 | using Base::m_sluA; |
| 565 | using Base::m_sluB; |
| 566 | using Base::m_sluX; |
| 567 | using Base::m_p; |
| 568 | using Base::m_q; |
| 569 | using Base::m_sluEtree; |
| 570 | using Base::m_sluEqued; |
| 571 | using Base::m_sluRscale; |
| 572 | using Base::m_sluCscale; |
| 573 | using Base::m_sluL; |
| 574 | using Base::m_sluU; |
| 575 | using Base::m_sluStat; |
| 576 | using Base::m_sluFerr; |
| 577 | using Base::m_sluBerr; |
| 578 | using Base::m_l; |
| 579 | using Base::m_u; |
| 580 | |
| 581 | using Base::m_analysisIsOk; |
| 582 | using Base::m_factorizationIsOk; |
| 583 | using Base::m_extractedDataAreDirty; |
| 584 | using Base::m_isInitialized; |
| 585 | using Base::m_info; |
| 586 | |
| 587 | void init() |
| 588 | { |
| 589 | Base::init(); |
| 590 | |
| 591 | set_default_options(&this->m_sluOptions); |
| 592 | m_sluOptions.PrintStat = NO; |
| 593 | m_sluOptions.ConditionNumber = NO; |
| 594 | m_sluOptions.Trans = NOTRANS; |
| 595 | m_sluOptions.ColPerm = COLAMD; |
| 596 | } |
| 597 | |
| 598 | |
| 599 | private: |
| 600 | SuperLU(SuperLU& ) { } |
| 601 | }; |
| 602 | |
| 603 | template<typename MatrixType> |
| 604 | void SuperLU<MatrixType>::factorize(const MatrixType& a) |
| 605 | { |
| 606 | eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); |
| 607 | if(!m_analysisIsOk) |
| 608 | { |
| 609 | m_info = InvalidInput; |
| 610 | return; |
| 611 | } |
| 612 | |
| 613 | this->initFactorization(a); |
| 614 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 615 | m_sluOptions.ColPerm = COLAMD; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 616 | int info = 0; |
| 617 | RealScalar recip_pivot_growth, rcond; |
| 618 | RealScalar ferr, berr; |
| 619 | |
| 620 | StatInit(&m_sluStat); |
| 621 | SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], |
| 622 | &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], |
| 623 | &m_sluL, &m_sluU, |
| 624 | NULL, 0, |
| 625 | &m_sluB, &m_sluX, |
| 626 | &recip_pivot_growth, &rcond, |
| 627 | &ferr, &berr, |
| 628 | &m_sluStat, &info, Scalar()); |
| 629 | StatFree(&m_sluStat); |
| 630 | |
| 631 | m_extractedDataAreDirty = true; |
| 632 | |
| 633 | // FIXME how to better check for errors ??? |
| 634 | m_info = info == 0 ? Success : NumericalIssue; |
| 635 | m_factorizationIsOk = true; |
| 636 | } |
| 637 | |
| 638 | template<typename MatrixType> |
| 639 | template<typename Rhs,typename Dest> |
| 640 | void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const |
| 641 | { |
| 642 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); |
| 643 | |
| 644 | const int size = m_matrix.rows(); |
| 645 | const int rhsCols = b.cols(); |
| 646 | eigen_assert(size==b.rows()); |
| 647 | |
| 648 | m_sluOptions.Trans = NOTRANS; |
| 649 | m_sluOptions.Fact = FACTORED; |
| 650 | m_sluOptions.IterRefine = NOREFINE; |
| 651 | |
| 652 | |
| 653 | m_sluFerr.resize(rhsCols); |
| 654 | m_sluBerr.resize(rhsCols); |
| 655 | m_sluB = SluMatrix::Map(b.const_cast_derived()); |
| 656 | m_sluX = SluMatrix::Map(x.derived()); |
| 657 | |
| 658 | typename Rhs::PlainObject b_cpy; |
| 659 | if(m_sluEqued!='N') |
| 660 | { |
| 661 | b_cpy = b; |
| 662 | m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); |
| 663 | } |
| 664 | |
| 665 | StatInit(&m_sluStat); |
| 666 | int info = 0; |
| 667 | RealScalar recip_pivot_growth, rcond; |
| 668 | SuperLU_gssvx(&m_sluOptions, &m_sluA, |
| 669 | m_q.data(), m_p.data(), |
| 670 | &m_sluEtree[0], &m_sluEqued, |
| 671 | &m_sluRscale[0], &m_sluCscale[0], |
| 672 | &m_sluL, &m_sluU, |
| 673 | NULL, 0, |
| 674 | &m_sluB, &m_sluX, |
| 675 | &recip_pivot_growth, &rcond, |
| 676 | &m_sluFerr[0], &m_sluBerr[0], |
| 677 | &m_sluStat, &info, Scalar()); |
| 678 | StatFree(&m_sluStat); |
| 679 | m_info = info==0 ? Success : NumericalIssue; |
| 680 | } |
| 681 | |
| 682 | // the code of this extractData() function has been adapted from the SuperLU's Matlab support code, |
| 683 | // |
| 684 | // Copyright (c) 1994 by Xerox Corporation. All rights reserved. |
| 685 | // |
| 686 | // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY |
| 687 | // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. |
| 688 | // |
| 689 | template<typename MatrixType, typename Derived> |
| 690 | void SuperLUBase<MatrixType,Derived>::extractData() const |
| 691 | { |
| 692 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()"); |
| 693 | if (m_extractedDataAreDirty) |
| 694 | { |
| 695 | int upper; |
| 696 | int fsupc, istart, nsupr; |
| 697 | int lastl = 0, lastu = 0; |
| 698 | SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store); |
| 699 | NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store); |
| 700 | Scalar *SNptr; |
| 701 | |
| 702 | const int size = m_matrix.rows(); |
| 703 | m_l.resize(size,size); |
| 704 | m_l.resizeNonZeros(Lstore->nnz); |
| 705 | m_u.resize(size,size); |
| 706 | m_u.resizeNonZeros(Ustore->nnz); |
| 707 | |
| 708 | int* Lcol = m_l.outerIndexPtr(); |
| 709 | int* Lrow = m_l.innerIndexPtr(); |
| 710 | Scalar* Lval = m_l.valuePtr(); |
| 711 | |
| 712 | int* Ucol = m_u.outerIndexPtr(); |
| 713 | int* Urow = m_u.innerIndexPtr(); |
| 714 | Scalar* Uval = m_u.valuePtr(); |
| 715 | |
| 716 | Ucol[0] = 0; |
| 717 | Ucol[0] = 0; |
| 718 | |
| 719 | /* for each supernode */ |
| 720 | for (int k = 0; k <= Lstore->nsuper; ++k) |
| 721 | { |
| 722 | fsupc = L_FST_SUPC(k); |
| 723 | istart = L_SUB_START(fsupc); |
| 724 | nsupr = L_SUB_START(fsupc+1) - istart; |
| 725 | upper = 1; |
| 726 | |
| 727 | /* for each column in the supernode */ |
| 728 | for (int j = fsupc; j < L_FST_SUPC(k+1); ++j) |
| 729 | { |
| 730 | SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)]; |
| 731 | |
| 732 | /* Extract U */ |
| 733 | for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i) |
| 734 | { |
| 735 | Uval[lastu] = ((Scalar*)Ustore->nzval)[i]; |
| 736 | /* Matlab doesn't like explicit zero. */ |
| 737 | if (Uval[lastu] != 0.0) |
| 738 | Urow[lastu++] = U_SUB(i); |
| 739 | } |
| 740 | for (int i = 0; i < upper; ++i) |
| 741 | { |
| 742 | /* upper triangle in the supernode */ |
| 743 | Uval[lastu] = SNptr[i]; |
| 744 | /* Matlab doesn't like explicit zero. */ |
| 745 | if (Uval[lastu] != 0.0) |
| 746 | Urow[lastu++] = L_SUB(istart+i); |
| 747 | } |
| 748 | Ucol[j+1] = lastu; |
| 749 | |
| 750 | /* Extract L */ |
| 751 | Lval[lastl] = 1.0; /* unit diagonal */ |
| 752 | Lrow[lastl++] = L_SUB(istart + upper - 1); |
| 753 | for (int i = upper; i < nsupr; ++i) |
| 754 | { |
| 755 | Lval[lastl] = SNptr[i]; |
| 756 | /* Matlab doesn't like explicit zero. */ |
| 757 | if (Lval[lastl] != 0.0) |
| 758 | Lrow[lastl++] = L_SUB(istart+i); |
| 759 | } |
| 760 | Lcol[j+1] = lastl; |
| 761 | |
| 762 | ++upper; |
| 763 | } /* for j ... */ |
| 764 | |
| 765 | } /* for k ... */ |
| 766 | |
| 767 | // squeeze the matrices : |
| 768 | m_l.resizeNonZeros(lastl); |
| 769 | m_u.resizeNonZeros(lastu); |
| 770 | |
| 771 | m_extractedDataAreDirty = false; |
| 772 | } |
| 773 | } |
| 774 | |
| 775 | template<typename MatrixType> |
| 776 | typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const |
| 777 | { |
| 778 | 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()"); |
| 779 | |
| 780 | if (m_extractedDataAreDirty) |
| 781 | this->extractData(); |
| 782 | |
| 783 | Scalar det = Scalar(1); |
| 784 | for (int j=0; j<m_u.cols(); ++j) |
| 785 | { |
| 786 | if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0) |
| 787 | { |
| 788 | int lastId = m_u.outerIndexPtr()[j+1]-1; |
| 789 | eigen_assert(m_u.innerIndexPtr()[lastId]<=j); |
| 790 | if (m_u.innerIndexPtr()[lastId]==j) |
| 791 | det *= m_u.valuePtr()[lastId]; |
| 792 | } |
| 793 | } |
| 794 | if(m_sluEqued!='N') |
| 795 | return det/m_sluRscale.prod()/m_sluCscale.prod(); |
| 796 | else |
| 797 | return det; |
| 798 | } |
| 799 | |
| 800 | #ifdef EIGEN_PARSED_BY_DOXYGEN |
| 801 | #define EIGEN_SUPERLU_HAS_ILU |
| 802 | #endif |
| 803 | |
| 804 | #ifdef EIGEN_SUPERLU_HAS_ILU |
| 805 | |
| 806 | /** \ingroup SuperLUSupport_Module |
| 807 | * \class SuperILU |
| 808 | * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library |
| 809 | * |
| 810 | * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization |
| 811 | * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers. |
| 812 | * |
| 813 | * \warning This class requires SuperLU 4 or later. |
| 814 | * |
| 815 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| 816 | * |
| 817 | * \sa \ref TutorialSparseDirectSolvers, class ConjugateGradient, class BiCGSTAB |
| 818 | */ |
| 819 | |
| 820 | template<typename _MatrixType> |
| 821 | class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > |
| 822 | { |
| 823 | public: |
| 824 | typedef SuperLUBase<_MatrixType,SuperILU> Base; |
| 825 | typedef _MatrixType MatrixType; |
| 826 | typedef typename Base::Scalar Scalar; |
| 827 | typedef typename Base::RealScalar RealScalar; |
| 828 | typedef typename Base::Index Index; |
| 829 | |
| 830 | public: |
| 831 | |
| 832 | SuperILU() : Base() { init(); } |
| 833 | |
| 834 | SuperILU(const MatrixType& matrix) : Base() |
| 835 | { |
| 836 | init(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 837 | Base::compute(matrix); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 838 | } |
| 839 | |
| 840 | ~SuperILU() |
| 841 | { |
| 842 | } |
| 843 | |
| 844 | /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| 845 | * |
| 846 | * This function is particularly useful when solving for several problems having the same structure. |
| 847 | * |
| 848 | * \sa factorize() |
| 849 | */ |
| 850 | void analyzePattern(const MatrixType& matrix) |
| 851 | { |
| 852 | Base::analyzePattern(matrix); |
| 853 | } |
| 854 | |
| 855 | /** Performs a numeric decomposition of \a matrix |
| 856 | * |
| 857 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| 858 | * |
| 859 | * \sa analyzePattern() |
| 860 | */ |
| 861 | void factorize(const MatrixType& matrix); |
| 862 | |
| 863 | #ifndef EIGEN_PARSED_BY_DOXYGEN |
| 864 | /** \internal */ |
| 865 | template<typename Rhs,typename Dest> |
| 866 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; |
| 867 | #endif // EIGEN_PARSED_BY_DOXYGEN |
| 868 | |
| 869 | protected: |
| 870 | |
| 871 | using Base::m_matrix; |
| 872 | using Base::m_sluOptions; |
| 873 | using Base::m_sluA; |
| 874 | using Base::m_sluB; |
| 875 | using Base::m_sluX; |
| 876 | using Base::m_p; |
| 877 | using Base::m_q; |
| 878 | using Base::m_sluEtree; |
| 879 | using Base::m_sluEqued; |
| 880 | using Base::m_sluRscale; |
| 881 | using Base::m_sluCscale; |
| 882 | using Base::m_sluL; |
| 883 | using Base::m_sluU; |
| 884 | using Base::m_sluStat; |
| 885 | using Base::m_sluFerr; |
| 886 | using Base::m_sluBerr; |
| 887 | using Base::m_l; |
| 888 | using Base::m_u; |
| 889 | |
| 890 | using Base::m_analysisIsOk; |
| 891 | using Base::m_factorizationIsOk; |
| 892 | using Base::m_extractedDataAreDirty; |
| 893 | using Base::m_isInitialized; |
| 894 | using Base::m_info; |
| 895 | |
| 896 | void init() |
| 897 | { |
| 898 | Base::init(); |
| 899 | |
| 900 | ilu_set_default_options(&m_sluOptions); |
| 901 | m_sluOptions.PrintStat = NO; |
| 902 | m_sluOptions.ConditionNumber = NO; |
| 903 | m_sluOptions.Trans = NOTRANS; |
| 904 | m_sluOptions.ColPerm = MMD_AT_PLUS_A; |
| 905 | |
| 906 | // no attempt to preserve column sum |
| 907 | m_sluOptions.ILU_MILU = SILU; |
| 908 | // only basic ILU(k) support -- no direct control over memory consumption |
| 909 | // better to use ILU_DropRule = DROP_BASIC | DROP_AREA |
| 910 | // and set ILU_FillFactor to max memory growth |
| 911 | m_sluOptions.ILU_DropRule = DROP_BASIC; |
| 912 | m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10; |
| 913 | } |
| 914 | |
| 915 | private: |
| 916 | SuperILU(SuperILU& ) { } |
| 917 | }; |
| 918 | |
| 919 | template<typename MatrixType> |
| 920 | void SuperILU<MatrixType>::factorize(const MatrixType& a) |
| 921 | { |
| 922 | eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); |
| 923 | if(!m_analysisIsOk) |
| 924 | { |
| 925 | m_info = InvalidInput; |
| 926 | return; |
| 927 | } |
| 928 | |
| 929 | this->initFactorization(a); |
| 930 | |
| 931 | int info = 0; |
| 932 | RealScalar recip_pivot_growth, rcond; |
| 933 | |
| 934 | StatInit(&m_sluStat); |
| 935 | SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], |
| 936 | &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], |
| 937 | &m_sluL, &m_sluU, |
| 938 | NULL, 0, |
| 939 | &m_sluB, &m_sluX, |
| 940 | &recip_pivot_growth, &rcond, |
| 941 | &m_sluStat, &info, Scalar()); |
| 942 | StatFree(&m_sluStat); |
| 943 | |
| 944 | // FIXME how to better check for errors ??? |
| 945 | m_info = info == 0 ? Success : NumericalIssue; |
| 946 | m_factorizationIsOk = true; |
| 947 | } |
| 948 | |
| 949 | template<typename MatrixType> |
| 950 | template<typename Rhs,typename Dest> |
| 951 | void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const |
| 952 | { |
| 953 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); |
| 954 | |
| 955 | const int size = m_matrix.rows(); |
| 956 | const int rhsCols = b.cols(); |
| 957 | eigen_assert(size==b.rows()); |
| 958 | |
| 959 | m_sluOptions.Trans = NOTRANS; |
| 960 | m_sluOptions.Fact = FACTORED; |
| 961 | m_sluOptions.IterRefine = NOREFINE; |
| 962 | |
| 963 | m_sluFerr.resize(rhsCols); |
| 964 | m_sluBerr.resize(rhsCols); |
| 965 | m_sluB = SluMatrix::Map(b.const_cast_derived()); |
| 966 | m_sluX = SluMatrix::Map(x.derived()); |
| 967 | |
| 968 | typename Rhs::PlainObject b_cpy; |
| 969 | if(m_sluEqued!='N') |
| 970 | { |
| 971 | b_cpy = b; |
| 972 | m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); |
| 973 | } |
| 974 | |
| 975 | int info = 0; |
| 976 | RealScalar recip_pivot_growth, rcond; |
| 977 | |
| 978 | StatInit(&m_sluStat); |
| 979 | SuperLU_gsisx(&m_sluOptions, &m_sluA, |
| 980 | m_q.data(), m_p.data(), |
| 981 | &m_sluEtree[0], &m_sluEqued, |
| 982 | &m_sluRscale[0], &m_sluCscale[0], |
| 983 | &m_sluL, &m_sluU, |
| 984 | NULL, 0, |
| 985 | &m_sluB, &m_sluX, |
| 986 | &recip_pivot_growth, &rcond, |
| 987 | &m_sluStat, &info, Scalar()); |
| 988 | StatFree(&m_sluStat); |
| 989 | |
| 990 | m_info = info==0 ? Success : NumericalIssue; |
| 991 | } |
| 992 | #endif |
| 993 | |
| 994 | namespace internal { |
| 995 | |
| 996 | template<typename _MatrixType, typename Derived, typename Rhs> |
| 997 | struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> |
| 998 | : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs> |
| 999 | { |
| 1000 | typedef SuperLUBase<_MatrixType,Derived> Dec; |
| 1001 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) |
| 1002 | |
| 1003 | template<typename Dest> void evalTo(Dest& dst) const |
| 1004 | { |
| 1005 | dec().derived()._solve(rhs(),dst); |
| 1006 | } |
| 1007 | }; |
| 1008 | |
| 1009 | template<typename _MatrixType, typename Derived, typename Rhs> |
| 1010 | struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> |
| 1011 | : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs> |
| 1012 | { |
| 1013 | typedef SuperLUBase<_MatrixType,Derived> Dec; |
| 1014 | EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) |
| 1015 | |
| 1016 | template<typename Dest> void evalTo(Dest& dst) const |
| 1017 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 1018 | this->defaultEvalTo(dst); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 1019 | } |
| 1020 | }; |
| 1021 | |
| 1022 | } // end namespace internal |
| 1023 | |
| 1024 | } // end namespace Eigen |
| 1025 | |
| 1026 | #endif // EIGEN_SUPERLUSUPPORT_H |