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 { SluMatrixSluMatrix92 SluMatrix() 93 { 94 Store = &storage; 95 } 96 SluMatrixSluMatrix97 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 setStorageTypeSluMatrix120 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> setScalarTypeSluMatrix133 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> MapSluMatrix150 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(); 163 res.storage.values = mat.data(); 164 return res; 165 } 166 167 template<typename MatrixType> MapSluMatrix168 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 */ 356 // template<typename Rhs> 357 // inline const internal::sparse_solve_retval<SuperLU, 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<SuperLU, Rhs>(*this, b.derived()); 363 // } 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> 380 void dumpMemory(Stream& s) 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 { 499 Base::init(); 500 compute(matrix); 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 615 int info = 0; 616 RealScalar recip_pivot_growth, rcond; 617 RealScalar ferr, berr; 618 619 StatInit(&m_sluStat); 620 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], 621 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], 622 &m_sluL, &m_sluU, 623 NULL, 0, 624 &m_sluB, &m_sluX, 625 &recip_pivot_growth, &rcond, 626 &ferr, &berr, 627 &m_sluStat, &info, Scalar()); 628 StatFree(&m_sluStat); 629 630 m_extractedDataAreDirty = true; 631 632 // FIXME how to better check for errors ??? 633 m_info = info == 0 ? Success : NumericalIssue; 634 m_factorizationIsOk = true; 635 } 636 637 template<typename MatrixType> 638 template<typename Rhs,typename Dest> 639 void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const 640 { 641 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); 642 643 const int size = m_matrix.rows(); 644 const int rhsCols = b.cols(); 645 eigen_assert(size==b.rows()); 646 647 m_sluOptions.Trans = NOTRANS; 648 m_sluOptions.Fact = FACTORED; 649 m_sluOptions.IterRefine = NOREFINE; 650 651 652 m_sluFerr.resize(rhsCols); 653 m_sluBerr.resize(rhsCols); 654 m_sluB = SluMatrix::Map(b.const_cast_derived()); 655 m_sluX = SluMatrix::Map(x.derived()); 656 657 typename Rhs::PlainObject b_cpy; 658 if(m_sluEqued!='N') 659 { 660 b_cpy = b; 661 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); 662 } 663 664 StatInit(&m_sluStat); 665 int info = 0; 666 RealScalar recip_pivot_growth, rcond; 667 SuperLU_gssvx(&m_sluOptions, &m_sluA, 668 m_q.data(), m_p.data(), 669 &m_sluEtree[0], &m_sluEqued, 670 &m_sluRscale[0], &m_sluCscale[0], 671 &m_sluL, &m_sluU, 672 NULL, 0, 673 &m_sluB, &m_sluX, 674 &recip_pivot_growth, &rcond, 675 &m_sluFerr[0], &m_sluBerr[0], 676 &m_sluStat, &info, Scalar()); 677 StatFree(&m_sluStat); 678 m_info = info==0 ? Success : NumericalIssue; 679 } 680 681 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code, 682 // 683 // Copyright (c) 1994 by Xerox Corporation. All rights reserved. 684 // 685 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 686 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 687 // 688 template<typename MatrixType, typename Derived> 689 void SuperLUBase<MatrixType,Derived>::extractData() const 690 { 691 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()"); 692 if (m_extractedDataAreDirty) 693 { 694 int upper; 695 int fsupc, istart, nsupr; 696 int lastl = 0, lastu = 0; 697 SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store); 698 NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store); 699 Scalar *SNptr; 700 701 const int size = m_matrix.rows(); 702 m_l.resize(size,size); 703 m_l.resizeNonZeros(Lstore->nnz); 704 m_u.resize(size,size); 705 m_u.resizeNonZeros(Ustore->nnz); 706 707 int* Lcol = m_l.outerIndexPtr(); 708 int* Lrow = m_l.innerIndexPtr(); 709 Scalar* Lval = m_l.valuePtr(); 710 711 int* Ucol = m_u.outerIndexPtr(); 712 int* Urow = m_u.innerIndexPtr(); 713 Scalar* Uval = m_u.valuePtr(); 714 715 Ucol[0] = 0; 716 Ucol[0] = 0; 717 718 /* for each supernode */ 719 for (int k = 0; k <= Lstore->nsuper; ++k) 720 { 721 fsupc = L_FST_SUPC(k); 722 istart = L_SUB_START(fsupc); 723 nsupr = L_SUB_START(fsupc+1) - istart; 724 upper = 1; 725 726 /* for each column in the supernode */ 727 for (int j = fsupc; j < L_FST_SUPC(k+1); ++j) 728 { 729 SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)]; 730 731 /* Extract U */ 732 for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i) 733 { 734 Uval[lastu] = ((Scalar*)Ustore->nzval)[i]; 735 /* Matlab doesn't like explicit zero. */ 736 if (Uval[lastu] != 0.0) 737 Urow[lastu++] = U_SUB(i); 738 } 739 for (int i = 0; i < upper; ++i) 740 { 741 /* upper triangle in the supernode */ 742 Uval[lastu] = SNptr[i]; 743 /* Matlab doesn't like explicit zero. */ 744 if (Uval[lastu] != 0.0) 745 Urow[lastu++] = L_SUB(istart+i); 746 } 747 Ucol[j+1] = lastu; 748 749 /* Extract L */ 750 Lval[lastl] = 1.0; /* unit diagonal */ 751 Lrow[lastl++] = L_SUB(istart + upper - 1); 752 for (int i = upper; i < nsupr; ++i) 753 { 754 Lval[lastl] = SNptr[i]; 755 /* Matlab doesn't like explicit zero. */ 756 if (Lval[lastl] != 0.0) 757 Lrow[lastl++] = L_SUB(istart+i); 758 } 759 Lcol[j+1] = lastl; 760 761 ++upper; 762 } /* for j ... */ 763 764 } /* for k ... */ 765 766 // squeeze the matrices : 767 m_l.resizeNonZeros(lastl); 768 m_u.resizeNonZeros(lastu); 769 770 m_extractedDataAreDirty = false; 771 } 772 } 773 774 template<typename MatrixType> 775 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const 776 { 777 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()"); 778 779 if (m_extractedDataAreDirty) 780 this->extractData(); 781 782 Scalar det = Scalar(1); 783 for (int j=0; j<m_u.cols(); ++j) 784 { 785 if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0) 786 { 787 int lastId = m_u.outerIndexPtr()[j+1]-1; 788 eigen_assert(m_u.innerIndexPtr()[lastId]<=j); 789 if (m_u.innerIndexPtr()[lastId]==j) 790 det *= m_u.valuePtr()[lastId]; 791 } 792 } 793 if(m_sluEqued!='N') 794 return det/m_sluRscale.prod()/m_sluCscale.prod(); 795 else 796 return det; 797 } 798 799 #ifdef EIGEN_PARSED_BY_DOXYGEN 800 #define EIGEN_SUPERLU_HAS_ILU 801 #endif 802 803 #ifdef EIGEN_SUPERLU_HAS_ILU 804 805 /** \ingroup SuperLUSupport_Module 806 * \class SuperILU 807 * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library 808 * 809 * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization 810 * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers. 811 * 812 * \warning This class requires SuperLU 4 or later. 813 * 814 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 815 * 816 * \sa \ref TutorialSparseDirectSolvers, class ConjugateGradient, class BiCGSTAB 817 */ 818 819 template<typename _MatrixType> 820 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > 821 { 822 public: 823 typedef SuperLUBase<_MatrixType,SuperILU> Base; 824 typedef _MatrixType MatrixType; 825 typedef typename Base::Scalar Scalar; 826 typedef typename Base::RealScalar RealScalar; 827 typedef typename Base::Index Index; 828 829 public: 830 831 SuperILU() : Base() { init(); } 832 833 SuperILU(const MatrixType& matrix) : Base() 834 { 835 init(); 836 compute(matrix); 837 } 838 839 ~SuperILU() 840 { 841 } 842 843 /** Performs a symbolic decomposition on the sparcity of \a matrix. 844 * 845 * This function is particularly useful when solving for several problems having the same structure. 846 * 847 * \sa factorize() 848 */ 849 void analyzePattern(const MatrixType& matrix) 850 { 851 Base::analyzePattern(matrix); 852 } 853 854 /** Performs a numeric decomposition of \a matrix 855 * 856 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 857 * 858 * \sa analyzePattern() 859 */ 860 void factorize(const MatrixType& matrix); 861 862 #ifndef EIGEN_PARSED_BY_DOXYGEN 863 /** \internal */ 864 template<typename Rhs,typename Dest> 865 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; 866 #endif // EIGEN_PARSED_BY_DOXYGEN 867 868 protected: 869 870 using Base::m_matrix; 871 using Base::m_sluOptions; 872 using Base::m_sluA; 873 using Base::m_sluB; 874 using Base::m_sluX; 875 using Base::m_p; 876 using Base::m_q; 877 using Base::m_sluEtree; 878 using Base::m_sluEqued; 879 using Base::m_sluRscale; 880 using Base::m_sluCscale; 881 using Base::m_sluL; 882 using Base::m_sluU; 883 using Base::m_sluStat; 884 using Base::m_sluFerr; 885 using Base::m_sluBerr; 886 using Base::m_l; 887 using Base::m_u; 888 889 using Base::m_analysisIsOk; 890 using Base::m_factorizationIsOk; 891 using Base::m_extractedDataAreDirty; 892 using Base::m_isInitialized; 893 using Base::m_info; 894 895 void init() 896 { 897 Base::init(); 898 899 ilu_set_default_options(&m_sluOptions); 900 m_sluOptions.PrintStat = NO; 901 m_sluOptions.ConditionNumber = NO; 902 m_sluOptions.Trans = NOTRANS; 903 m_sluOptions.ColPerm = MMD_AT_PLUS_A; 904 905 // no attempt to preserve column sum 906 m_sluOptions.ILU_MILU = SILU; 907 // only basic ILU(k) support -- no direct control over memory consumption 908 // better to use ILU_DropRule = DROP_BASIC | DROP_AREA 909 // and set ILU_FillFactor to max memory growth 910 m_sluOptions.ILU_DropRule = DROP_BASIC; 911 m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10; 912 } 913 914 private: 915 SuperILU(SuperILU& ) { } 916 }; 917 918 template<typename MatrixType> 919 void SuperILU<MatrixType>::factorize(const MatrixType& a) 920 { 921 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 922 if(!m_analysisIsOk) 923 { 924 m_info = InvalidInput; 925 return; 926 } 927 928 this->initFactorization(a); 929 930 int info = 0; 931 RealScalar recip_pivot_growth, rcond; 932 933 StatInit(&m_sluStat); 934 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], 935 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], 936 &m_sluL, &m_sluU, 937 NULL, 0, 938 &m_sluB, &m_sluX, 939 &recip_pivot_growth, &rcond, 940 &m_sluStat, &info, Scalar()); 941 StatFree(&m_sluStat); 942 943 // FIXME how to better check for errors ??? 944 m_info = info == 0 ? Success : NumericalIssue; 945 m_factorizationIsOk = true; 946 } 947 948 template<typename MatrixType> 949 template<typename Rhs,typename Dest> 950 void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const 951 { 952 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); 953 954 const int size = m_matrix.rows(); 955 const int rhsCols = b.cols(); 956 eigen_assert(size==b.rows()); 957 958 m_sluOptions.Trans = NOTRANS; 959 m_sluOptions.Fact = FACTORED; 960 m_sluOptions.IterRefine = NOREFINE; 961 962 m_sluFerr.resize(rhsCols); 963 m_sluBerr.resize(rhsCols); 964 m_sluB = SluMatrix::Map(b.const_cast_derived()); 965 m_sluX = SluMatrix::Map(x.derived()); 966 967 typename Rhs::PlainObject b_cpy; 968 if(m_sluEqued!='N') 969 { 970 b_cpy = b; 971 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); 972 } 973 974 int info = 0; 975 RealScalar recip_pivot_growth, rcond; 976 977 StatInit(&m_sluStat); 978 SuperLU_gsisx(&m_sluOptions, &m_sluA, 979 m_q.data(), m_p.data(), 980 &m_sluEtree[0], &m_sluEqued, 981 &m_sluRscale[0], &m_sluCscale[0], 982 &m_sluL, &m_sluU, 983 NULL, 0, 984 &m_sluB, &m_sluX, 985 &recip_pivot_growth, &rcond, 986 &m_sluStat, &info, Scalar()); 987 StatFree(&m_sluStat); 988 989 m_info = info==0 ? Success : NumericalIssue; 990 } 991 #endif 992 993 namespace internal { 994 995 template<typename _MatrixType, typename Derived, typename Rhs> 996 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> 997 : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs> 998 { 999 typedef SuperLUBase<_MatrixType,Derived> Dec; 1000 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 1001 1002 template<typename Dest> void evalTo(Dest& dst) const 1003 { 1004 dec().derived()._solve(rhs(),dst); 1005 } 1006 }; 1007 1008 template<typename _MatrixType, typename Derived, typename Rhs> 1009 struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> 1010 : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs> 1011 { 1012 typedef SuperLUBase<_MatrixType,Derived> Dec; 1013 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 1014 1015 template<typename Dest> void evalTo(Dest& dst) const 1016 { 1017 dec().derived()._solve(rhs(),dst); 1018 } 1019 }; 1020 1021 } // end namespace internal 1022 1023 } // end namespace Eigen 1024 1025 #endif // EIGEN_SUPERLUSUPPORT_H 1026