dune-istl 3.0-git
colcompmatrix.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_COLCOMPMATRIX_HH
4#define DUNE_ISTL_COLCOMPMATRIX_HH
5#include "bcrsmatrix.hh"
6#include "bvector.hh"
7#include <dune/common/fmatrix.hh>
8#include <dune/common/fvector.hh>
9#include <dune/common/typetraits.hh>
10#include <dune/common/unused.hh>
11#include <limits>
12
13namespace Dune
14{
20 template<class M>
22 {
23 public:
25 typedef M Matrix;
28
34 : m_(m)
35 {}
36
39 {
40 return m_.begin();
41 }
44 {
45 return m_.end();
46 }
47 private:
48 const Matrix& m_;
49 };
50
58 template<class M, class S>
60 {
61 public:
63 typedef M Matrix;
65 typedef S RowIndexSet;
66
73 : m_(m), s_(s)
74 {}
75
76 const Matrix& matrix() const
77 {
78 return m_;
79 }
80
81 const RowIndexSet& rowIndexSet() const
82 {
83 return s_;
84 }
85
88 : public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
89 {
90 public:
91 const_iterator(typename Matrix::const_iterator firstRow,
92 typename RowIndexSet::const_iterator pos)
93 : firstRow_(firstRow), pos_(pos)
94 {}
95
96
97 const typename Matrix::row_type& dereference() const
98 {
99 return *(firstRow_+ *pos_);
100 }
101 bool equals(const const_iterator& o) const
102 {
103 return pos_==o.pos_;
104 }
106 {
107 ++pos_;
108 }
109 typename RowIndexSet::value_type index() const
110 {
111 return *pos_;
112 }
113
114 private:
116 typename Matrix::const_iterator firstRow_;
118 typename RowIndexSet::const_iterator pos_;
119 };
120
123 {
124 return const_iterator(m_.begin(), s_.begin());
125 }
128 {
129 return const_iterator(m_.begin(), s_.end());
130 }
131
132 private:
134 const Matrix& m_;
136 const RowIndexSet& s_;
137 };
138
143 template<class M>
145 {};
146
152 template<class M>
155
156 template<class M, class X, class TM, class TD, class T1>
158
159 template<class T, bool flag>
161
166 template<class B, class TA, int n, int m>
168 {
169 friend struct ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<B,n,m>,TA> >;
170
171 public:
175
176 typedef typename Matrix::size_type size_type;
177
182 explicit ColCompMatrix(const Matrix& mat);
183
185
187 virtual ~ColCompMatrix();
188
193 size_type N() const
194 {
195 return N_;
196 }
197
199 {
200 return Nnz_/n/m;
201 }
202
207 size_type M() const
208 {
209 return M_;
210 }
211
212 B* getValues() const
213 {
214 return values;
215 }
216
217 int* getRowIndex() const
218 {
219 return rowindex;
220 }
221
222 int* getColStart() const
223 {
224 return colstart;
225 }
226
229
235 virtual void setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs);
237 virtual void free();
238
240 virtual void setMatrix(const Matrix& mat);
241
242 public:
243 int N_, M_, Nnz_;
247 };
248
249 template<class T, class A, int n, int m>
251 {
252 template<class I, class S, class D>
254 public:
258 typedef typename Matrix::size_type size_type;
259
261
263
265
266 template<typename Iter>
267 void addRowNnz(const Iter& row) const;
268
269 template<typename Iter, typename Set>
270 void addRowNnz(const Iter& row, const Set& s) const;
271
272 void allocate();
273
274 template<typename Iter>
275 void countEntries(const Iter& row, const CIter& col) const;
276
277 void countEntries(size_type colidx) const;
278
279 void calcColstart() const;
280
281 template<typename Iter>
282 void copyValue(const Iter& row, const CIter& col) const;
283
284 void copyValue(const CIter& col, size_type rowindex, size_type colidx) const;
285
286 virtual void createMatrix() const;
287
288 protected:
289
290 void allocateMatrixStorage() const;
291
292 void allocateMarker();
293
297 };
298
299 template<class T, class A, int n, int m>
301 : mat(&mat_), cols(mat_.M()), marker(0)
302 {
303 mat->Nnz_=0;
304 }
305
306 template<class T, class A, int n, int m>
308 : mat(0), cols(0), marker(0)
309 {}
310
311 template<class T, class A, int n, int m>
313 {
314 if(marker)
315 delete[] marker;
316 }
317
318 template<class T, class A, int n, int m>
319 template<typename Iter>
320 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row) const
321 {
322 mat->Nnz_+=row->getsize();
323 }
324
325 template<class T, class A, int n, int m>
326 template<typename Iter, typename Map>
327 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row,
328 const Map& indices) const
329 {
330 typedef typename Iter::value_type::const_iterator RIter;
331 typedef typename Map::const_iterator MIter;
332 MIter siter =indices.begin();
333 for(RIter entry=row->begin(); entry!=row->end(); ++entry)
334 {
335 for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
336 if(siter==indices.end())
337 break;
338 if(*siter==entry.index())
339 // index is in subdomain
340 ++mat->Nnz_;
341 }
342 }
343
344 template<class T, class A, int n, int m>
346 {
347 allocateMatrixStorage();
348 allocateMarker();
349 }
350
351 template<class T, class A, int n, int m>
352 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage() const
353 {
354 mat->Nnz_*=n*m;
355 // initialize data
356 mat->values=new T[mat->Nnz_];
357 mat->rowindex=new int[mat->Nnz_];
358 mat->colstart=new int[cols+1];
359 }
360
361 template<class T, class A, int n, int m>
363 {
364 marker = new typename Matrix::size_type[cols];
365
366 for(size_type i=0; i < cols; ++i)
367 marker[i]=0;
368 }
369
370 template<class T, class A, int n, int m>
371 template<typename Iter>
373 {
375 countEntries(col.index());
376 }
377
378 template<class T, class A, int n, int m>
380 {
381 for(size_type i=0; i < m; ++i)
382 {
383 assert(colindex*m+i<cols);
384 marker[colindex*m+i]+=n;
385 }
386 }
387
388 template<class T, class A, int n, int m>
390 {
391 mat->colstart[0]=0;
392 for(size_type i=0; i < cols; ++i) {
393 assert(i<cols);
394 mat->colstart[i+1]=mat->colstart[i]+marker[i];
395 marker[i]=mat->colstart[i];
396 }
397 }
398
399 template<class T, class A, int n, int m>
400 template<typename Iter>
401 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const Iter& row, const CIter& col) const
402 {
403 copyValue(col, row.index(), col.index());
404 }
405
406 template<class T, class A, int n, int m>
408 {
409 for(size_type i=0; i<n; i++) {
410 for(size_type j=0; j<m; j++) {
411 assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
412 assert((int)marker[colindex*m+j]<mat->Nnz_);
413 mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
414 mat->values[marker[colindex*m+j]]=(*col)[i][j];
415 ++marker[colindex*m+j]; // index for next entry in column
416 }
417 }
418 }
419
420 template<class T, class A, int n, int m>
422 {
423 delete[] marker;
424 marker=0;
425 }
426
427 template<class F, class MRS>
429 {
430 typedef typename MRS::const_iterator Iter;
431 typedef typename std::iterator_traits<Iter>::value_type::const_iterator CIter;
432 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
433 initializer.addRowNnz(row);
434
435 initializer.allocate();
436
437 for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
438
439 for(CIter col=row->begin(); col != row->end(); ++col)
440 initializer.countEntries(row, col);
441 }
442
443 initializer.calcColstart();
444
445 for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
446 for(CIter col=row->begin(); col != row->end(); ++col) {
447 initializer.copyValue(row, col);
448 }
449
450 }
451 initializer.createMatrix();
452 }
453
454 template<class F, class M,class S>
456 {
458 typedef typename MRS::RowIndexSet SIS;
459 typedef typename SIS::const_iterator SIter;
460 typedef typename MRS::const_iterator Iter;
461 typedef typename std::iterator_traits<Iter>::value_type row_type;
462 typedef typename row_type::const_iterator CIter;
463
464 // Calculate upper Bound for nonzeros
465 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
466 initializer.addRowNnz(row, mrs.rowIndexSet());
467
468 initializer.allocate();
469
470 typedef typename MRS::Matrix::size_type size_type;
471
472 // A vector containing the corresponding indices in
473 // the to create submatrix.
474 // If an entry is the maximum of size_type then this index will not appear in
475 // the submatrix.
476 std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
477 std::numeric_limits<size_type>::max());
478 size_type s=0;
479 for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
480 subMatrixIndex[*index]=s++;
481
482 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
483 for(CIter col=row->begin(); col != row->end(); ++col) {
484 if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
485 // This column is in our subset (use submatrix column index)
486 initializer.countEntries(subMatrixIndex[col.index()]);
487 }
488
489 initializer.calcColstart();
490
491 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
492 for(CIter col=row->begin(); col != row->end(); ++col) {
493 if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
494 // This value is in our submatrix -> copy (use submatrix indices
495 initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
496 }
497 initializer.createMatrix();
498 }
499
500#ifndef DOXYGEN
501
502 template<class B, class TA, int n, int m>
503 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::ColCompMatrix()
504 : N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
505 {}
506
507 template<class B, class TA, int n, int m>
508 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
509 ::ColCompMatrix(const Matrix& mat)
510 : N_(n*mat.N()), M_(m*mat.M()), Nnz_(n*m*mat.nonzeroes())
511 {}
512
513 template<class B, class TA, int n, int m>
514 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
515 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const Matrix& mat)
516 {
517 if(N_+M_+Nnz_!=0)
518 free();
519 setMatrix(mat);
520 return *this;
521 }
522
523 template<class B, class TA, int n, int m>
524 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
525 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const ColCompMatrix& mat)
526 {
527 if(N_+M_+Nnz_!=0)
528 free();
529 N_=mat.N_;
530 M_=mat.M_;
531 Nnz_= mat.Nnz_;
532 if(M_>0) {
533 colstart=new int[M_+1];
534 for(int i=0; i<=M_; ++i)
535 colstart[i]=mat.colstart[i];
536 }
537
538 if(Nnz_>0) {
539 values = new B[Nnz_];
540 rowindex = new int[Nnz_];
541
542 for(int i=0; i<Nnz_; ++i)
543 values[i]=mat.values[i];
544
545 for(int i=0; i<Nnz_; ++i)
546 rowindex[i]=mat.rowindex[i];
547 }
548 return *this;
549 }
550
551 template<class B, class TA, int n, int m>
552 void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
553 ::setMatrix(const Matrix& mat)
554 {
555 N_=n*mat.N();
556 M_=m*mat.M();
557 ColCompMatrixInitializer<Matrix> initializer(*this);
558
559 copyToColCompMatrix(initializer, MatrixRowSet<Matrix>(mat));
560 }
561
562 template<class B, class TA, int n, int m>
563 void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
564 ::setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs)
565 {
566 if(N_+M_+Nnz_!=0)
567 free();
568 N_=mrs.size()*n;
569 M_=mrs.size()*m;
570 ColCompMatrixInitializer<Matrix> initializer(*this);
571
572 copyToColCompMatrix(initializer, MatrixRowSubset<Matrix,std::set<std::size_t> >(mat,mrs));
573 }
574
575 template<class B, class TA, int n, int m>
576 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::~ColCompMatrix()
577 {
578 if(N_+M_+Nnz_!=0)
579 free();
580 }
581
582 template<class B, class TA, int n, int m>
583 void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::free()
584 {
585 delete[] values;
586 delete[] rowindex;
587 delete[] colstart;
588 N_ = 0;
589 M_ = 0;
590 Nnz_ = 0;
591 }
592
593#endif // DOXYGEN
594
595}
596#endif
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Implementation of the BCRSMatrix class.
Col col
Definition matrixmatrix.hh:347
Matrix & mat
Definition matrixmatrix.hh:343
std::size_t countEntries(const BlockVector< FieldVector< T, i >, A > &vector)
Definition matrixmarket.hh:886
Definition basearray.hh:19
void copyToColCompMatrix(F &initializer, const MRS &mrs)
Definition colcompmatrix.hh:428
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:412
A::size_type size_type
The type for the index access and the size.
Definition bcrsmatrix.hh:446
Iterator access to matrix rows
Definition bcrsmatrix.hh:527
Definition bvector.hh:660
Provides access to an iterator over all matrix rows.
Definition colcompmatrix.hh:22
Matrix::ConstRowIterator const_iterator
The matrix row iterator type.
Definition colcompmatrix.hh:27
const_iterator begin() const
Get the row iterator at the first row.
Definition colcompmatrix.hh:38
M Matrix
The type of the matrix.
Definition colcompmatrix.hh:25
MatrixRowSet(const Matrix &m)
Construct an row set over all matrix rows.
Definition colcompmatrix.hh:33
const_iterator end() const
Get the row iterator at the end of all rows.
Definition colcompmatrix.hh:43
Provides access to an iterator over an arbitrary subset of matrix rows.
Definition colcompmatrix.hh:60
const RowIndexSet & rowIndexSet() const
Definition colcompmatrix.hh:81
const Matrix & matrix() const
Definition colcompmatrix.hh:76
S RowIndexSet
the type of the set of valid row indices.
Definition colcompmatrix.hh:65
const_iterator begin() const
Get the row iterator at the first row.
Definition colcompmatrix.hh:122
M Matrix
the type of the matrix class.
Definition colcompmatrix.hh:63
MatrixRowSubset(const Matrix &m, const RowIndexSet &s)
Construct an row set over all matrix rows.
Definition colcompmatrix.hh:72
const_iterator end() const
Get the row iterator at the end of all rows.
Definition colcompmatrix.hh:127
The matrix row iterator type.
Definition colcompmatrix.hh:89
bool equals(const const_iterator &o) const
Definition colcompmatrix.hh:101
const Matrix::row_type & dereference() const
Definition colcompmatrix.hh:97
const_iterator(typename Matrix::const_iterator firstRow, typename RowIndexSet::const_iterator pos)
Definition colcompmatrix.hh:91
RowIndexSet::value_type index() const
Definition colcompmatrix.hh:109
void increment()
Definition colcompmatrix.hh:105
Utility class for converting an ISTL Matrix into a column-compressed matrix.
Definition colcompmatrix.hh:145
Inititializer for the ColCompMatrix as needed by OverlappingSchwarz.
Definition colcompmatrix.hh:154
Sequential overlapping Schwarz preconditioner.
Definition overlappingschwarz.hh:742
Definition overlappingschwarz.hh:683
BCRSMatrix< FieldMatrix< B, n, m >, TA > Matrix
The type of the matrix to convert.
Definition colcompmatrix.hh:173
B * getValues() const
Definition colcompmatrix.hh:212
int * getRowIndex() const
Definition colcompmatrix.hh:217
Matrix::size_type size_type
Definition colcompmatrix.hh:176
int * getColStart() const
Definition colcompmatrix.hh:222
size_type nnz() const
Definition colcompmatrix.hh:198
ColCompMatrix(const Matrix &mat)
Constructor that initializes the data.
virtual void setMatrix(const Matrix &mat, const std::set< std::size_t > &mrs)
Initialize data from a given set of matrix rows and columns.
ColCompMatrix & operator=(const ColCompMatrix &mat)
size_type N() const
Get the number of rows.
Definition colcompmatrix.hh:193
virtual void setMatrix(const Matrix &mat)
Initialize data from given matrix.
size_type M() const
Get the number of columns.
Definition colcompmatrix.hh:207
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A > Matrix
Definition colcompmatrix.hh:255
Dune::ColCompMatrix< Matrix > ColCompMatrix
Definition colcompmatrix.hh:256
Matrix::size_type size_type
Definition colcompmatrix.hh:258
Matrix::row_type::const_iterator CIter
Definition colcompmatrix.hh:257
ConstIterator class for sequential access.
Definition matrix.hh:398
A generic dynamic dense matrix.
Definition matrix.hh:555
Definition matrixutils.hh:25
Initializer for SuperLU Matrices representing the subdomains.
Definition overlappingschwarz.hh:43