Flow123d  master-f44eb46
matrix_coo.hpp
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file matrix_coo.hpp
15  * @brief
16  */
17 
18 #ifndef la_matrix_coo_h
19 #define la_matrix_coo_h
20 
21 #define ALLOWED_MEM_OVERHEAD 2.e+6 // entries, this corresponds to 32 MB of
22  // overhead in <int, int, double> matrixCoo
23 //------------------------------------------------------------------------------
24 #include <vector>
25 #include <cassert>
26 #include <cmath>
27 #include <algorithm>
28 #include <iomanip>
29 #include <numeric>
30 
31 //------------------------------------------------------------------------------
32 namespace la{
33  template< typename INDT, typename VALT >
34  class MatrixCoo;
35 }
36 
37 //------------------------------------------------------------------------------
38 /** \brief Temporary storage for a sparse matrix
39  * \details This object stores and handles matrix in coordinate format, i.e.,
40  * <row_index,col_index,value>-combinations for each matrix entry of a sparse matrix.
41  * Functionalities are:
42  * - preparing assembly of elements - allocating proper space to all elements
43  * in unassembled form
44  * - inserting an element - it is user's responsibility not to insert zeros
45  * - finishing assembly of elements - this is the core routine which
46  * - sorts triples, merges entries with the same pair of indices,
47  * truncates the array to assembled size
48  * - 'zeroing' a row and setting its diagonal entry to a given value
49  * - query of number of entries - this counts also zeros if they were inserted
50  * - get inside/outide-a-range counts of elements - for PETSc preallocation
51  * - filling the arrays for a row-compressed data storage for external
52  * solver libraries
53  * - computing inf norm
54  * - cleanup and writing
55  */
56 template< typename INDT, typename VALT >
58 {
59 public:
60  typedef INDT IndexType;
61  typedef VALT ValueType;
62 
63 public:
64  MatrixCoo( unsigned nRows = 0, unsigned nCols = 0 )
65  : nRows_( nRows ),
66  nCols_( nCols ),
67  sorted_ ( false )
68  {
69  sortedSize_ = 0;
70  }
71 
72  ~MatrixCoo( ) {
73  matrixCooVec_.clear();
74  }
75 
76  //! Return number of non-zeros
77  IndexType nnz() const { return matrixCooVec_.size(); }
78 
79  //! Prepare matrix for assembly - allocate buffer for UNASSEMBLED entries
80  //! serial copies of element matrices to memory
81  //! length = numElements * size( elementMatrix )
82  void prepareAssembly( const IndexType length )
83  {
84  matrixCooVec_.reserve( length );
85 
86  return;
87  }
88 
89  //! Insert a value to unassembled vector
90  void insert( const IndexType i, const IndexType j, const ValueType value )
91  {
92  // combine index pair
93  IndexPair_ ip( i, j );
94 
95  // combine triples
96  Triple_ trip( ip, value );
97 
98  // insert it to buffer in any case
99  matrixCooVec_.push_back( trip );
100 
101  sorted_ = false;
102 
103  // resort if number of new unsorted entries excess a critical value
104  if ( (matrixCooVec_.size() - sortedSize_) > ALLOWED_MEM_OVERHEAD ) {
105  this -> finishAssembly();
106  }
107 
108  return;
109  }
110  //! Finalize matrix assembly - build the assembled input and truncate matrix
112  {
113  // sort entries to make sure
114  std::sort( matrixCooVec_.begin( ), matrixCooVec_.end( ), TripleLessThan_( ) );
115 
116  // combine the sorted matrix into one vector
117  typename MatrixCooVec_::const_iterator itGlobale = matrixCooVec_.end();
118 
119  typename MatrixCooVec_::const_iterator itSelectb = matrixCooVec_.begin();
120  typename MatrixCooVec_::const_iterator itSelecte = itSelectb;
121 
122  Triple_ tmpTriple;
123 
124  typename MatrixCooVec_::iterator itStore = matrixCooVec_.begin();
125 
126  // for each pair of indices
127  while (itSelectb != itGlobale) {
128 
129  // find any change in indices
130  itSelecte = std::upper_bound( itSelectb, itGlobale, *itSelectb, TripleLessThan_( ) );
131 
132  // sum values in the range
133  tmpTriple = combineTriples_ ( itSelectb, itSelecte );
134 
135  // store new
136  *itStore = tmpTriple;
137  itStore++;
138 
139  // shift ranges
140  itSelectb = itSelecte;
141  }
142 
143  // determine number of nonzeros in assembled matrix
144  unsigned nnz = std::distance( matrixCooVec_.begin( ), itStore);
145 
146  // truncate the vector from the end
147  matrixCooVec_.resize( nnz );
148 
149  // mark the amount of sorted entries
150  sortedSize_ = nnz;
151 
152  sorted_ = true;
153 
154  return;
155  }
156 
157  //! Determine number of nonzeros within a block [low,high) - used for PETSc preallocation
159  std::vector<IndexType> & nnzRowsIn, std::vector<IndexType> & nnzRowsOut)
160  {
161  // sort entries to make sure
162  if ( not sorted_ ) {
163  this -> finishAssembly( );
164  }
165 
166  // run through all entries
167  typename MatrixCooVec_::const_iterator iter = matrixCooVec_.begin();
168  typename MatrixCooVec_::const_iterator itere = matrixCooVec_.end();
169 
170  for ( ; iter != itere; ++iter ) {
171 
172  IndexType row = (*iter).first.first;
173  IndexType col = (*iter).first.second;
174 
175  // check that it belongs to the interval
176  if ( row >= low and row < high ) {
177 
178  IndexType rowLoc = row - low;
179  // it is on my row, now where to put it
180 
181  if ( col >= low and col < high ) {
182  // it is inside the diagonal block
183  nnzRowsIn[rowLoc]++;
184  }
185  else {
186  // it is outside the diagonal block
187  nnzRowsOut[rowLoc]++;
188  }
189  }
190  }
191  }
192 
193  //! Get specific entry
194  void getEntry( IndexType index, IndexType & row, IndexType & col, ValueType & val )
195  {
196  row = matrixCooVec_[ index ].first.first;
197  col = matrixCooVec_[ index ].first.second;
198  val = matrixCooVec_[ index ].second;
199  }
200 
201  //! Zero a row and set diagonal to scalar
202  void zeroRow( const IndexType row, const ValueType factor )
203  {
204 
205  // sort entries to make sure
206  if ( not sorted_ ) {
207  this -> finishAssembly( );
208  }
209 
210  // erasing iterator
211  typename MatrixCooVec_::iterator tripEraseIter;
212  typename MatrixCooVec_::iterator tripEraseIterE;
213 
214  // find constrained row in matrix and set iterator there
215  tripEraseIter = std::lower_bound ( matrixCooVec_.begin(), matrixCooVec_.end(),
216  row, TripleRowIndexLessThan_( ) );
217 
218  // find last entry at the row
219  tripEraseIterE = std::upper_bound( tripEraseIter, matrixCooVec_.end(),
220  *tripEraseIter, TripleRowLessThan_( ) );
221 
222  // find diagonal
223  typename MatrixCooVec_::iterator tripEraseIterDiag;
224  tripEraseIterDiag = std::find_if( tripEraseIter, tripEraseIterE, TripleDiagEntry_( ) );
225 
226  // set row entries to zero
227  for ( ; tripEraseIter != tripEraseIterE; ++tripEraseIter ) {
228  (*tripEraseIter).second = static_cast< ValueType >( 0.0 );
229  }
230 
231  // put value from constraint on diagonal
232  Triple_ tmpTrip = std::make_pair( std::make_pair ( row, row ), factor );
233 
234  if ( tripEraseIterDiag != tripEraseIterE ) {
235  // diagonal entry already present
236  *tripEraseIterDiag = tmpTrip;
237  }
238  else {
239  // nothing on diagonal, add it
240  matrixCooVec_.push_back( tmpTrip );
241  sorted_ = false;
242  }
243  }
244 
245  //! Fill arrays with column indices, row pointers and values for
246  //! row-compressed data storage
248  std::vector<IndexType> & rowPtrs,
249  std::vector<ValueType> & values )
250  {
251 
252  // sort entries to make sure
253  if ( not sorted_ ) {
254  this -> finishAssembly( );
255  }
256 
257  // iterate over all elements and fill the nonzero,
258  // column index, and row pointer arrys
259  typename MatrixCooVec_::const_iterator tvIter = matrixCooVec_.begin();
260  typename MatrixCooVec_::const_iterator tvEnd = matrixCooVec_.end();
261 
262  // reserve space
263  colIndices.reserve(matrixCooVec_.size());
264  values.reserve(matrixCooVec_.size());
265 
266  unsigned nr = this -> deduceNumRows_( );
267  rowPtrs.reserve( nr + 1 );
268 
269  int previousRow = -1;
270  IndexType index = 0;
271 
272  for ( ; tvIter != tvEnd; ++tvIter ) {
273  const IndexType row = (*tvIter).first.first;
274  const IndexType col = (*tvIter).first.second;
275 
276  values.push_back( tvIter -> second );
277  colIndices.push_back( col );
278  if ( row == previousRow ) index++;
279  else rowPtrs.push_back( index++ );
280 
281  previousRow = row;
282  }
283 
284  rowPtrs.push_back( values.size() );
285 
286  return;
287  }
288 
289  //! Fill independent arrays with column indices, row pointers and values
291  std::vector<IndexType> & colIndices,
292  std::vector<ValueType> & values,
293  bool onlyUpperTriangle = false )
294  {
295 
296  // sort entries to make sure
297  if ( not sorted_ ) {
298  this -> finishAssembly( );
299  }
300 
301  // iterate over all elements and fill the nonzero,
302  // column index, and row index arrys
303  typename MatrixCooVec_::const_iterator tvIter = matrixCooVec_.begin();
304  typename MatrixCooVec_::const_iterator tvEnd = matrixCooVec_.end();
305 
306  // reserve space
307  colIndices.reserve(matrixCooVec_.size());
308  rowIndices.reserve(matrixCooVec_.size());
309  values.reserve(matrixCooVec_.size());
310 
311  for ( ; tvIter != tvEnd; ++tvIter ) {
312  const IndexType row = (*tvIter).first.first ;
313  const IndexType col = (*tvIter).first.second ;
314 
315  if ( not onlyUpperTriangle or col >= row ) {
316  values.push_back( (*tvIter).second );
317  rowIndices.push_back( row );
318  colIndices.push_back( col );
319  }
320  }
321  return;
322  }
323 
324  //! destroy all contents
325  void clear()
326  {
327  matrixCooVec_.clear();
328  // force deallocation of the vector
329  MatrixCooVec_().swap( matrixCooVec_ );
330  return;
331  }
332 
333  //! erase all contents
334  void erase()
335  {
336  // erasing iterator
337  typename MatrixCooVec_::iterator tripEraseIter = matrixCooVec_.begin( );
338  typename MatrixCooVec_::iterator tripEraseIterE = matrixCooVec_.end( );
339 
340  // loop to erase values
341  for ( ; tripEraseIter != tripEraseIterE; ++tripEraseIter ) {
342  (*tripEraseIter).second = static_cast< ValueType >(0.0);
343  }
344 
345  return;
346  }
347 
348  //! debug-output
349  std::ostream & write( std::ostream & out )
350  {
351  // sort entries to make sure
352  if ( not sorted_ ) {
353  this -> finishAssembly( );
354  }
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;
360 
361  out << row << " " << col << " "
362  << std::setprecision(12)
363  << (*tvIter).second << std::endl;
364 
365  }
366  return out;
367  }
368 
369  //----------------------------------------------------------------------------
370  //! Compute maximal row sum (infinity-norm) of a matrix in coordinate format.
371  double infNorm()
372  {
373  // sort entries to make sure
374  if ( not sorted_ ) {
375  this -> finishAssembly( );
376  }
377  unsigned nr = this -> deduceNumRows_( );
378  std::vector<double> rowSums( nr, 0.0 );
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 );
384  }
385 
386  return *std::max_element( rowSums.begin(), rowSums.end() );
387  }
388 
389  //--------------------------------------------------------------------------
390  //! Return Frobenius norm of matrix \f$||a_{ij}||_F=\sqrt{\sum_i\sum_j |a_{ij}|^2}\f$
392  {
393  // sort entries to make sure
394  if ( not sorted_ ) {
395  this -> finishAssembly( );
396  }
397 
398  ValueType sumOfSquares = 0.0;
399 
400  typename MatrixCooVec_::const_iterator iter = matrixCooVec_.begin();
401  typename MatrixCooVec_::const_iterator last = matrixCooVec_.end();
402  for ( ; iter != last; ++iter ) {
403 
404  ValueType val = (*iter).second;
405  sumOfSquares += val * val;
406  }
407 
408  return std::sqrt( sumOfSquares );
409  }
410 
411  //--------------------------------------------------------------------------
412  //! Get scalar suitable to be put on diagonal while fixing BC
413  //! This is the mean value of extremal entries on diagonal - the aim of this guess is to
414  //! be within the spectra of the original matrix - this value immediately corresponds to an
415  //! eigenvalue
416  double getDiagScalar( )
417  {
418 
419  ValueType val;
420  // initialize maximum and minimum
421  if ( matrixCooVec_.size() > 0 ) {
422  val = std::fabs( matrixCooVec_[0].second );
423  }
424  else {
425  val = 0.0;
426  }
427  ValueType maxVal = val;
428  ValueType minVal = val;
429  IndexType row, col;
430  for ( IndexType i = 1; i < static_cast<IndexType>( matrixCooVec_.size() ); i++ ) {
431 
432  row = matrixCooVec_[i].first.first;
433  col = matrixCooVec_[i].first.second;
434 
435  if ( row == col ) {
436  val = std::fabs( matrixCooVec_[i].second );
437 
438  // update maximum
439  if ( val > maxVal ) {
440  maxVal = val;
441  }
442 
443  // update minimum
444  if ( val < minVal ) {
445  minVal = val;
446  }
447  }
448  }
449 
450  // recommend value in between extremal absolute values along diagonal
451  ValueType scalar = ( minVal + maxVal ) / 2.0;
452 
453  return scalar;
454 
455  }
456 
457 private:
458  typedef std::pair<IndexType, IndexType> IndexPair_; //!< pair of indices
459  typedef std::pair<IndexPair_,ValueType> Triple_; //!< triple
460  typedef std::vector<Triple_> MatrixCooVec_;//!< vector of triples
461  MatrixCooVec_ matrixCooVec_;//!< matrix in coordinate format <i,j, values>
462 
463  IndexType nRows_; //!< number of rows of matrix
464  IndexType nCols_; //!< number of columns of matrix
465 
466  unsigned sortedSize_; //!< number of sorted entries
467 
468  bool sorted_; //!< state variable - is matrix sorted?
469 
470  //--------------------------------------------------------------------------
471  // private functions
472 
473  //! deduce number of rows - if it was given by user, use it, otherwise,
474  //! deduce it from the index of the last entry
476  {
477  IndexType nr;
478 
479  if ( nRows_ > 0 ) {
480  // number of rows was given by user
481  nr = nRows_;
482  }
483  else {
484  // deduce number of rows from sorted matrix
485  nr = (*matrixCooVec_.rbegin()).first.first + 1;
486  }
487 
488  return nr;
489  }
490 
491 
492  //! returns true if the index pair of the first entry is smaller than of the second entry
493  struct TripleLessThan_ : std::binary_function< Triple_, Triple_, bool >
494  {
495  bool operator() ( Triple_ firstTriple, Triple_ secondTriple) {
496 
497  IndexPair_ firstTripleIndxPair = firstTriple.first;
498  IndexPair_ secondTripleIndxPair = secondTriple.first;
499 
500  return (firstTripleIndxPair < secondTripleIndxPair);
501  }
502  };
503 
504  //! returns true if the first entry has smaller ROW index than the second one
505  struct TripleRowLessThan_ : std::binary_function< Triple_, Triple_, bool >
506  {
507  bool operator() ( Triple_ firstTriple, Triple_ secondTriple) {
508 
509  IndexType firstRowIndex = firstTriple.first.first;
510  IndexType secondRowIndex = secondTriple.first.first;
511 
512  return (firstRowIndex < secondRowIndex);
513  }
514  };
515 
516  //! returns true if the row of an entry is smaller than a prescribed index
518  {
519  bool operator() ( Triple_ triple, IndexType row ) {
520 
521  IndexType tripleRowIndex = triple.first.first;
522 
523  return (tripleRowIndex < row);
524  }
525  };
526 
527  //! returns true if the entry is on diagonal, i.e. row == column
529  {
530  bool operator() ( Triple_ triple ) {
531 
532  IndexType tripleRowIndex = triple.first.first;
533  IndexType tripleColIndex = triple.first.second;
534 
535  return (tripleRowIndex == tripleColIndex);
536  }
537  };
538 
539  //! merge range of triples by addition
540  Triple_ combineTriples_ ( typename MatrixCooVec_::const_iterator itb, typename MatrixCooVec_::const_iterator ite ) {
541 
542  IndexPair_ tmpPair = (*itb).first;
543  ValueType val = 0.0;
544 
545  for ( ; itb != ite; ++itb) {
546  val += (*itb).second;
547  }
548 
549  Triple_ combinedTriple = std::make_pair(tmpPair, val);
550 
551  return combinedTriple;
552 
553  }
554 };
555 
556 #endif
Temporary storage for a sparse matrix.
Definition: matrix_coo.hpp:58
unsigned sortedSize_
number of sorted entries
Definition: matrix_coo.hpp:466
IndexType deduceNumRows_()
Definition: matrix_coo.hpp:475
void clear()
destroy all contents
Definition: matrix_coo.hpp:325
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.
Definition: matrix_coo.hpp:158
bool sorted_
state variable - is matrix sorted?
Definition: matrix_coo.hpp:468
std::vector< Triple_ > MatrixCooVec_
vector of triples
Definition: matrix_coo.hpp:460
void prepareAssembly(const IndexType length)
Definition: matrix_coo.hpp:82
void zeroRow(const IndexType row, const ValueType factor)
Zero a row and set diagonal to scalar.
Definition: matrix_coo.hpp:202
IndexType nRows_
number of rows of matrix
Definition: matrix_coo.hpp:463
void insert(const IndexType i, const IndexType j, const ValueType value)
Insert a value to unassembled vector.
Definition: matrix_coo.hpp:90
std::ostream & write(std::ostream &out)
debug-output
Definition: matrix_coo.hpp:349
void erase()
erase all contents
Definition: matrix_coo.hpp:334
IndexType nnz() const
Return number of non-zeros.
Definition: matrix_coo.hpp:77
IndexType nCols_
number of columns of matrix
Definition: matrix_coo.hpp:464
ValueType frobeniusNorm()
Return Frobenius norm of matrix .
Definition: matrix_coo.hpp:391
double infNorm()
Compute maximal row sum (infinity-norm) of a matrix in coordinate format.
Definition: matrix_coo.hpp:371
MatrixCooVec_ matrixCooVec_
matrix in coordinate format <i,j, values>
Definition: matrix_coo.hpp:461
Triple_ combineTriples_(typename MatrixCooVec_::const_iterator itb, typename MatrixCooVec_::const_iterator ite)
merge range of triples by addition
Definition: matrix_coo.hpp:540
std::pair< IndexType, IndexType > IndexPair_
pair of indices
Definition: matrix_coo.hpp:458
std::pair< IndexPair_, ValueType > Triple_
triple
Definition: matrix_coo.hpp:459
double getDiagScalar()
Definition: matrix_coo.hpp:416
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.
Definition: matrix_coo.hpp:290
MatrixCoo(unsigned nRows=0, unsigned nCols=0)
Definition: matrix_coo.hpp:64
void getEntry(IndexType index, IndexType &row, IndexType &col, ValueType &val)
Get specific entry.
Definition: matrix_coo.hpp:194
void setCompressedRowArrays(std::vector< IndexType > &colIndices, std::vector< IndexType > &rowPtrs, std::vector< ValueType > &values)
Definition: matrix_coo.hpp:247
void finishAssembly()
Finalize matrix assembly - build the assembled input and truncate matrix.
Definition: matrix_coo.hpp:111
static constexpr bool value
Definition: json.hpp:87
#define ALLOWED_MEM_OVERHEAD
Definition: matrix_coo.hpp:21
returns true if the entry is on diagonal, i.e. row == column
Definition: matrix_coo.hpp:529
bool operator()(Triple_ triple)
Definition: matrix_coo.hpp:530
returns true if the index pair of the first entry is smaller than of the second entry
Definition: matrix_coo.hpp:494
bool operator()(Triple_ firstTriple, Triple_ secondTriple)
Definition: matrix_coo.hpp:495
returns true if the row of an entry is smaller than a prescribed index
Definition: matrix_coo.hpp:518
bool operator()(Triple_ triple, IndexType row)
Definition: matrix_coo.hpp:519
returns true if the first entry has smaller ROW index than the second one
Definition: matrix_coo.hpp:506
bool operator()(Triple_ firstTriple, Triple_ secondTriple)
Definition: matrix_coo.hpp:507