$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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()[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 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