1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 5 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H 12 #define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H 13 14 namespace Eigen { 15 namespace internal { 16 17 /** \ingroup SparseLU_Module 18 * \brief a class to manipulate the L supernodal factor from the SparseLU factorization 19 * 20 * This class contain the data to easily store 21 * and manipulate the supernodes during the factorization and solution phase of Sparse LU. 22 * Only the lower triangular matrix has supernodes. 23 * 24 * NOTE : This class corresponds to the SCformat structure in SuperLU 25 * 26 */ 27 /* TODO 28 * InnerIterator as for sparsematrix 29 * SuperInnerIterator to iterate through all supernodes 30 * Function for triangular solve 31 */ 32 template <typename _Scalar, typename _StorageIndex> 33 class MappedSuperNodalMatrix 34 { 35 public: 36 typedef _Scalar Scalar; 37 typedef _StorageIndex StorageIndex; 38 typedef Matrix<StorageIndex,Dynamic,1> IndexVector; 39 typedef Matrix<Scalar,Dynamic,1> ScalarVector; 40 public: 41 MappedSuperNodalMatrix() 42 { 43 44 } 45 MappedSuperNodalMatrix(Index m, Index n, ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind, 46 IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col ) 47 { 48 setInfos(m, n, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col); 49 } 50 51 ~MappedSuperNodalMatrix() 52 { 53 54 } 55 /** 56 * Set appropriate pointers for the lower triangular supernodal matrix 57 * These infos are available at the end of the numerical factorization 58 * FIXME This class will be modified such that it can be use in the course 59 * of the factorization. 60 */ 61 void setInfos(Index m, Index n, ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind, 62 IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col ) 63 { 64 m_row = m; 65 m_col = n; 66 m_nzval = nzval.data(); 67 m_nzval_colptr = nzval_colptr.data(); 68 m_rowind = rowind.data(); 69 m_rowind_colptr = rowind_colptr.data(); 70 m_nsuper = col_to_sup(n); 71 m_col_to_sup = col_to_sup.data(); 72 m_sup_to_col = sup_to_col.data(); 73 } 74 75 /** 76 * Number of rows 77 */ 78 Index rows() { return m_row; } 79 80 /** 81 * Number of columns 82 */ 83 Index cols() { return m_col; } 84 85 /** 86 * Return the array of nonzero values packed by column 87 * 88 * The size is nnz 89 */ 90 Scalar* valuePtr() { return m_nzval; } 91 92 const Scalar* valuePtr() const 93 { 94 return m_nzval; 95 } 96 /** 97 * Return the pointers to the beginning of each column in \ref valuePtr() 98 */ 99 StorageIndex* colIndexPtr() 100 { 101 return m_nzval_colptr; 102 } 103 104 const StorageIndex* colIndexPtr() const 105 { 106 return m_nzval_colptr; 107 } 108 109 /** 110 * Return the array of compressed row indices of all supernodes 111 */ 112 StorageIndex* rowIndex() { return m_rowind; } 113 114 const StorageIndex* rowIndex() const 115 { 116 return m_rowind; 117 } 118 119 /** 120 * Return the location in \em rowvaluePtr() which starts each column 121 */ 122 StorageIndex* rowIndexPtr() { return m_rowind_colptr; } 123 124 const StorageIndex* rowIndexPtr() const 125 { 126 return m_rowind_colptr; 127 } 128 129 /** 130 * Return the array of column-to-supernode mapping 131 */ 132 StorageIndex* colToSup() { return m_col_to_sup; } 133 134 const StorageIndex* colToSup() const 135 { 136 return m_col_to_sup; 137 } 138 /** 139 * Return the array of supernode-to-column mapping 140 */ 141 StorageIndex* supToCol() { return m_sup_to_col; } 142 143 const StorageIndex* supToCol() const 144 { 145 return m_sup_to_col; 146 } 147 148 /** 149 * Return the number of supernodes 150 */ 151 Index nsuper() const 152 { 153 return m_nsuper; 154 } 155 156 class InnerIterator; 157 template<typename Dest> 158 void solveInPlace( MatrixBase<Dest>&X) const; 159 160 161 162 163 protected: 164 Index m_row; // Number of rows 165 Index m_col; // Number of columns 166 Index m_nsuper; // Number of supernodes 167 Scalar* m_nzval; //array of nonzero values packed by column 168 StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j 169 StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes 170 StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j 171 StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs 172 StorageIndex* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode 173 174 private : 175 }; 176 177 /** 178 * \brief InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L 179 * 180 */ 181 template<typename Scalar, typename StorageIndex> 182 class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator 183 { 184 public: 185 InnerIterator(const MappedSuperNodalMatrix& mat, Index outer) 186 : m_matrix(mat), 187 m_outer(outer), 188 m_supno(mat.colToSup()[outer]), 189 m_idval(mat.colIndexPtr()[outer]), 190 m_startidval(m_idval), 191 m_endidval(mat.colIndexPtr()[outer+1]), 192 m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]), 193 m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1]) 194 {} 195 inline InnerIterator& operator++() 196 { 197 m_idval++; 198 m_idrow++; 199 return *this; 200 } 201 inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; } 202 203 inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); } 204 205 inline Index index() const { return m_matrix.rowIndex()[m_idrow]; } 206 inline Index row() const { return index(); } 207 inline Index col() const { return m_outer; } 208 209 inline Index supIndex() const { return m_supno; } 210 211 inline operator bool() const 212 { 213 return ( (m_idval < m_endidval) && (m_idval >= m_startidval) 214 && (m_idrow < m_endidrow) ); 215 } 216 217 protected: 218 const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix 219 const Index m_outer; // Current column 220 const Index m_supno; // Current SuperNode number 221 Index m_idval; // Index to browse the values in the current column 222 const Index m_startidval; // Start of the column value 223 const Index m_endidval; // End of the column value 224 Index m_idrow; // Index to browse the row indices 225 Index m_endidrow; // End index of row indices of the current column 226 }; 227 228 /** 229 * \brief Solve with the supernode triangular matrix 230 * 231 */ 232 template<typename Scalar, typename Index_> 233 template<typename Dest> 234 void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) const 235 { 236 /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */ 237 // eigen_assert(X.rows() <= NumTraits<Index>::highest()); 238 // eigen_assert(X.cols() <= NumTraits<Index>::highest()); 239 Index n = int(X.rows()); 240 Index nrhs = Index(X.cols()); 241 const Scalar * Lval = valuePtr(); // Nonzero values 242 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector 243 work.setZero(); 244 for (Index k = 0; k <= nsuper(); k ++) 245 { 246 Index fsupc = supToCol()[k]; // First column of the current supernode 247 Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column 248 Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode 249 Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode 250 Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode 251 Index irow; //Current index row 252 253 if (nsupc == 1 ) 254 { 255 for (Index j = 0; j < nrhs; j++) 256 { 257 InnerIterator it(*this, fsupc); 258 ++it; // Skip the diagonal element 259 for (; it; ++it) 260 { 261 irow = it.row(); 262 X(irow, j) -= X(fsupc, j) * it.value(); 263 } 264 } 265 } 266 else 267 { 268 // The supernode has more than one column 269 Index luptr = colIndexPtr()[fsupc]; 270 Index lda = colIndexPtr()[fsupc+1] - luptr; 271 272 // Triangular solve 273 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 274 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 275 U = A.template triangularView<UnitLower>().solve(U); 276 277 // Matrix-vector product 278 new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); 279 work.topRows(nrow).noalias() = A * U; 280 281 //Begin Scatter 282 for (Index j = 0; j < nrhs; j++) 283 { 284 Index iptr = istart + nsupc; 285 for (Index i = 0; i < nrow; i++) 286 { 287 irow = rowIndex()[iptr]; 288 X(irow, j) -= work(i, j); // Scatter operation 289 work(i, j) = Scalar(0); 290 iptr++; 291 } 292 } 293 } 294 } 295 } 296 297 } // end namespace internal 298 299 } // end namespace Eigen 300 301 #endif // EIGEN_SPARSELU_MATRIX_H 302
