SparseLU_SupernodalMatrix.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00005 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
00012 #define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
00013 
00014 namespace Eigen {
00015 namespace internal {
00016 
00027 /* TODO
00028  * InnerIterator as for sparsematrix 
00029  * SuperInnerIterator to iterate through all supernodes 
00030  * Function for triangular solve
00031  */
00032 template <typename _Scalar, typename _Index>
00033 class MappedSuperNodalMatrix
00034 {
00035   public:
00036     typedef _Scalar Scalar; 
00037     typedef _Index Index;
00038     typedef Matrix<Index,Dynamic,1> IndexVector; 
00039     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
00040   public:
00041     MappedSuperNodalMatrix()
00042     {
00043       
00044     }
00045     MappedSuperNodalMatrix(Index m, Index n,  ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind, 
00046              IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
00047     {
00048       setInfos(m, n, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col);
00049     }
00050     
00051     ~MappedSuperNodalMatrix()
00052     {
00053       
00054     }
00061     void setInfos(Index m, Index n, ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind, 
00062              IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
00063     {
00064       m_row = m;
00065       m_col = n; 
00066       m_nzval = nzval.data(); 
00067       m_nzval_colptr = nzval_colptr.data(); 
00068       m_rowind = rowind.data(); 
00069       m_rowind_colptr = rowind_colptr.data(); 
00070       m_nsuper = col_to_sup(n); 
00071       m_col_to_sup = col_to_sup.data(); 
00072       m_sup_to_col = sup_to_col.data(); 
00073     }
00074     
00078     Index rows() { return m_row; }
00079     
00083     Index cols() { return m_col; }
00084     
00090     Scalar* valuePtr() {  return m_nzval; }
00091     
00092     const Scalar* valuePtr() const 
00093     {
00094       return m_nzval; 
00095     }
00099     Index* colIndexPtr()
00100     {
00101       return m_nzval_colptr; 
00102     }
00103     
00104     const Index* colIndexPtr() const
00105     {
00106       return m_nzval_colptr; 
00107     }
00108     
00112     Index* rowIndex()  { return m_rowind; }
00113     
00114     const Index* rowIndex() const
00115     {
00116       return m_rowind; 
00117     }
00118     
00122     Index* rowIndexPtr() { return m_rowind_colptr; }
00123     
00124     const Index* rowIndexPtr() const 
00125     {
00126       return m_rowind_colptr; 
00127     }
00128     
00132     Index* colToSup()  { return m_col_to_sup; }
00133     
00134     const Index* colToSup() const
00135     {
00136       return m_col_to_sup;       
00137     }
00141     Index* supToCol() { return m_sup_to_col; }
00142     
00143     const Index* supToCol() const 
00144     {
00145       return m_sup_to_col;
00146     }
00147     
00151     Index nsuper() const 
00152     {
00153       return m_nsuper; 
00154     }
00155     
00156     class InnerIterator; 
00157     template<typename Dest>
00158     void solveInPlace( MatrixBase<Dest>&X) const;
00159     
00160       
00161       
00162     
00163   protected:
00164     Index m_row; // Number of rows
00165     Index m_col; // Number of columns 
00166     Index m_nsuper; // Number of supernodes 
00167     Scalar* m_nzval; //array of nonzero values packed by column
00168     Index* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j 
00169     Index* m_rowind; // Array of compressed row indices of rectangular supernodes
00170     Index* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
00171     Index* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
00172     Index* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
00173     
00174   private :
00175 };
00176 
00181 template<typename Scalar, typename Index>
00182 class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
00183 {
00184   public:
00185      InnerIterator(const MappedSuperNodalMatrix& mat, Index outer)
00186       : m_matrix(mat),
00187         m_outer(outer), 
00188         m_supno(mat.colToSup()[outer]),
00189         m_idval(mat.colIndexPtr()[outer]),
00190         m_startidval(m_idval),
00191         m_endidval(mat.colIndexPtr()[outer+1]),
00192         m_idrow(mat.rowIndexPtr()[outer]),
00193         m_endidrow(mat.rowIndexPtr()[outer+1])
00194     {}
00195     inline InnerIterator& operator++()
00196     { 
00197       m_idval++; 
00198       m_idrow++;
00199       return *this;
00200     }
00201     inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; }
00202     
00203     inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); }
00204     
00205     inline Index index() const { return m_matrix.rowIndex()[m_idrow]; }
00206     inline Index row() const { return index(); }
00207     inline Index col() const { return m_outer; }
00208     
00209     inline Index supIndex() const { return m_supno; }
00210     
00211     inline operator bool() const 
00212     { 
00213       return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
00214                 && (m_idrow < m_endidrow) );
00215     }
00216     
00217   protected:
00218     const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix 
00219     const Index m_outer;                    // Current column 
00220     const Index m_supno;                    // Current SuperNode number
00221     Index m_idval;                          // Index to browse the values in the current column
00222     const Index m_startidval;               // Start of the column value
00223     const Index m_endidval;                 // End of the column value
00224     Index m_idrow;                          // Index to browse the row indices 
00225     Index m_endidrow;                       // End index of row indices of the current column
00226 };
00227 
00232 template<typename Scalar, typename Index>
00233 template<typename Dest>
00234 void MappedSuperNodalMatrix<Scalar,Index>::solveInPlace( MatrixBase<Dest>&X) const
00235 {
00236     Index n = X.rows(); 
00237     Index nrhs = X.cols(); 
00238     const Scalar * Lval = valuePtr();                 // Nonzero values 
00239     Matrix<Scalar,Dynamic,Dynamic> work(n, nrhs);     // working vector
00240     work.setZero();
00241     for (Index k = 0; k <= nsuper(); k ++)
00242     {
00243       Index fsupc = supToCol()[k];                    // First column of the current supernode 
00244       Index istart = rowIndexPtr()[fsupc];            // Pointer index to the subscript of the current column
00245       Index nsupr = rowIndexPtr()[fsupc+1] - istart;  // Number of rows in the current supernode
00246       Index nsupc = supToCol()[k+1] - fsupc;          // Number of columns in the current supernode
00247       Index nrow = nsupr - nsupc;                     // Number of rows in the non-diagonal part of the supernode
00248       Index irow;                                     //Current index row
00249       
00250       if (nsupc == 1 )
00251       {
00252         for (Index j = 0; j < nrhs; j++)
00253         {
00254           InnerIterator it(*this, fsupc);
00255           ++it; // Skip the diagonal element
00256           for (; it; ++it)
00257           {
00258             irow = it.row();
00259             X(irow, j) -= X(fsupc, j) * it.value();
00260           }
00261         }
00262       }
00263       else
00264       {
00265         // The supernode has more than one column 
00266         Index luptr = colIndexPtr()[fsupc]; 
00267         Index lda = colIndexPtr()[fsupc+1] - luptr;
00268         
00269         // Triangular solve 
00270         Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 
00271         Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 
00272         U = A.template triangularView<UnitLower>().solve(U); 
00273         
00274         // Matrix-vector product 
00275         new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); 
00276         work.block(0, 0, nrow, nrhs) = A * U; 
00277         
00278         //Begin Scatter 
00279         for (Index j = 0; j < nrhs; j++)
00280         {
00281           Index iptr = istart + nsupc; 
00282           for (Index i = 0; i < nrow; i++)
00283           {
00284             irow = rowIndex()[iptr]; 
00285             X(irow, j) -= work(i, j); // Scatter operation
00286             work(i, j) = Scalar(0); 
00287             iptr++;
00288           }
00289         }
00290       }
00291     } 
00292 }
00293 
00294 } // end namespace internal
00295 
00296 } // end namespace Eigen
00297 
00298 #endif // EIGEN_SPARSELU_MATRIX_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:39:23