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 Derived& derived() 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *static_cast<Derived*>(this); 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Derived& derived() const 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *static_cast<const Derived*>(this); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns a reference to the integer vector IPARM of PaStiX parameters 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to modify the default parameters. 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The statistics related to the different phases of factorization and solve are saved here as well 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Array<Index,IPARM_SIZE,1>& iparm() 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm; 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Return a reference to a particular index parameter of the IPARM vector 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa iparm() 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int& iparm(int idxparam) 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm(idxparam); 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns a reference to the double vector DPARM of PaStiX parameters 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The statistics related to the different phases of factorization and solve are saved here as well 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Array<RealScalar,IPARM_SIZE,1>& dparm() 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_dparm; 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Return a reference to a particular index parameter of the DPARM vector 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa dparm() 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double& dparm(int idxparam) 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_dparm(idxparam); 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_size; } 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_size; } 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the PaStiX reports a problem 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c InvalidInput if the input matrix is invalid 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa iparm() 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Decomposition is not initialized."); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::sparse_solve_retval<PastixBase, Rhs> 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const SparseMatrixBase<Rhs>& b) const 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized."); 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows()==b.rows() 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived()); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Initialize the Pastix data structure, check the matrix 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init(); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the ordering and the symbolic factorization 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(ColSpMatrix& mat); 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the numerical factorization 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(ColSpMatrix& mat); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Free all the data allocated by Pastix 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void clean() 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_CLEAN; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute(ColSpMatrix& mat); 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_initisOk; 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_analysisIsOk; 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int m_factorizationIsOk; 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable ComputationInfo m_info; 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable pastix_data_t *m_pastixdata; // Data structure for pastix 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable int m_comm; // The MPI communicator identifier 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mutable int m_size; // Size of the matrix 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initialize the PaStiX data structure. 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *A first call to this function fills iparm and dparm with the default PaStiX parameters 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa iparm() dparm() 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::init() 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = 0; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm.setZero(IPARM_SIZE); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_dparm.setZero(DPARM_SIZE); 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_MODIFY_PARAMETER) = API_NO; 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pastix(&m_pastixdata, MPI_COMM_WORLD, 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0, 0, 0, 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_VERBOSE] = 2; 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_INCOMPLETE] = API_NO; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_OOC_LIMIT] = 2000; 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_RHS_MAKING] = API_RHS_B; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_INIT; 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_INIT; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0, 0, 0, m_iparm.data(), m_dparm.data()); 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm(IPARM_ERROR_NUMBER)) { 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = InvalidInput; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_initisOk = false; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else { 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_initisOk = true; 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::compute(ColSpMatrix& mat) 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared"); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath analyzePattern(mat); 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath factorize(mat); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = m_factorizationIsOk; 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_initisOk && "The initialization of PaSTiX failed"); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // clean previous calls 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_size>0) 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath clean(); 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = mat.rows(); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.resize(m_size); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_invp.resize(m_size); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_ORDERING; 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm(IPARM_ERROR_NUMBER)) 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = false; 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_analysisIsOk = true; 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <class Derived> 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid PastixBase<Derived>::factorize(ColSpMatrix& mat) 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// if(&m_cpyMat != &mat) m_cpyMat = mat; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase"); 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT; 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT; 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_size = mat.rows(); 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_iparm(IPARM_ERROR_NUMBER)) 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = false; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = false; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_factorizationIsOk = true; 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* Solve the system */ 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Base> 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Rhs,typename Dest> 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "The matrix should be factorized first"); 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int rhs = 1; 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = b; /* on return, x is overwritten by the computed solution */ 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i < b.cols(); i++){ 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_START_TASK] = API_TASK_SOLVE; 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm[IPARM_END_TASK] = API_TASK_REFINE; 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0, 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the returned error 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_iparm(IPARM_ERROR_NUMBER)==0; 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PastixLU 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Sparse direct LU solver based on PaStiX library 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B with a supernodal LU 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * factorization in the PaStiX library. The matrix A should be squared and nonsingular 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * PaStiX requires that the matrix A has a symmetric structural pattern. 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This interface can symmetrize the input matrix otherwise. 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse. 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam IsStrSym Indicates if the input matrix has a symmetric pattern, default is false 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * NOTE : Note that if the analysis and factorization phase are called separately, 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the input matrix will be symmetrized at each call, hence it is advised to 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * symmetrize the matrix in a end-user program and set \p IsStrSym to true 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa \ref TutorialSparseDirectSolvers 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, bool IsStrSym> 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixLU : public PastixBase< PastixLU<_MatrixType> > 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<PastixLU<MatrixType> > Base; 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::ColSpMatrix ColSpMatrix; 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLU() : Base() 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLU(const MatrixType& matrix):Base() 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LU supernodal factorization of \p matrix. 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * iparm and dparm can be used to tune the PaStiX parameters. 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * see the PaStiX user's manual 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute (const MatrixType& matrix) 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = false; 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(temp); 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Several ordering methods can be used at this step. See the PaStiX user's manual. 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result of this operation can be used with successive matrices having the same pattern as \p matrix 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& matrix) 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = false; 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(temp); 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LU supernodal factorization of \p matrix 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING The matrix \p matrix should have the same structural pattern 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * as the same used in the analysis phase. 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& matrix) 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::factorize(temp); 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = false; 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_SYM) = API_SYM_NO; 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(IsStrSym) 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out = matrix; 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!m_structureIsUptodate) 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // update the transposed structure 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transposedStructure = matrix.transpose(); 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Set the elements of the matrix to zero 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it) 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath it.valueRef() = 0.0; 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_structureIsUptodate = true; 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out = m_transposedStructure + matrix; 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::c_to_fortran_numbering(out); 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iparm; 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_dparm; 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix m_transposedStructure; 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_structureIsUptodate; 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PastixLLT 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B via a LL^T supernodal Cholesky factorization 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * available in the PaStiX library. The matrix A should be symmetric and positive definite 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa \ref TutorialSparseDirectSolvers 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base; 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::ColSpMatrix ColSpMatrix; 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLLT() : Base() 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLLT(const MatrixType& matrix):Base() 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the L factor of the LL^T supernodal factorization of \p matrix 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute (const MatrixType& matrix) 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(temp); 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LL^T symbolic factorization of \p matrix using its sparsity pattern 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result of this operation can be used with successive matrices having the same pattern as \p matrix 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& matrix) 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(temp); 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LL^T supernodal numerical factorization of \p matrix 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& matrix) 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::factorize(temp); 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iparm; 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_SYM) = API_SYM_YES; 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Pastix supports only lower, column-major matrices 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::c_to_fortran_numbering(out); 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup PaStiXSupport_Module 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PastixLDLT 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class is used to solve the linear systems A.X = B via a LDL^T supernodal Cholesky factorization 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * available in the PaStiX library. The matrix A should be symmetric and positive definite 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The vectors or matrices X and B can be either dense or sparse 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa \ref TutorialSparseDirectSolvers 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Base::ColSpMatrix ColSpMatrix; 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { UpLo = _UpLo }; 631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLDLT():Base() 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PastixLDLT(const MatrixType& matrix):Base() 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath init(); 639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the L and D factors of the LDL^T factorization of \p matrix 643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa analyzePattern() factorize() 644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute (const MatrixType& matrix) 646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(temp); 650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LDL^T symbolic factorization of \p matrix using its sparsity pattern 653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result of this operation can be used with successive matrices having the same pattern as \p matrix 654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa factorize() 655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void analyzePattern(const MatrixType& matrix) 657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::analyzePattern(temp); 661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Compute the LDL^T supernodal numerical factorization of \p matrix 663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void factorize(const MatrixType& matrix) 666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColSpMatrix temp; 668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath grabMatrix(matrix, temp); 669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::factorize(temp); 670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iparm; 674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void init() 676c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 677c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_SYM) = API_SYM_YES; 678c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; 679c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 680c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 681c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 682c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 683c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Pastix supports only lower, column-major matrices 684c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::c_to_fortran_numbering(out); 686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 690c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 691c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 692c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<PastixBase<_MatrixType>, Rhs> 693c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<PastixBase<_MatrixType>, Rhs> 694c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 695c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<_MatrixType> Dec; 696c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 697c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 698c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 699c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 700c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve(rhs(),dst); 701c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 702c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 703c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 704c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 705c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> 706c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs> 707c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 708c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PastixBase<_MatrixType> Dec; 709c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 710c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 711c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 712c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 7137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->defaultEvalTo(dst); 714c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 715c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 716c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 717c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 718c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 719c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 720c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 721c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 722