1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_PASTIXSUPPORT_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_PASTIXSUPPORT_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Interface to the PaStix solver 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B via the PaStix library. 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix can be either real or complex, symmetric or not. 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa TutorialSparseDirectSolvers 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, bool IsStrSym = false> class PastixLU; 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int Options> class PastixLLT; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int Options> class PastixLDLT; 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<class Pastix> struct pastix_traits; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType> 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pastix_traits< PastixLU<_MatrixType> > 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Scalar Scalar; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::RealScalar RealScalar; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Index Index; 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType, int Options> 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pastix_traits< PastixLLT<_MatrixType,Options> > 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Scalar Scalar; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::RealScalar RealScalar; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Index Index; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType, int Options> 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath struct pastix_traits< PastixLDLT<_MatrixType,Options> > 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Scalar Scalar; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::RealScalar RealScalar; 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename _MatrixType::Index Index; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, float *vals, int *perm, int * invp, float *x, int nbrhs, int *iparm, double *dparm) 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nbrhs == 0) {x = NULL; nbrhs=1;} 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, double *vals, int *perm, int * invp, double *x, int nbrhs, int *iparm, double *dparm) 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nbrhs == 0) {x = NULL; nbrhs=1;} 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<float> *vals, int *perm, int * invp, std::complex<float> *x, int nbrhs, int *iparm, double *dparm) 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nbrhs == 0) {x = NULL; nbrhs=1;} 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm) 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nbrhs == 0) {x = NULL; nbrhs=1;} 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Convert the matrix to Fortran-style Numbering 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename MatrixType> 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void c_to_fortran_numbering (MatrixType& mat) 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( !(mat.outerIndexPtr()[0]) ) 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(i = 0; i <= mat.rows(); ++i) 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++mat.outerIndexPtr()[i]; 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(i = 0; i < mat.nonZeros(); ++i) 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++mat.innerIndexPtr()[i]; 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Convert to C-style Numbering 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename MatrixType> 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void fortran_to_c_numbering (MatrixType& mat) 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the Numbering 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ( mat.outerIndexPtr()[0] == 1 ) 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { // Convert to C-style numbering 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(i = 0; i <= mat.rows(); ++i) 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --mat.outerIndexPtr()[i]; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(i = 0; i < mat.nonZeros(); ++i) 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --mat.innerIndexPtr()[i]; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This is the base class to interface with PaStiX functions. 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Users should not used this class directly. 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixBase : internal::noncopyable 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType; 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> Vector; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix; 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0) 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ~PastixBase() 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath clean(); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::solve_retval<PastixBase, Rhs> 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Pastix solver is not initialized."); 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==b.rows() 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<PastixBase, Rhs>(*this, b.derived()); 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==b.rows()); 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const int NbColsAtOnce = 1; 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int rhsCols = b.cols(); 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int size = b.rows(); 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int k=0; k<rhsCols; k+=NbColsAtOnce) 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.leftCols(actualCols) = b.middleCols(k,actualCols); 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols)); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView(); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Derived& derived() 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *static_cast<Derived*>(this); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Derived& derived() const 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *static_cast<const Derived*>(this); 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns a reference to the integer vector IPARM of PaStiX parameters 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to modify the default parameters. 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The statistics related to the different phases of factorization and solve are saved here as well 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Array<Index,IPARM_SIZE,1>& iparm() 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Return a reference to a particular index parameter of the IPARM vector 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa iparm() 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int& iparm(int idxparam) 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm(idxparam); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns a reference to the double vector DPARM of PaStiX parameters 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The statistics related to the different phases of factorization and solve are saved here as well 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Array<RealScalar,IPARM_SIZE,1>& dparm() 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_dparm; 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Return a reference to a particular index parameter of the DPARM vector 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa dparm() 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double& dparm(int idxparam) 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_dparm(idxparam); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_size; } 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_size; } 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the PaStiX reports a problem 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c InvalidInput if the input matrix is invalid 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa iparm() 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Decomposition is not initialized."); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::sparse_solve_retval<PastixBase, Rhs> 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const SparseMatrixBase<Rhs>& b) const 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized."); 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==b.rows() 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived()); 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Initialize the Pastix data structure, check the matrix 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init(); 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the ordering and the symbolic factorization 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(ColSpMatrix& mat); 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the numerical factorization 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(ColSpMatrix& mat); 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Free all the data allocated by Pastix 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void clean() 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_CLEAN; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute(ColSpMatrix& mat); 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_initisOk; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_analysisIsOk; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_factorizationIsOk; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable ComputationInfo m_info; 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable pastix_data_t *m_pastixdata; // Data structure for pastix 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable int m_comm; // The MPI communicator identifier 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable int m_size; // Size of the matrix 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initialize the PaStiX data structure. 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *A first call to this function fills iparm and dparm with the default PaStiX parameters 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa iparm() dparm() 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::init() 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = 0; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm.setZero(IPARM_SIZE); 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_dparm.setZero(DPARM_SIZE); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_MODIFY_PARAMETER) = API_NO; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pastix(&m_pastixdata, MPI_COMM_WORLD, 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0, 0, 0, 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_VERBOSE] = 2; 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_INCOMPLETE] = API_NO; 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_OOC_LIMIT] = 2000; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_RHS_MAKING] = API_RHS_B; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_INIT; 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_INIT; 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0, 0, 0, m_iparm.data(), m_dparm.data()); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm(IPARM_ERROR_NUMBER)) { 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = InvalidInput; 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_initisOk = false; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else { 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_initisOk = true; 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::compute(ColSpMatrix& mat) 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared"); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath analyzePattern(mat); 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath factorize(mat); 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = m_factorizationIsOk; 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_initisOk && "The initialization of PaSTiX failed"); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // clean previous calls 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_size>0) 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath clean(); 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = mat.rows(); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.resize(m_size); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_invp.resize(m_size); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_ORDERING; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE; 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm(IPARM_ERROR_NUMBER)) 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = false; 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::factorize(ColSpMatrix& mat) 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// if(&m_cpyMat != &mat) m_cpyMat = mat; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase"); 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT; 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT; 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = mat.rows(); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm(IPARM_ERROR_NUMBER)) 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = false; 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = false; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* Solve the system */ 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Base> 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Rhs,typename Dest> 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "The matrix should be factorized first"); 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int rhs = 1; 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = b; /* on return, x is overwritten by the computed solution */ 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i < b.cols(); i++){ 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_START_TASK] = API_TASK_SOLVE; 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_END_TASK] = API_TASK_REFINE; 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0, 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue; 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm(IPARM_ERROR_NUMBER)==0; 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PastixLU 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Sparse direct LU solver based on PaStiX library 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B with a supernodal LU 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * factorization in the PaStiX library. The matrix A should be squared and nonsingular 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * PaStiX requires that the matrix A has a symmetric structural pattern. 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This interface can symmetrize the input matrix otherwise. 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse. 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam IsStrSym Indicates if the input matrix has a symmetric pattern, default is false 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * NOTE : Note that if the analysis and factorization phase are called separately, 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the input matrix will be symmetrized at each call, hence it is advised to 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * symmetrize the matrix in a end-user program and set \p IsStrSym to true 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa \ref TutorialSparseDirectSolvers 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, bool IsStrSym> 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixLU : public PastixBase< PastixLU<_MatrixType> > 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<PastixLU<MatrixType> > Base; 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::ColSpMatrix ColSpMatrix; 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLU() : Base() 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLU(const MatrixType& matrix):Base() 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LU supernodal factorization of \p matrix. 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * iparm and dparm can be used to tune the PaStiX parameters. 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * see the PaStiX user's manual 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute (const MatrixType& matrix) 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = false; 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(temp); 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Several ordering methods can be used at this step. See the PaStiX user's manual. 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result of this operation can be used with successive matrices having the same pattern as \p matrix 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& matrix) 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = false; 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(temp); 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LU supernodal factorization of \p matrix 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING The matrix \p matrix should have the same structural pattern 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * as the same used in the analysis phase. 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& matrix) 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::factorize(temp); 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = false; 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_SYM) = API_SYM_NO; 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(IsStrSym) 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out = matrix; 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!m_structureIsUptodate) 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // update the transposed structure 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transposedStructure = matrix.transpose(); 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Set the elements of the matrix to zero 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it) 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath it.valueRef() = 0.0; 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = true; 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out = m_transposedStructure + matrix; 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::c_to_fortran_numbering(out); 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iparm; 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_dparm; 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix m_transposedStructure; 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_structureIsUptodate; 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PastixLLT 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B via a LL^T supernodal Cholesky factorization 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * available in the PaStiX library. The matrix A should be symmetric and positive definite 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa \ref TutorialSparseDirectSolvers 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base; 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::ColSpMatrix ColSpMatrix; 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLLT() : Base() 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLLT(const MatrixType& matrix):Base() 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the L factor of the LL^T supernodal factorization of \p matrix 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute (const MatrixType& matrix) 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(temp); 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LL^T symbolic factorization of \p matrix using its sparsity pattern 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result of this operation can be used with successive matrices having the same pattern as \p matrix 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& matrix) 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(temp); 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LL^T supernodal numerical factorization of \p matrix 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& matrix) 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::factorize(temp); 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iparm; 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_SYM) = API_SYM_YES; 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Pastix supports only lower, column-major matrices 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::c_to_fortran_numbering(out); 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PastixLDLT 630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B via a LDL^T supernodal Cholesky factorization 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * available in the PaStiX library. The matrix A should be symmetric and positive definite 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa \ref TutorialSparseDirectSolvers 641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> 643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > 644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::ColSpMatrix ColSpMatrix; 649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLDLT():Base() 653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLDLT(const MatrixType& matrix):Base() 658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the L and D factors of the LDL^T factorization of \p matrix 664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute (const MatrixType& matrix) 667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(temp); 671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LDL^T symbolic factorization of \p matrix using its sparsity pattern 674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result of this operation can be used with successive matrices having the same pattern as \p matrix 675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 676c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 677c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& matrix) 678c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 679c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 680c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 681c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(temp); 682c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 683c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LDL^T supernodal numerical factorization of \p matrix 684c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& matrix) 687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 690c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::factorize(temp); 691c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 692c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 693c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 694c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iparm; 695c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 696c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 697c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 698c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_SYM) = API_SYM_YES; 699c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; 700c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 701c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 702c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 703c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 704c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Pastix supports only lower, column-major matrices 705c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 706c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::c_to_fortran_numbering(out); 707c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 708c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 709c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 710c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 711c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 712c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 713c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<PastixBase<_MatrixType>, Rhs> 714c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<PastixBase<_MatrixType>, Rhs> 715c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 716c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<_MatrixType> Dec; 717c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 718c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 719c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 720c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 721c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve(rhs(),dst); 722c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 723c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 724c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 725c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 726c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> 727c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs> 728c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 729c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<_MatrixType> Dec; 730c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 731c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 732c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 733c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 734c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve_sparse(rhs(),dst); 735c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 736c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 737c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 738c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 739c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 740c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 741c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 742c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 743