![]() |
Eigen
3.3.3
|
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 _StorageIndex> 00033 class MappedSuperNodalMatrix 00034 { 00035 public: 00036 typedef _Scalar Scalar; 00037 typedef _StorageIndex StorageIndex; 00038 typedef Matrix<StorageIndex,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 StorageIndex* colIndexPtr() 00100 { 00101 return m_nzval_colptr; 00102 } 00103 00104 const StorageIndex* colIndexPtr() const 00105 { 00106 return m_nzval_colptr; 00107 } 00108 00112 StorageIndex* rowIndex() { return m_rowind; } 00113 00114 const StorageIndex* rowIndex() const 00115 { 00116 return m_rowind; 00117 } 00118 00122 StorageIndex* rowIndexPtr() { return m_rowind_colptr; } 00123 00124 const StorageIndex* rowIndexPtr() const 00125 { 00126 return m_rowind_colptr; 00127 } 00128 00132 StorageIndex* colToSup() { return m_col_to_sup; } 00133 00134 const StorageIndex* colToSup() const 00135 { 00136 return m_col_to_sup; 00137 } 00141 StorageIndex* supToCol() { return m_sup_to_col; } 00142 00143 const StorageIndex* 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 StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j 00169 StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes 00170 StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j 00171 StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs 00172 StorageIndex* 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 StorageIndex> 00182 class MappedSuperNodalMatrix<Scalar,StorageIndex>::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()[mat.supToCol()[mat.colToSup()[outer]]]), 00193 m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[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 /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */ 00237 // eigen_assert(X.rows() <= NumTraits<Index>::highest()); 00238 // eigen_assert(X.cols() <= NumTraits<Index>::highest()); 00239 Index n = int(X.rows()); 00240 Index nrhs = Index(X.cols()); 00241 const Scalar * Lval = valuePtr(); // Nonzero values 00242 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector 00243 work.setZero(); 00244 for (Index k = 0; k <= nsuper(); k ++) 00245 { 00246 Index fsupc = supToCol()[k]; // First column of the current supernode 00247 Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column 00248 Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode 00249 Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode 00250 Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode 00251 Index irow; //Current index row 00252 00253 if (nsupc == 1 ) 00254 { 00255 for (Index j = 0; j < nrhs; j++) 00256 { 00257 InnerIterator it(*this, fsupc); 00258 ++it; // Skip the diagonal element 00259 for (; it; ++it) 00260 { 00261 irow = it.row(); 00262 X(irow, j) -= X(fsupc, j) * it.value(); 00263 } 00264 } 00265 } 00266 else 00267 { 00268 // The supernode has more than one column 00269 Index luptr = colIndexPtr()[fsupc]; 00270 Index lda = colIndexPtr()[fsupc+1] - luptr; 00271 00272 // Triangular solve 00273 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 00274 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 00275 U = A.template triangularView<UnitLower>().solve(U); 00276 00277 // Matrix-vector product 00278 new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); 00279 work.topRows(nrow).noalias() = A * U; 00280 00281 //Begin Scatter 00282 for (Index j = 0; j < nrhs; j++) 00283 { 00284 Index iptr = istart + nsupc; 00285 for (Index i = 0; i < nrow; i++) 00286 { 00287 irow = rowIndex()[iptr]; 00288 X(irow, j) -= work(i, j); // Scatter operation 00289 work(i, j) = Scalar(0); 00290 iptr++; 00291 } 00292 } 00293 } 00294 } 00295 } 00296 00297 } // end namespace internal 00298 00299 } // end namespace Eigen 00300 00301 #endif // EIGEN_SPARSELU_MATRIX_H