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#ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen {
157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal {
167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \ingroup SparseLU_Module
187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief a class to manipulate the L supernodal factor from the SparseLU factorization
197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This class  contain the data to easily store
217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and manipulate the supernodes during the factorization and solution phase of Sparse LU.
227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Only the lower triangular matrix has supernodes.
237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE : This class corresponds to the SCformat structure in SuperLU
257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */
277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/* TODO
287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * InnerIterator as for sparsematrix
297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * SuperInnerIterator to iterate through all supernodes
307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Function for triangular solve
317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */
327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate <typename _Scalar, typename _Index>
337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass MappedSuperNodalMatrix
347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  public:
367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef _Scalar Scalar;
377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef _Index Index;
387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef Matrix<Index,Dynamic,1> IndexVector;
397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    typedef Matrix<Scalar,Dynamic,1> ScalarVector;
407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  public:
417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    MappedSuperNodalMatrix()
427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    MappedSuperNodalMatrix(Index m, Index n,  ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind,
467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez             IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      setInfos(m, n, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col);
497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    ~MappedSuperNodalMatrix()
527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Set appropriate pointers for the lower triangular supernodal matrix
577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * These infos are available at the end of the numerical factorization
587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * FIXME This class will be modified such that it can be use in the course
597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * of the factorization.
607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void setInfos(Index m, Index n, ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind,
627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez             IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_row = m;
657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_col = n;
667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_nzval = nzval.data();
677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_nzval_colptr = nzval_colptr.data();
687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_rowind = rowind.data();
697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_rowind_colptr = rowind_colptr.data();
707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_nsuper = col_to_sup(n);
717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_col_to_sup = col_to_sup.data();
727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_sup_to_col = sup_to_col.data();
737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Number of rows
777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index rows() { return m_row; }
797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Number of columns
827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index cols() { return m_col; }
847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the array of nonzero values packed by column
877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     *
887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * The size is nnz
897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Scalar* valuePtr() {  return m_nzval; }
917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Scalar* valuePtr() const
937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_nzval;
957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the pointers to the beginning of each column in \ref valuePtr()
987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* colIndexPtr()
1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_nzval_colptr;
1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index* colIndexPtr() const
1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_nzval_colptr;
1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the array of compressed row indices of all supernodes
1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* rowIndex()  { return m_rowind; }
1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index* rowIndex() const
1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_rowind;
1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the location in \em rowvaluePtr() which starts each column
1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* rowIndexPtr() { return m_rowind_colptr; }
1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index* rowIndexPtr() const
1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_rowind_colptr;
1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the array of column-to-supernode mapping
1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* colToSup()  { return m_col_to_sup; }
1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index* colToSup() const
1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_col_to_sup;
1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the array of supernode-to-column mapping
1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* supToCol() { return m_sup_to_col; }
1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index* supToCol() const
1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_sup_to_col;
1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    /**
1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     * Return the number of supernodes
1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     */
1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index nsuper() const
1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return m_nsuper;
1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    class InnerIterator;
1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    template<typename Dest>
1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    void solveInPlace( MatrixBase<Dest>&X) const;
1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  protected:
1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index m_row; // Number of rows
1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index m_col; // Number of columns
1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index m_nsuper; // Number of supernodes
1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Scalar* m_nzval; //array of nonzero values packed by column
1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* m_rowind; // Array of compressed row indices of rectangular supernodes
1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  private :
1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/**
1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  * \brief InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L
1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  *
1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  */
1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Scalar, typename Index>
1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  public:
1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez     InnerIterator(const MappedSuperNodalMatrix& mat, Index outer)
1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      : m_matrix(mat),
1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_outer(outer),
1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_supno(mat.colToSup()[outer]),
1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_idval(mat.colIndexPtr()[outer]),
1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_startidval(m_idval),
1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        m_endidval(mat.colIndexPtr()[outer+1]),
192615d816d068b4d0f5e8df601930b5f160bf7eda1Tim Murray        m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
193615d816d068b4d0f5e8df601930b5f160bf7eda1Tim Murray        m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {}
1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline InnerIterator& operator++()
1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_idval++;
1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      m_idrow++;
1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return *this;
2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; }
2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); }
2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Index index() const { return m_matrix.rowIndex()[m_idrow]; }
2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Index row() const { return index(); }
2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Index col() const { return m_outer; }
2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline Index supIndex() const { return m_supno; }
2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    inline operator bool() const
2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez                && (m_idrow < m_endidrow) );
2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  protected:
2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix
2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index m_outer;                    // Current column
2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index m_supno;                    // Current SuperNode number
2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index m_idval;                          // Index to browse the values in the current column
2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index m_startidval;               // Start of the column value
2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Index m_endidval;                 // End of the column value
2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index m_idrow;                          // Index to browse the row indices
2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index m_endidrow;                       // End index of row indices of the current column
2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez};
2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/**
2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Solve with the supernode triangular matrix
2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *
2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */
2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Scalar, typename Index>
2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Dest>
2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid MappedSuperNodalMatrix<Scalar,Index>::solveInPlace( MatrixBase<Dest>&X) const
2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{
2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index n = X.rows();
2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Index nrhs = X.cols();
2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    const Scalar * Lval = valuePtr();                 // Nonzero values
2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    Matrix<Scalar,Dynamic,Dynamic> work(n, nrhs);     // working vector
2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    work.setZero();
2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    for (Index k = 0; k <= nsuper(); k ++)
2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    {
2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index fsupc = supToCol()[k];                    // First column of the current supernode
2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index istart = rowIndexPtr()[fsupc];            // Pointer index to the subscript of the current column
2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index nsupr = rowIndexPtr()[fsupc+1] - istart;  // Number of rows in the current supernode
2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index nsupc = supToCol()[k+1] - fsupc;          // Number of columns in the current supernode
2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index nrow = nsupr - nsupc;                     // Number of rows in the non-diagonal part of the supernode
2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      Index irow;                                     //Current index row
2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      if (nsupc == 1 )
2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      {
2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        for (Index j = 0; j < nrhs; j++)
2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        {
2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          InnerIterator it(*this, fsupc);
2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          ++it; // Skip the diagonal element
2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          for (; it; ++it)
2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          {
2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez            irow = it.row();
2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez            X(irow, j) -= X(fsupc, j) * it.value();
2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          }
2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        }
2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      }
2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      else
2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      {
2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // The supernode has more than one column
2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        Index luptr = colIndexPtr()[fsupc];
2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        Index lda = colIndexPtr()[fsupc+1] - luptr;
2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // Triangular solve
2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        U = A.template triangularView<UnitLower>().solve(U);
2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        // Matrix-vector product
2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        work.block(0, 0, nrow, nrhs) = A * U;
2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        //Begin Scatter
2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        for (Index j = 0; j < nrhs; j++)
2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        {
2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          Index iptr = istart + nsupc;
2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          for (Index i = 0; i < nrow; i++)
2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          {
2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez            irow = rowIndex()[iptr];
2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez            X(irow, j) -= work(i, j); // Scatter operation
2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez            work(i, j) = Scalar(0);
2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez            iptr++;
2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez          }
2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez        }
2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez      }
2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez    }
2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}
2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal
2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen
2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez
2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // EIGEN_SPARSELU_MATRIX_H
299