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