SparseLU_SupernodalMatrix.h (12837B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 5 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@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() const { return m_row; } 79 80 /** 81 * Number of columns 82 */ 83 Index cols() const { 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 template<bool Conjugate, typename Dest> 160 void solveTransposedInPlace( MatrixBase<Dest>&X) const; 161 162 163 164 165 166 protected: 167 Index m_row; // Number of rows 168 Index m_col; // Number of columns 169 Index m_nsuper; // Number of supernodes 170 Scalar* m_nzval; //array of nonzero values packed by column 171 StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j 172 StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes 173 StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j 174 StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs 175 StorageIndex* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode 176 177 private : 178 }; 179 180 /** 181 * \brief InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L 182 * 183 */ 184 template<typename Scalar, typename StorageIndex> 185 class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator 186 { 187 public: 188 InnerIterator(const MappedSuperNodalMatrix& mat, Index outer) 189 : m_matrix(mat), 190 m_outer(outer), 191 m_supno(mat.colToSup()[outer]), 192 m_idval(mat.colIndexPtr()[outer]), 193 m_startidval(m_idval), 194 m_endidval(mat.colIndexPtr()[outer+1]), 195 m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]), 196 m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1]) 197 {} 198 inline InnerIterator& operator++() 199 { 200 m_idval++; 201 m_idrow++; 202 return *this; 203 } 204 inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; } 205 206 inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); } 207 208 inline Index index() const { return m_matrix.rowIndex()[m_idrow]; } 209 inline Index row() const { return index(); } 210 inline Index col() const { return m_outer; } 211 212 inline Index supIndex() const { return m_supno; } 213 214 inline operator bool() const 215 { 216 return ( (m_idval < m_endidval) && (m_idval >= m_startidval) 217 && (m_idrow < m_endidrow) ); 218 } 219 220 protected: 221 const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix 222 const Index m_outer; // Current column 223 const Index m_supno; // Current SuperNode number 224 Index m_idval; // Index to browse the values in the current column 225 const Index m_startidval; // Start of the column value 226 const Index m_endidval; // End of the column value 227 Index m_idrow; // Index to browse the row indices 228 Index m_endidrow; // End index of row indices of the current column 229 }; 230 231 /** 232 * \brief Solve with the supernode triangular matrix 233 * 234 */ 235 template<typename Scalar, typename Index_> 236 template<typename Dest> 237 void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) const 238 { 239 /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */ 240 // eigen_assert(X.rows() <= NumTraits<Index>::highest()); 241 // eigen_assert(X.cols() <= NumTraits<Index>::highest()); 242 Index n = int(X.rows()); 243 Index nrhs = Index(X.cols()); 244 const Scalar * Lval = valuePtr(); // Nonzero values 245 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector 246 work.setZero(); 247 for (Index k = 0; k <= nsuper(); k ++) 248 { 249 Index fsupc = supToCol()[k]; // First column of the current supernode 250 Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column 251 Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode 252 Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode 253 Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode 254 Index irow; //Current index row 255 256 if (nsupc == 1 ) 257 { 258 for (Index j = 0; j < nrhs; j++) 259 { 260 InnerIterator it(*this, fsupc); 261 ++it; // Skip the diagonal element 262 for (; it; ++it) 263 { 264 irow = it.row(); 265 X(irow, j) -= X(fsupc, j) * it.value(); 266 } 267 } 268 } 269 else 270 { 271 // The supernode has more than one column 272 Index luptr = colIndexPtr()[fsupc]; 273 Index lda = colIndexPtr()[fsupc+1] - luptr; 274 275 // Triangular solve 276 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 277 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 278 U = A.template triangularView<UnitLower>().solve(U); 279 280 // Matrix-vector product 281 new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); 282 work.topRows(nrow).noalias() = A * U; 283 284 //Begin Scatter 285 for (Index j = 0; j < nrhs; j++) 286 { 287 Index iptr = istart + nsupc; 288 for (Index i = 0; i < nrow; i++) 289 { 290 irow = rowIndex()[iptr]; 291 X(irow, j) -= work(i, j); // Scatter operation 292 work(i, j) = Scalar(0); 293 iptr++; 294 } 295 } 296 } 297 } 298 } 299 300 template<typename Scalar, typename Index_> 301 template<bool Conjugate, typename Dest> 302 void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<Dest>&X) const 303 { 304 using numext::conj; 305 Index n = int(X.rows()); 306 Index nrhs = Index(X.cols()); 307 const Scalar * Lval = valuePtr(); // Nonzero values 308 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector 309 work.setZero(); 310 for (Index k = nsuper(); k >= 0; k--) 311 { 312 Index fsupc = supToCol()[k]; // First column of the current supernode 313 Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column 314 Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode 315 Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode 316 Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode 317 Index irow; //Current index row 318 319 if (nsupc == 1 ) 320 { 321 for (Index j = 0; j < nrhs; j++) 322 { 323 InnerIterator it(*this, fsupc); 324 ++it; // Skip the diagonal element 325 for (; it; ++it) 326 { 327 irow = it.row(); 328 X(fsupc,j) -= X(irow, j) * (Conjugate?conj(it.value()):it.value()); 329 } 330 } 331 } 332 else 333 { 334 // The supernode has more than one column 335 Index luptr = colIndexPtr()[fsupc]; 336 Index lda = colIndexPtr()[fsupc+1] - luptr; 337 338 //Begin Gather 339 for (Index j = 0; j < nrhs; j++) 340 { 341 Index iptr = istart + nsupc; 342 for (Index i = 0; i < nrow; i++) 343 { 344 irow = rowIndex()[iptr]; 345 work.topRows(nrow)(i,j)= X(irow,j); // Gather operation 346 iptr++; 347 } 348 } 349 350 // Matrix-vector product with transposed submatrix 351 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); 352 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 353 if(Conjugate) 354 U = U - A.adjoint() * work.topRows(nrow); 355 else 356 U = U - A.transpose() * work.topRows(nrow); 357 358 // Triangular solve (of transposed diagonal block) 359 new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 360 if(Conjugate) 361 U = A.adjoint().template triangularView<UnitUpper>().solve(U); 362 else 363 U = A.transpose().template triangularView<UnitUpper>().solve(U); 364 365 } 366 367 } 368 } 369 370 371 } // end namespace internal 372 373 } // end namespace Eigen 374 375 #endif // EIGEN_SPARSELU_MATRIX_H