Dune Core Modules (2.4.1)

colcompmatrix.hh
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"
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 }
105 void increment()
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>
154 {};
155
156 template<class M, class X, class TM, class TD, class T1>
158
159 template<class T>
160 struct SeqOverlappingSchwarzAssembler;
161
166 template<class B, class TA, int n, int m>
168 {
169 template<class M, class X, class TM, class TD, class T1>
170 friend class SeqOverlappingSchwarz;
171 friend struct ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<B,n,m>,TA> >;
172
173 public:
176 friend struct SeqOverlappingSchwarzAssembler<ColCompMatrix<Matrix> >;
177
178 typedef typename Matrix::size_type size_type;
179
184 explicit ColCompMatrix(const Matrix& mat);
185
187
189 virtual ~ColCompMatrix();
190
195 size_type N() const
196 {
197 return N_;
198 }
199
200 size_type nnz() const
201 {
202 return Nnz_/n/m;
203 }
204
209 size_type M() const
210 {
211 return M_;
212 }
213
214 B* getValues() const
215 {
216 return values;
217 }
218
219 int* getRowIndex() const
220 {
221 return rowindex;
222 }
223
224 int* getColStart() const
225 {
226 return colstart;
227 }
228
229 ColCompMatrix& operator=(const Matrix& mat);
230 ColCompMatrix& operator=(const ColCompMatrix& mat);
231
238 virtual void setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs);
240 virtual void free();
241
243 virtual void setMatrix(const Matrix& mat);
244
245 public:
246 int N_, M_, Nnz_;
247 B* values;
248 int* rowindex;
249 int* colstart;
250 };
251
252 template<class T, class A, int n, int m>
254 {
255 template<class I, class S, class D>
257 public:
260 typedef typename Matrix::row_type::const_iterator CIter;
261 typedef typename Matrix::size_type size_type;
262
264
266
268
269 template<typename Iter>
270 void addRowNnz(const Iter& row) const;
271
272 template<typename Iter, typename Set>
273 void addRowNnz(const Iter& row, const Set& s) const;
274
275 void allocate();
276
277 template<typename Iter>
278 void countEntries(const Iter& row, const CIter& col) const;
279
280 void countEntries(size_type colidx) const;
281
282 void calcColstart() const;
283
284 template<typename Iter>
285 void copyValue(const Iter& row, const CIter& col) const;
286
287 void copyValue(const CIter& col, size_type rowindex, size_type colidx) const;
288
289 virtual void createMatrix() const;
290
291 protected:
292
293 void allocateMatrixStorage() const;
294
295 void allocateMarker();
296
297 ColCompMatrix* mat;
298 size_type cols;
299 mutable size_type *marker;
300 };
301
302 template<class T, class A, int n, int m>
303 ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::ColCompMatrixInitializer(ColCompMatrix& mat_)
304 : mat(&mat_), cols(mat_.N()), marker(0)
305 {
306 mat->Nnz_=0;
307 }
308
309 template<class T, class A, int n, int m>
310 ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::ColCompMatrixInitializer()
311 : mat(0), cols(0), marker(0)
312 {}
313
314 template<class T, class A, int n, int m>
315 ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::~ColCompMatrixInitializer()
316 {
317 if(marker)
318 delete[] marker;
319 }
320
321 template<class T, class A, int n, int m>
322 template<typename Iter>
323 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row) const
324 {
325 mat->Nnz_+=row->getsize();
326 }
327
328 template<class T, class A, int n, int m>
329 template<typename Iter, typename Map>
330 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row,
331 const Map& indices) const
332 {
333 typedef typename Iter::value_type::const_iterator RIter;
334 typedef typename Map::const_iterator MIter;
335 MIter siter =indices.begin();
336 for(RIter entry=row->begin(); entry!=row->end(); ++entry)
337 {
338 for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
339 if(siter==indices.end())
340 break;
341 if(*siter==entry.index())
342 // index is in subdomain
343 ++mat->Nnz_;
344 }
345 }
346
347 template<class T, class A, int n, int m>
348 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocate()
349 {
350 allocateMatrixStorage();
351 allocateMarker();
352 }
353
354 template<class T, class A, int n, int m>
355 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage() const
356 {
357 mat->Nnz_*=n*m;
358 // initialize data
359 mat->values=new T[mat->Nnz_];
360 mat->rowindex=new int[mat->Nnz_];
361 mat->colstart=new int[cols+1];
362 }
363
364 template<class T, class A, int n, int m>
365 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMarker()
366 {
367 marker = new typename Matrix::size_type[cols];
368
369 for(size_type i=0; i < cols; ++i)
370 marker[i]=0;
371 }
372
373 template<class T, class A, int n, int m>
374 template<typename Iter>
375 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(const Iter& row, const CIter& col) const
376 {
378 countEntries(col.index());
379 }
380
381 template<class T, class A, int n, int m>
382 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(size_type colindex) const
383 {
384 for(size_type i=0; i < m; ++i)
385 {
386 assert(colindex*m+i<cols);
387 marker[colindex*m+i]+=n;
388 }
389 }
390
391 template<class T, class A, int n, int m>
392 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::calcColstart() const
393 {
394 mat->colstart[0]=0;
395 for(size_type i=0; i < cols; ++i) {
396 assert(i<cols);
397 mat->colstart[i+1]=mat->colstart[i]+marker[i];
398 marker[i]=mat->colstart[i];
399 }
400 }
401
402 template<class T, class A, int n, int m>
403 template<typename Iter>
404 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const Iter& row, const CIter& col) const
405 {
406 copyValue(col, row.index(), col.index());
407 }
408
409 template<class T, class A, int n, int m>
410 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const CIter& col, size_type rowindex, size_type colindex) const
411 {
412 for(size_type i=0; i<n; i++) {
413 for(size_type j=0; j<m; j++) {
414 assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
415 assert((int)marker[colindex*m+j]<mat->Nnz_);
416 mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
417 mat->values[marker[colindex*m+j]]=(*col)[i][j];
418 ++marker[colindex*m+j]; // index for next entry in column
419 }
420 }
421 }
422
423 template<class T, class A, int n, int m>
424 void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::createMatrix() const
425 {
426 delete[] marker;
427 marker=0;
428 }
429
430 template<class F, class MRS>
431 void copyToColCompMatrix(F& initializer, const MRS& mrs)
432 {
433 typedef typename MRS::const_iterator Iter;
434 typedef typename std::iterator_traits<Iter>::value_type::const_iterator CIter;
435 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
436 initializer.addRowNnz(row);
437
438 initializer.allocate();
439
440 for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
441
442 for(CIter col=row->begin(); col != row->end(); ++col)
443 initializer.countEntries(row, col);
444 }
445
446 initializer.calcColstart();
447
448 for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
449 for(CIter col=row->begin(); col != row->end(); ++col) {
450 initializer.copyValue(row, col);
451 }
452
453 }
454 initializer.createMatrix();
455 }
456
457 template<class F, class M,class S>
458 void copyToColCompMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
459 {
460 typedef MatrixRowSubset<M,S> MRS;
461 typedef typename MRS::RowIndexSet SIS;
462 typedef typename SIS::const_iterator SIter;
463 typedef typename MRS::const_iterator Iter;
464 typedef typename std::iterator_traits<Iter>::value_type row_type;
465 typedef typename row_type::const_iterator CIter;
466
467 // Calculate upper Bound for nonzeros
468 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
469 initializer.addRowNnz(row, mrs.rowIndexSet());
470
471 initializer.allocate();
472
473 typedef typename MRS::Matrix::size_type size_type;
474
475 // A vector containing the corresponding indices in
476 // the to create submatrix.
477 // If an entry is the maximum of size_type then this index will not appear in
478 // the submatrix.
479 std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
480 std::numeric_limits<size_type>::max());
481 size_type s=0;
482 for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
483 subMatrixIndex[*index]=s++;
484
485 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
486 for(CIter col=row->begin(); col != row->end(); ++col) {
487 if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
488 // This column is in our subset (use submatrix column index)
489 initializer.countEntries(subMatrixIndex[col.index()]);
490 }
491
492 initializer.calcColstart();
493
494 for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
495 for(CIter col=row->begin(); col != row->end(); ++col) {
496 if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
497 // This value is in our submatrix -> copy (use submatrix indices
498 initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
499 }
500 initializer.createMatrix();
501 }
502
503#ifndef DOXYGEN
504
505 template<class B, class TA, int n, int m>
506 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::ColCompMatrix()
507 : N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
508 {}
509
510 template<class B, class TA, int n, int m>
511 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
512 ::ColCompMatrix(const Matrix& mat)
513 : N_(n*mat.N()), M_(m*mat.M()), Nnz_(n*m*mat.nonzeroes())
514 {}
515
516 template<class B, class TA, int n, int m>
517 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
518 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const Matrix& mat)
519 {
520 if(N_+M_+Nnz_!=0)
521 free();
522 setMatrix(mat);
523 return *this;
524 }
525
526 template<class B, class TA, int n, int m>
527 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
528 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const ColCompMatrix& mat)
529 {
530 if(N_+M_+Nnz_!=0)
531 free();
532 N_=mat.N_;
533 M_=mat.M_;
534 Nnz_= mat.Nnz_;
535 if(M_>0) {
536 colstart=new int[M_+1];
537 for(int i=0; i<=M_; ++i)
538 colstart[i]=mat.colstart[i];
539 }
540
541 if(Nnz_>0) {
542 values = new B[Nnz_];
543 rowindex = new int[Nnz_];
544
545 for(int i=0; i<Nnz_; ++i)
546 values[i]=mat.values[i];
547
548 for(int i=0; i<Nnz_; ++i)
549 rowindex[i]=mat.rowindex[i];
550 }
551 return *this;
552 }
553
554 template<class B, class TA, int n, int m>
555 void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
556 ::setMatrix(const Matrix& mat)
557 {
558 N_=n*mat.N();
559 M_=m*mat.M();
560 ColCompMatrixInitializer<Matrix> initializer(*this);
561
562 copyToColCompMatrix(initializer, MatrixRowSet<Matrix>(mat));
563 }
564
565 template<class B, class TA, int n, int m>
566 void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
567 ::setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs)
568 {
569 if(N_+M_+Nnz_!=0)
570 free();
571 N_=mrs.size()*n;
572 M_=mrs.size()*m;
573 ColCompMatrixInitializer<Matrix> initializer(*this);
574
575 copyToColCompMatrix(initializer, MatrixRowSubset<Matrix,std::set<std::size_t> >(mat,mrs));
576 }
577
578 template<class B, class TA, int n, int m>
579 ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::~ColCompMatrix()
580 {
581 if(N_+M_+Nnz_!=0)
582 free();
583 }
584
585 template<class B, class TA, int n, int m>
586 void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::free()
587 {
588 delete[] values;
589 delete[] rowindex;
590 delete[] colstart;
591 N_ = 0;
592 M_ = 0;
593 Nnz_ = 0;
594 }
595
596#endif // DOXYGEN
597
598}
599#endif
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:447
Definition: bvector.hh:584
BCRSMatrix< FieldMatrix< B, n, m >, TA > Matrix
The type of the matrix to convert.
Definition: colcompmatrix.hh:175
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.
size_type N() const
Get the number of rows.
Definition: colcompmatrix.hh:195
virtual void setMatrix(const Matrix &mat)
Initialize data from given matrix.
size_type M() const
Get the number of columns.
Definition: colcompmatrix.hh:209
A dense n x m matrix.
Definition: fmatrix.hh:67
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:141
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
The matrix row iterator type.
Definition: colcompmatrix.hh:89
Provides access to an iterator over an arbitrary subset of matrix rows.
Definition: colcompmatrix.hh:60
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
A generic dynamic dense matrix.
Definition: matrix.hh:25
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:41
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
size_type N() const
Return the number of rows.
Definition: matrix.hh:160
Initializer for SuperLU Matrices representing the subdomains.
Definition: overlappingschwarz.hh:43
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
ConstIterator class for sequential access.
Definition: vbvector.hh:647
RealIterator< const B > const_iterator
iterator class for sequential access
Definition: basearray.hh:195
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignment.hh:10
Inititializer for the ColCompMatrix as needed by OverlappingSchwarz.
Definition: colcompmatrix.hh:154
Utility class for converting an ISTL Matrix into a column-compressed matrix.
Definition: colcompmatrix.hh:145
Traits for type conversions and type information.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)