18 #ifndef la_matrix_coo_h
19 #define la_matrix_coo_h
21 #define ALLOWED_MEM_OVERHEAD 2.e+6
33 template<
typename INDT,
typename VALT >
56 template<
typename INDT,
typename VALT >
117 typename MatrixCooVec_::const_iterator itGlobale =
matrixCooVec_.end();
119 typename MatrixCooVec_::const_iterator itSelectb =
matrixCooVec_.begin();
120 typename MatrixCooVec_::const_iterator itSelecte = itSelectb;
124 typename MatrixCooVec_::iterator itStore =
matrixCooVec_.begin();
127 while (itSelectb != itGlobale) {
130 itSelecte = std::upper_bound( itSelectb, itGlobale, *itSelectb,
TripleLessThan_( ) );
136 *itStore = tmpTriple;
140 itSelectb = itSelecte;
167 typename MatrixCooVec_::const_iterator iter =
matrixCooVec_.begin();
168 typename MatrixCooVec_::const_iterator itere =
matrixCooVec_.end();
170 for ( ; iter != itere; ++iter ) {
176 if ( row >= low and row < high ) {
181 if ( col >= low and col < high ) {
187 nnzRowsOut[rowLoc]++;
211 typename MatrixCooVec_::iterator tripEraseIter;
212 typename MatrixCooVec_::iterator tripEraseIterE;
219 tripEraseIterE = std::upper_bound( tripEraseIter,
matrixCooVec_.end(),
223 typename MatrixCooVec_::iterator tripEraseIterDiag;
224 tripEraseIterDiag = std::find_if( tripEraseIter, tripEraseIterE,
TripleDiagEntry_( ) );
227 for ( ; tripEraseIter != tripEraseIterE; ++tripEraseIter ) {
228 (*tripEraseIter).second =
static_cast< ValueType >( 0.0 );
232 Triple_ tmpTrip = std::make_pair( std::make_pair ( row, row ), factor );
234 if ( tripEraseIterDiag != tripEraseIterE ) {
236 *tripEraseIterDiag = tmpTrip;
259 typename MatrixCooVec_::const_iterator tvIter =
matrixCooVec_.begin();
260 typename MatrixCooVec_::const_iterator tvEnd =
matrixCooVec_.end();
267 rowPtrs.reserve( nr + 1 );
269 int previousRow = -1;
272 for ( ; tvIter != tvEnd; ++tvIter ) {
273 const IndexType row = (*tvIter).first.first;
274 const IndexType col = (*tvIter).first.second;
276 values.push_back( tvIter -> second );
277 colIndices.push_back( col );
278 if ( row == previousRow ) index++;
279 else rowPtrs.push_back( index++ );
284 rowPtrs.push_back( values.size() );
293 bool onlyUpperTriangle =
false )
303 typename MatrixCooVec_::const_iterator tvIter =
matrixCooVec_.begin();
304 typename MatrixCooVec_::const_iterator tvEnd =
matrixCooVec_.end();
311 for ( ; tvIter != tvEnd; ++tvIter ) {
312 const IndexType row = (*tvIter).first.first ;
313 const IndexType col = (*tvIter).first.second ;
315 if ( not onlyUpperTriangle or col >= row ) {
316 values.push_back( (*tvIter).second );
317 rowIndices.push_back( row );
318 colIndices.push_back( col );
337 typename MatrixCooVec_::iterator tripEraseIter =
matrixCooVec_.begin( );
338 typename MatrixCooVec_::iterator tripEraseIterE =
matrixCooVec_.end( );
341 for ( ; tripEraseIter != tripEraseIterE; ++tripEraseIter ) {
342 (*tripEraseIter).second =
static_cast< ValueType >(0.0);
349 std::ostream &
write( std::ostream & out )
355 typename MatrixCooVec_::const_iterator tvIter =
matrixCooVec_.begin();
356 typename MatrixCooVec_::const_iterator tvEnd =
matrixCooVec_.end();
357 for ( ; tvIter != tvEnd; ++tvIter ) {
358 const IndexType row = (*tvIter).first.first;
359 const IndexType col = (*tvIter).first.second;
361 out << row <<
" " << col <<
" "
362 << std::setprecision(12)
363 << (*tvIter).second << std::endl;
379 typename MatrixCooVec_::const_iterator iter =
matrixCooVec_.begin();
380 typename MatrixCooVec_::const_iterator last =
matrixCooVec_.end();
381 for ( ; iter != last; ++iter ) {
382 const IndexType iRow = (*iter).first.first;
383 rowSums[ iRow ] += std::fabs( (*iter).second );
386 return *std::max_element( rowSums.begin(), rowSums.end() );
400 typename MatrixCooVec_::const_iterator iter =
matrixCooVec_.begin();
401 typename MatrixCooVec_::const_iterator last =
matrixCooVec_.end();
402 for ( ; iter != last; ++iter ) {
405 sumOfSquares += val * val;
408 return std::sqrt( sumOfSquares );
439 if ( val > maxVal ) {
444 if ( val < minVal ) {
451 ValueType scalar = ( minVal + maxVal ) / 2.0;
459 typedef std::pair<IndexPair_,ValueType>
Triple_;
497 IndexPair_ firstTripleIndxPair = firstTriple.first;
498 IndexPair_ secondTripleIndxPair = secondTriple.first;
500 return (firstTripleIndxPair < secondTripleIndxPair);
509 IndexType firstRowIndex = firstTriple.first.first;
510 IndexType secondRowIndex = secondTriple.first.first;
512 return (firstRowIndex < secondRowIndex);
521 IndexType tripleRowIndex = triple.first.first;
523 return (tripleRowIndex < row);
532 IndexType tripleRowIndex = triple.first.first;
533 IndexType tripleColIndex = triple.first.second;
535 return (tripleRowIndex == tripleColIndex);
545 for ( ; itb != ite; ++itb) {
546 val += (*itb).second;
549 Triple_ combinedTriple = std::make_pair(tmpPair, val);
551 return combinedTriple;
Temporary storage for a sparse matrix.
unsigned sortedSize_
number of sorted entries
IndexType deduceNumRows_()
void clear()
destroy all contents
void getIntervalCounts(IndexType low, IndexType high, std::vector< IndexType > &nnzRowsIn, std::vector< IndexType > &nnzRowsOut)
Determine number of nonzeros within a block [low,high) - used for PETSc preallocation.
bool sorted_
state variable - is matrix sorted?
std::vector< Triple_ > MatrixCooVec_
vector of triples
void prepareAssembly(const IndexType length)
void zeroRow(const IndexType row, const ValueType factor)
Zero a row and set diagonal to scalar.
IndexType nRows_
number of rows of matrix
void insert(const IndexType i, const IndexType j, const ValueType value)
Insert a value to unassembled vector.
std::ostream & write(std::ostream &out)
debug-output
void erase()
erase all contents
IndexType nnz() const
Return number of non-zeros.
IndexType nCols_
number of columns of matrix
ValueType frobeniusNorm()
Return Frobenius norm of matrix .
double infNorm()
Compute maximal row sum (infinity-norm) of a matrix in coordinate format.
MatrixCooVec_ matrixCooVec_
matrix in coordinate format <i,j, values>
Triple_ combineTriples_(typename MatrixCooVec_::const_iterator itb, typename MatrixCooVec_::const_iterator ite)
merge range of triples by addition
std::pair< IndexType, IndexType > IndexPair_
pair of indices
std::pair< IndexPair_, ValueType > Triple_
triple
void extractArrays(std::vector< IndexType > &rowIndices, std::vector< IndexType > &colIndices, std::vector< ValueType > &values, bool onlyUpperTriangle=false)
Fill independent arrays with column indices, row pointers and values.
MatrixCoo(unsigned nRows=0, unsigned nCols=0)
void getEntry(IndexType index, IndexType &row, IndexType &col, ValueType &val)
Get specific entry.
void setCompressedRowArrays(std::vector< IndexType > &colIndices, std::vector< IndexType > &rowPtrs, std::vector< ValueType > &values)
void finishAssembly()
Finalize matrix assembly - build the assembled input and truncate matrix.
static constexpr bool value
#define ALLOWED_MEM_OVERHEAD
returns true if the entry is on diagonal, i.e. row == column
bool operator()(Triple_ triple)
returns true if the index pair of the first entry is smaller than of the second entry
bool operator()(Triple_ firstTriple, Triple_ secondTriple)
returns true if the row of an entry is smaller than a prescribed index
bool operator()(Triple_ triple, IndexType row)
returns true if the first entry has smaller ROW index than the second one
bool operator()(Triple_ firstTriple, Triple_ secondTriple)