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 = (void*)(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<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 } 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 init(); 500 Base::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 m_sluOptions.ColPerm = COLAMD; 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(); 837 Base::compute(matrix); 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 { 1018 this->defaultEvalTo(dst); 1019 } 1020 }; 1021 1022 } // end namespace internal 1023 1024 } // end namespace Eigen 1025 1026 #endif // EIGEN_SUPERLUSUPPORT_H 1027