17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_SPARSE_LU_H 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_SPARSE_LU_H 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::Index> > class SparseLU; 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType; 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType; 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \ingroup SparseLU_Module 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class SparseLU 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Sparse supernodal LU factorization for general matrices 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This class implements the supernodal LU factorization for general matrices. 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It uses the main techniques from the sequential SuperLU package 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles transparently real 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and complex arithmetics with single and double precision, depending on the 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * scalar type of your input matrix. 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The code has been optimized to provide BLAS-3 operations during supernode-panel updates. 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It benefits directly from the built-in high-performant Eigen BLAS routines. 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Moreover, when the size of a supernode is very small, the BLAS calls are avoided to 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * enable a better optimization from the compiler. For best performance, 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors. 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * An important parameter of this class is the ordering method. It is used to reorder the columns 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * (and eventually the rows) of the matrix to reduce the number of new elements that are created during 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * numerical factorization. The cheapest method available is COLAMD. 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * See \link OrderingMethods_Module the OrderingMethods module \endlink for the list of 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * built-in and external ordering methods. 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Simple example with key steps 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \code 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * VectorXd x(n), b(n); 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * SparseMatrix<double, ColMajor> A; 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<Index> > solver; 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * // fill A and b; 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * // Compute the ordering permutation vector from the structural pattern of A 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * solver.analyzePattern(A); 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * // Compute the numerical factorization 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * solver.factorize(A); 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * //Use the factors to solve the linear system 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * x = solver.solve(b); 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \endcode 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \warning The input matrix A should be in a \b compressed and \b column-major form. 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix. 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \note Unlike the initial SuperLU implementation, there is no step to equilibrate the matrix. 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * For badly scaled matrices, this step can be useful to reduce the pivoting during factorization. 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If this is the case for your matrices, you can try the basic scaling method at 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * "unsupported/Eigen/src/IterativeSolvers/Scaling.h" 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<> 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS. Default is COLMAD 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa \ref TutorialSparseDirectSolvers 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa \ref OrderingMethods_Module 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename _MatrixType, typename _OrderingType> 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index> 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _MatrixType MatrixType; 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _OrderingType OrderingType; 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::RealScalar RealScalar; 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix; 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef internal::MappedSuperNodalMatrix<Scalar, Index> SCMatrix; 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,Dynamic,1> ScalarVector; 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Index,Dynamic,1> IndexVector; 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef internal::SparseLUImpl<Scalar, Index> Base; 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez initperfvalues(); 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLU(const MatrixType& matrix):m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez initperfvalues(); 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(matrix); 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ~SparseLU() 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Free all explicit dynamic pointers 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void analyzePattern (const MatrixType& matrix); 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void factorize (const MatrixType& matrix); 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void simplicialfactorize(const MatrixType& matrix); 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Compute the symbolic and numeric factorization of the input sparse matrix. 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The input matrix should be in column-major storage. 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void compute (const MatrixType& matrix) 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Analyze 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez analyzePattern(matrix); 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Factorize 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez factorize(matrix); 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index rows() const { return m_mat.rows(); } 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index cols() const { return m_mat.cols(); } 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Indicate that the pattern of the input matrix is symmetric */ 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void isSymmetric(bool sym) 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_symmetricmode = sym; 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns an expression of the matrix L, internally stored as supernodes 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The only operation available with this expression is the triangular solve 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \code 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * y = b; matrixL().solveInPlace(y); 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \endcode 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLUMatrixLReturnType<SCMatrix> matrixL() const 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore); 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns an expression of the matrix U, 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The only operation available with this expression is the triangular solve 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \code 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * y = b; matrixU().solveInPlace(y); 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \endcode 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,Index> > matrixU() const 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,Index> >(m_Lstore, m_Ustore); 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa colsPermutation() 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline const PermutationType& rowsPermutation() const 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_perm_r; 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa rowsPermutation() 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline const PermutationType& colsPermutation() const 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_perm_c; 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Set the threshold used for a diagonal entry to be an acceptable pivot. */ 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void setPivotThreshold(const RealScalar& thresh) 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_diagpivotthresh = thresh; 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \warning the destination matrix X in X = this->solve(B) must be colmun-major. 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute() 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Rhs> 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline const internal::solve_retval<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_factorizationIsOk && "SparseLU is not initialized."); 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(rows()==B.rows() 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && "SparseLU::solve(): invalid number of rows of the right hand side matrix B"); 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return internal::solve_retval<SparseLU, Rhs>(*this, B.derived()); 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute() 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Rhs> 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline const internal::sparse_solve_retval<SparseLU, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_factorizationIsOk && "SparseLU is not initialized."); 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(rows()==B.rows() 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && "SparseLU::solve(): invalid number of rows of the right hand side matrix B"); 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived()); 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Reports whether previous computation was successful. 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns \c Success if computation was succesful, 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \c InvalidInput if the input matrix is invalid 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa iparm() 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo info() const 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "Decomposition is not initialized."); 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_info; 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A string describing the type of error 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::string lastErrorMessage() const 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_lastError; 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Rhs, typename Dest> 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Dest& X(X_base.derived()); 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_factorizationIsOk && "The matrix should be factorized first"); 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Permute the right hand side to form X = Pr*B 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // on return, X is overwritten by the computed solution 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez X.resize(B.rows(),B.cols()); 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // this ugly const_cast_derived() helps to detect aliasing when applying the permutations 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index j = 0; j < B.cols(); ++j) 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez X.col(j) = rowsPermutation() * B.const_cast_derived().col(j); 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Forward substitution with L 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->matrixL().solveInPlace(X); 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->matrixU().solveInPlace(X); 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Permute back the solution 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j = 0; j < B.cols(); ++j) 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez X.col(j) = colsPermutation().inverse() * X.col(j); 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns the absolute value of the determinant of the matrix of which 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * *this is the QR decomposition. 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \warning a determinant can be very big or small, so for matrices 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * of large enough dimension, there is a risk of overflow/underflow. 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * One way to work around that is to use logAbsDeterminant() instead. 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa logAbsDeterminant(), signDeterminant() 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar absDeterminant() 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Initialize with the determinant of the row matrix 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar det = Scalar(1.); 263615d816d068b4d0f5e8df601930b5f160bf7eda1Tim Murray // Note that the diagonal blocks of U are stored in supernodes, 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // which are available in the L part :) 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j = 0; j < this->cols(); ++j) 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 269615d816d068b4d0f5e8df601930b5f160bf7eda1Tim Murray if(it.index() == j) 2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez det *= (std::abs)(it.value()); 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return det; 2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns the natural log of the absolute value of the determinant of the matrix 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * of which **this is the QR decomposition 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \note This method is useful to work around the risk of overflow/underflow that's 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * inherent to the determinant computation. 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa absDeterminant(), signDeterminant() 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar logAbsDeterminant() const 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar det = Scalar(0.); 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j = 0; j < this->cols(); ++j) 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(it.row() < j) continue; 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(it.row() == j) 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez det += (std::log)((std::abs)(it.value())); 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return det; 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns A number representing the sign of the determinant 3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa absDeterminant(), logAbsDeterminant() 3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar signDeterminant() 3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return Scalar(m_detPermR); 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez protected: 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Functions 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void initperfvalues() 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perfv.panel_size = 1; 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perfv.relax = 1; 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perfv.maxsuper = 128; 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perfv.rowblk = 16; 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perfv.colblk = 8; 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perfv.fillfactor = 20; 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Variables 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mutable ComputationInfo m_info; 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_isInitialized; 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_factorizationIsOk; 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_analysisIsOk; 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::string m_lastError; 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez NCMatrix m_mat; // The input (permuted ) matrix 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SCMatrix m_Lstore; // The lower triangular matrix (supernodal) 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore; // The upper triangular matrix 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez PermutationType m_perm_c; // Column permutation 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez PermutationType m_perm_r ; // Row permutation 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector m_etree; // Column elimination tree 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename Base::GlobalLU_t m_glu; 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // SparseLU options 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_symmetricmode; 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // values for performance 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::perfvalues<Index> m_perfv; 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m_nnzL, m_nnzU; // Nonzeros in L and U factors 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m_detPermR; // Determinant of the coefficient matrix 3507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez private: 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Disable copy constructor 3527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLU (const SparseLU& ); 3537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; // End class SparseLU 3557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Functions needed by the anaysis phase 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** 3607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Compute the column permutation to minimize the fill-in 3617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 3627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * - Apply this permutation to the input matrix - 3637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 3647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * - Compute the column elimination tree on the permuted matrix 3657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 3667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * - Postorder the elimination tree and the column permutation 3677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 3687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename MatrixType, typename OrderingType> 3707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) 3717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. 3747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OrderingType ord; 3767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ord(mat,m_perm_c); 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Apply the permutation to the column of the input matrix 3797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //First copy the whole input matrix. 3807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat = mat; 3817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_perm_c.size()) { 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. 3837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Then, permute only the column pointers 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index * outerIndexPtr; 3857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr(); 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index *outerIndexPtr_t = new Index[mat.cols()+1]; 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez outerIndexPtr = outerIndexPtr_t; 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i = 0; i < mat.cols(); i++) 3937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; 3967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(!mat.isCompressed()) delete[] outerIndexPtr; 3987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Compute the column elimination tree of the permuted matrix 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector firstRowElt; 4017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::coletree(m_mat, m_etree,firstRowElt); 4027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // In symmetric mode, do not do postorder here 4047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (!m_symmetricmode) { 4057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector post, iwork; 4067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Post order etree 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::treePostorder(m_mat.cols(), m_etree, post); 4087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Renumber etree in postorder 4117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m = m_mat.cols(); 4127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez iwork.resize(m+1); 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i)); 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_etree = iwork; 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez PermutationType post_perm(m); 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i = 0; i < m; i++) 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez post_perm.indices()(i) = post(i); 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Combine the two permutations : postorder the permutation for future use 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(m_perm_c.size()) { 4237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perm_c = post_perm * m_perm_c; 4247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end postordering 4277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_analysisIsOk = true; 4297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 4307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Functions needed by the numerical factorization phase 4327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * - Numerical factorization 4367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * - Interleaved with the symbolic factorization 4377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * On exit, info is 4387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * = 0: successful factorization 4407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * > 0: if info = i, and i is 4427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * <= A->ncol: U(i,i) is exactly zero. The factorization has 4447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * been completed, but the factor U is exactly singular, 4457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and division by zero will occur if it is used to solve a 4467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * system of equations. 4477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * > A->ncol: number of bytes allocated when memory allocation 4497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * failure occurred, plus A->ncol. If lwork = -1, it is 4507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * the estimated amount of space needed, plus A->ncol. 4517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 4527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename MatrixType, typename OrderingType> 4537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) 4547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 4557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using internal::emptyIdxLU; 4567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 4577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices"); 4587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename IndexVector::Scalar Index; 4607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Apply the column permutation computed in analyzepattern() 4637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // m_mat = matrix * m_perm_c.inverse(); 4647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat = matrix; 4657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_perm_c.size()) 4667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. 4687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Then, permute only the column pointers 4697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index * outerIndexPtr; 4707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr(); 4717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 4727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index* outerIndexPtr_t = new Index[matrix.cols()+1]; 4747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; 4757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez outerIndexPtr = outerIndexPtr_t; 4767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i = 0; i < matrix.cols(); i++) 4787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; 4807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; 4817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(!matrix.isCompressed()) delete[] outerIndexPtr; 4837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 4857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { //FIXME This should not be needed if the empty permutation is handled transparently 4867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perm_c.resize(matrix.cols()); 4877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i; 4887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m = m_mat.rows(); 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index n = m_mat.cols(); 4927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nnz = m_mat.nonZeros(); 4937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index maxpanel = m_perfv.panel_size * m; 4947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Allocate working storage common to the factor routines 4957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index lwork = 0; 4967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 4977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (info) 4987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ; 5007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_factorizationIsOk = false; 5017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return ; 5027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Set up pointers for integer working arrays 5057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector segrep(m); segrep.setZero(); 5067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector parent(m); parent.setZero(); 5077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector xplore(m); xplore.setZero(); 5087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector repfnz(maxpanel); 5097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector panel_lsub(maxpanel); 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector xprune(n); xprune.setZero(); 5117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector marker(m*internal::LUNoMarker); marker.setZero(); 5127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez repfnz.setConstant(-1); 5147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez panel_lsub.setConstant(-1); 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Set up pointers for scalar working arrays 5177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ScalarVector dense; 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dense.setZero(maxpanel); 5197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ScalarVector tempv; 5207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) ); 5217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Compute the inverse of perm_c 5237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez PermutationType iperm_c(m_perm_c.inverse()); 5247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Identify initial relaxed snodes 5267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector relax_end(n); 5277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if ( m_symmetricmode == true ) 5287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); 5297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 5307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); 5317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perm_r.resize(m); 5347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_perm_r.indices().setConstant(-1); 5357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez marker.setConstant(-1); 5367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_detPermR = 1; // Record the determinant of the row permutation 5377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0); 5397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0); 5407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Work on one 'panel' at a time. A panel is one of the following : 5427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // (a) a relaxed supernode at the bottom of the etree, or 5437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // (b) panel_size contiguous columns, <panel_size> defined by the user 5447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index jcol; 5457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez IndexVector panel_histo(n); 5467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index pivrow; // Pivotal row number in the original row matrix 5477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nseg1; // Number of segments in U-column above panel row jcol 5487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nseg; // Number of segments in each U-column 5497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index irep; 5507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index i, k, jj; 5517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (jcol = 0; jcol < n; ) 5527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Adjust panel size so that a panel won't overlap with the next relaxed snode. 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index panel_size = m_perfv.panel_size; // upper bound on panel width 5557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++) 5567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (relax_end(k) != emptyIdxLU) 5587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez panel_size = k - jcol; 5607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 5617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (k == n) 5647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez panel_size = n - jcol; 5657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Symbolic outer factorization on a panel of columns 5677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu); 5687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Numeric sup-panel updates in topological order 5707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 5717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Sparse LU within the panel, and below the panel diagonal 5737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for ( jj = jcol; jj< jcol + panel_size; jj++) 5747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez k = (jj - jcol) * m; // Column index for w-wide arrays 5767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez nseg = nseg1; // begin after all the panel segments 5787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Depth-first-search for the current column 5797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); 5807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorBlock<IndexVector> repfnz_k(repfnz, k, m); 5817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); 5827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if ( info ) 5837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() "; 5857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NumericalIssue; 5867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_factorizationIsOk = false; 5877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 5887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Numeric updates to this column 5907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorBlock<ScalarVector> dense_k(dense, k, m); 5917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1); 5927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 5937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if ( info ) 5947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() "; 5967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NumericalIssue; 5977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_factorizationIsOk = false; 5987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 5997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Copy the U-segments to ucol(*) 6027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 6037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if ( info ) 6047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() "; 6067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NumericalIssue; 6077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_factorizationIsOk = false; 6087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 6097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Form the L-segment 6127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); 6137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if ( info ) 6147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT "; 6167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::ostringstream returnInfo; 6177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez returnInfo << info; 6187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_lastError += returnInfo.str(); 6197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = NumericalIssue; 6207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_factorizationIsOk = false; 6217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 6227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Update the determinant of the row permutation matrix 6257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (pivrow != jj) m_detPermR *= -1; 6267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Prune columns (0:jj-1) using column jj 6287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 6297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Reset repfnz for this column 6317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (i = 0; i < nseg; i++) 6327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez irep = segrep(i); 6347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez repfnz_k(irep) = emptyIdxLU; 6357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end SparseLU within the panel 6377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez jcol += panel_size; // Move to the next panel 6387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end for -- end elimination 6397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Count the number of nonzeros in factors 6417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::countnz(n, m_nnzL, m_nnzU, m_glu); 6427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Apply permutation to the L subscripts 6437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Base::fixupL(n, m_perm_r.indices(), m_glu); 6447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Create supernode matrix L 6467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup); 6477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Create the column major upper sparse matrix U; 6487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez new (&m_Ustore) MappedSparseMatrix<Scalar, ColMajor, Index> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() ); 6497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = Success; 6517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_factorizationIsOk = true; 6527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 6537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MappedSupernodalType> 6557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SparseLUMatrixLReturnType : internal::no_assignment_operator 6567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MappedSupernodalType::Index Index; 6587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MappedSupernodalType::Scalar Scalar; 6597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL) 6607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { } 6617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows() { return m_mapL.rows(); } 6627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index cols() { return m_mapL.cols(); } 6637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> 6647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void solveInPlace( MatrixBase<Dest> &X) const 6657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_mapL.solveInPlace(X); 6677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MappedSupernodalType& m_mapL; 6697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 6707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixLType, typename MatrixUType> 6727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct SparseLUMatrixUReturnType : internal::no_assignment_operator 6737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixLType::Index Index; 6757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixLType::Scalar Scalar; 6767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU) 6777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_mapL(mapL),m_mapU(mapU) 6787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { } 6797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows() { return m_mapL.rows(); } 6807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index cols() { return m_mapL.cols(); } 6817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const 6837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nrhs = X.cols(); 6857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index n = X.rows(); 6867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Backward solve with U 6877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index k = m_mapL.nsuper(); k >= 0; k--) 6887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index fsupc = m_mapL.supToCol()[k]; 6907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension 6917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nsupc = m_mapL.supToCol()[k+1] - fsupc; 6927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index luptr = m_mapL.colIndexPtr()[fsupc]; 6937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (nsupc == 1) 6957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j = 0; j < nrhs; j++) 6977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez X(fsupc, j) /= m_mapL.valuePtr()[luptr]; 6997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 7027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 7047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 7057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez U = A.template triangularView<Upper>().solve(U); 7067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j = 0; j < nrhs; ++j) 7097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) 7117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename MatrixUType::InnerIterator it(m_mapU, jcol); 7137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for ( ; it; ++it) 7147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index irow = it.index(); 7167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez X(irow, j) -= X(jcol, j) * it.value(); 7177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // End For U-solve 7217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixLType& m_mapL; 7237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixUType& m_mapU; 7247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 7277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, typename Derived, typename Rhs> 7297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct solve_retval<SparseLU<_MatrixType,Derived>, Rhs> 7307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs> 7317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SparseLU<_MatrixType,Derived> Dec; 7337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 7347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void evalTo(Dest& dst) const 7367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dec()._solve(rhs(),dst); 7387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, typename Derived, typename Rhs> 7427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs> 7437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs> 7447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SparseLU<_MatrixType,Derived> Dec; 7467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 7477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void evalTo(Dest& dst) const 7497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->defaultEvalTo(dst); 7517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 7547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // End namespace Eigen 7567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif 758