Loading...
Searching...
No Matches
bcrsmatrix.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_BCRSMATRIX_HH
#define DUNE_ISTL_BCRSMATRIX_HH
#include <cmath>
#include <complex>
#include <set>
#include <iostream>
#include <algorithm>
#include <numeric>
#include <vector>
#include <map>
#include <memory>
#include "istlexception.hh"
#include "bvector.hh"
#include "matrixutils.hh"
#include "matrixindexset.hh"
#include <dune/common/stdstreams.hh>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
#include <dune/istl/blocklevel.hh>
namespace Dune {
template<typename M>
struct MatrixDimension;
template<typename size_type>
struct CompressionStatistics
{
double avg;
size_type maximum;
size_type overflow_total;
double mem_ratio;
};
template<class M_>
class ImplicitMatrixBuilder
{
public:
typedef M_ Matrix;
typedef typename Matrix::block_type block_type;
typedef typename Matrix::size_type size_type;
class row_object
{
public:
block_type& operator[](size_type j) const
{
return _m.entry(_i,j);
}
#ifndef DOXYGEN
row_object(Matrix& m, size_type i)
: _m(m)
, _i(i)
{}
#endif
private:
Matrix& _m;
size_type _i;
};
ImplicitMatrixBuilder(Matrix& m)
: _m(m)
{
if (m.buildMode() != Matrix::implicit)
DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
if (m.buildStage() != Matrix::building)
DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
}
ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
: _m(m)
{
if (m.buildStage() != Matrix::notAllocated)
DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
m.setBuildMode(Matrix::implicit);
m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
m.setSize(rows,cols);
}
row_object operator[](size_type i) const
{
return row_object(_m,i);
}
size_type N() const
{
return _m.N();
}
size_type M() const
{
return _m.M();
}
private:
Matrix& _m;
};
template<class B, class A=std::allocator<B> >
class BCRSMatrix
{
friend struct MatrixDimension<BCRSMatrix>;
public:
enum BuildStage {
notbuilt=0,
notAllocated=0,
building=1,
rowSizesBuilt=2,
built=3
};
//===== type definitions and constants
using field_type = typename Imp::BlockTraits<B>::field_type;
typedef B block_type;
typedef A allocator_type;
typedef typename A::size_type size_type;
typedef Imp::CompressedBlockVectorWindow<B,size_type> row_type;
typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
enum BuildMode {
row_wise,
random,
implicit,
unknown
};
//===== random access interface to rows of the matrix
row_type& operator[] (size_type i)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (build_mode == implicit && ready != built)
DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
return r[i];
}
const row_type& operator[] (size_type i) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (build_mode == implicit && ready != built)
DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
return r[i];
}
//===== iterator interface to rows of the matrix
template<class T>
class RealRowIterator
: public RandomAccessIteratorFacade<RealRowIterator<T>, T>
{
public:
typedef typename std::remove_const<T>::type ValueType;
friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>;
friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>;
friend class RealRowIterator<const ValueType>;
friend class RealRowIterator<ValueType>;
RealRowIterator (row_type* _p, size_type _i)
: p(_p), i(_i)
{}
RealRowIterator ()
: p(0), i(0)
{}
RealRowIterator(const RealRowIterator<ValueType>& it)
: p(it.p), i(it.i)
{}
size_type index () const
{
return i;
}
std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
{
assert(other.p==p);
return (other.i-i);
}
std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
{
assert(other.p==p);
return (other.i-i);
}
bool equals (const RealRowIterator<ValueType>& other) const
{
assert(other.p==p);
return i==other.i;
}
bool equals (const RealRowIterator<const ValueType>& other) const
{
assert(other.p==p);
return i==other.i;
}
private:
void increment()
{
++i;
}
void decrement()
{
--i;
}
void advance(std::ptrdiff_t diff)
{
i+=diff;
}
T& elementAt(std::ptrdiff_t diff) const
{
return p[i+diff];
}
row_type& dereference () const
{
return p[i];
}
row_type* p;
size_type i;
};
typedef RealRowIterator<row_type> iterator;
typedef RealRowIterator<row_type> Iterator;
Iterator begin ()
{
return Iterator(r,0);
}
Iterator end ()
{
return Iterator(r,n);
}
Iterator beforeEnd ()
{
return Iterator(r,n-1);
}
Iterator beforeBegin ()
{
return Iterator(r,-1);
}
typedef Iterator RowIterator;
typedef typename row_type::Iterator ColIterator;
typedef RealRowIterator<const row_type> const_iterator;
typedef RealRowIterator<const row_type> ConstIterator;
ConstIterator begin () const
{
return ConstIterator(r,0);
}
ConstIterator end () const
{
return ConstIterator(r,n);
}
ConstIterator beforeEnd() const
{
return ConstIterator(r,n-1);
}
ConstIterator beforeBegin () const
{
return ConstIterator(r,-1);
}
typedef ConstIterator ConstRowIterator;
typedef typename row_type::ConstIterator ConstColIterator;
//===== constructors & resizers
// we use a negative compressionBufferSize to indicate that the implicit
// mode parameters have not been set yet
BCRSMatrix ()
: build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
allocationSize_(0), r(0), a(0),
avg(0), compressionBufferSize_(-1.0)
{}
BCRSMatrix (size_type _n, size_type _m, size_type _nnz, BuildMode bm)
: build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
allocationSize_(0), r(0), a(0),
avg(0), compressionBufferSize_(-1.0)
{
allocate(_n, _m, _nnz,true,false);
}
BCRSMatrix (size_type _n, size_type _m, BuildMode bm)
: build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
allocationSize_(0), r(0), a(0),
avg(0), compressionBufferSize_(-1.0)
{
allocate(_n, _m,0,true,false);
}
BCRSMatrix (size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
: build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
allocationSize_(0), r(0), a(0),
avg(_avg), compressionBufferSize_(compressionBufferSize)
{
if (bm != implicit)
DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
// Prevent user from setting a negative compression buffer size:
// 1) It doesn't make sense
// 2) We use a negative value to indicate that the parameters
// have not been set yet
if (compressionBufferSize_ < 0.0)
DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
implicit_allocate(_n,_m);
}
BCRSMatrix (const BCRSMatrix& Mat)
: build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
allocationSize_(0), r(0), a(0),
avg(Mat.avg), compressionBufferSize_(Mat.compressionBufferSize_)
{
if (!(Mat.ready == notAllocated || Mat.ready == built))
DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
// deep copy in global array
size_type _nnz = Mat.nonzeroes();
j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
allocate(Mat.n, Mat.m, _nnz, true, true);
// build window structure
copyWindowStructure(Mat);
}
~BCRSMatrix ()
{
deallocate();
}
void setBuildMode(BuildMode bm)
{
if (ready == notAllocated)
{
build_mode = bm;
return;
}
if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
build_mode = bm;
else
DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
}
void setSize(size_type rows, size_type columns, size_type nnz=0)
{
// deallocate already setup memory
deallocate();
if (build_mode == implicit)
{
if (nnz>0)
DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
// implicit allocates differently
implicit_allocate(rows,columns);
}
else
{
// allocate matrix memory
allocate(rows, columns, nnz, true, false);
}
}
void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
{
// Prevent user from setting a negative compression buffer size:
// 1) It doesn't make sense
// 2) We use a negative value to indicate that the parameters
// have not been set yet
if (compressionBufferSize < 0.0)
DUNE_THROW(BCRSMatrixError,"You cannot set a negative compressionBufferSize value");
// make sure the parameters aren't changed after memory has been allocated
if (ready != notAllocated)
DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
avg = _avg;
compressionBufferSize_ = compressionBufferSize;
}
BCRSMatrix& operator= (const BCRSMatrix& Mat)
{
// return immediately when self-assignment
if (&Mat==this) return *this;
if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
// make it simple: ALWAYS throw away memory for a and j_
// and deallocate rows only if n != Mat.n
deallocate(n!=Mat.n);
// reallocate the rows if required
if (n>0 && n!=Mat.n) {
// free rows
for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
rowAllocator_.deallocate(r,n);
}
nnz_ = Mat.nonzeroes();
// allocate a, share j_
j_ = Mat.j_;
allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
// build window structure
copyWindowStructure(Mat);
return *this;
}
BCRSMatrix& operator= (const field_type& k)
{
if (!(ready == notAllocated || ready == built))
DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
for (size_type i=0; i<n; i++) r[i] = k;
return *this;
}
//===== row-wise creation interface
class CreateIterator
{
public:
CreateIterator (BCRSMatrix& _Mat, size_type _i)
: Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
{
if (Mat.build_mode == unknown && Mat.ready == building)
{
Mat.build_mode = row_wise;
}
if (i==0 && Mat.ready != building)
DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
if(Mat.build_mode!=row_wise)
DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
if(i==0 && _Mat.N()==0)
// empty Matrix is always built.
Mat.ready = built;
}
CreateIterator& operator++()
{
// this should only be called if matrix is in creation
if (Mat.ready != building)
DUNE_THROW(BCRSMatrixError,"matrix already built up");
// row i is defined through the pattern
// get memory for the row and initialize the j_ array
// this depends on the allocation mode
// compute size of the row
size_type s = pattern.size();
if(s>0) {
// update number of nonzeroes including this row
nnz += s;
// alloc memory / set window
if (Mat.nnz_ > 0)
{
// memory is allocated in one long array
// check if that memory is sufficient
if (nnz > Mat.nnz_)
DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
// set row i
Mat.r[i].set(s,nullptr,current_row.getindexptr());
current_row.setindexptr(current_row.getindexptr()+s);
}else{
// memory is allocated individually per row
// allocate and set row i
B* b = Mat.allocator_.allocate(s);
// use placement new to call constructor that allocates
// additional memory.
new (b) B[s];
size_type* j = Mat.sizeAllocator_.allocate(s);
Mat.r[i].set(s,b,j);
}
}else
// setup empty row
Mat.r[i].set(0,nullptr,nullptr);
// initialize the j array for row i from pattern
std::visit([&](const auto& row_indices){
std::copy(row_indices.cbegin(), row_indices.cend(), Mat.r[i].getindexptr());
}, pattern.storage());
// now go to next row
i++;
pattern.clear();
// check if this was last row
if (i==Mat.n)
{
Mat.ready = built;
if(Mat.nnz_ > 0)
{
// Set nnz to the exact number of nonzero blocks inserted
// as some methods rely on it
Mat.nnz_ = nnz;
// allocate data array
Mat.allocateData();
Mat.setDataPointers();
}
}
// done
return *this;
}
bool operator!= (const CreateIterator& it) const
{
return (i!=it.i) || (&Mat!=&it.Mat);
}
bool operator== (const CreateIterator& it) const
{
return (i==it.i) && (&Mat==&it.Mat);
}
size_type index () const
{
return i;
}
void insert (size_type j)
{
pattern.insert(j);
}
bool contains (size_type j)
{
return pattern.contains(j);
}
size_type size() const
{
return pattern.size();
}
private:
BCRSMatrix& Mat; // the matrix we are defining
size_type i; // current row to be defined
size_type nnz; // count total number of nonzeros
using PatternType = typename Impl::RowIndexSet;
PatternType pattern; // used to compile entries in a row
row_type current_row; // row pointing to the current row to setup
};
friend class CreateIterator;
CreateIterator createbegin ()
{
return CreateIterator(*this,0);
}
CreateIterator createend ()
{
return CreateIterator(*this,n);
}
//===== random creation interface
void setrowsize (size_type i, size_type s)
{
if (build_mode!=random)
DUNE_THROW(BCRSMatrixError,"requires random build mode");
if (ready != building)
DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
r[i].setsize(s);
}
size_type getrowsize (size_type i) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
return r[i].getsize();
}
void incrementrowsize (size_type i, size_type s = 1)
{
if (build_mode!=random)
DUNE_THROW(BCRSMatrixError,"requires random build mode");
if (ready != building)
DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
r[i].setsize(r[i].getsize()+s);
}
void endrowsizes ()
{
if (build_mode!=random)
DUNE_THROW(BCRSMatrixError,"requires random build mode");
if (ready != building)
DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
// compute total size, check positivity
size_type total=0;
for (size_type i=0; i<n; i++)
{
total += r[i].getsize();
}
if(nnz_ == 0)
// allocate/check memory
allocate(n,m,total,false,false);
else if(nnz_ < total)
DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
<<"sufficient for calculated nonzeros ("<<total<<"! ");
// set the window pointers correctly
setColumnPointers(begin());
// initialize j_ array with m (an invalid column index)
// this indicates an unused entry
for (size_type k=0; k<nnz_; k++)
j_.get()[k] = m;
ready = rowSizesBuilt;
}
void addindex (size_type row, size_type col)
{
if (build_mode!=random)
DUNE_THROW(BCRSMatrixError,"requires random build mode");
if (ready==built)
DUNE_THROW(BCRSMatrixError,"matrix already built up");
if (ready==building)
DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
if (ready==notAllocated)
DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
if (col >= m)
DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
// get row range
size_type* const first = r[row].getindexptr();
size_type* const last = first + r[row].getsize();
// find correct insertion position for new column index
size_type* pos = std::lower_bound(first,last,col);
// check if index is already in row
if (pos!=last && *pos == col) return;
// find end of already inserted column indices
size_type* end = std::lower_bound(pos,last,m);
if (end==last)
DUNE_THROW(BCRSMatrixError,"row is too small");
// insert new column index at correct position
std::copy_backward(pos,end,end+1);
*pos = col;
}
template<typename It>
void setIndicesNoSort(size_type row, It begin, It end)
{
size_type row_size = r[row].size();
size_type* col_begin = r[row].getindexptr();
size_type* col_end;
// consistency check between allocated row size and number of passed column indices
if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
DUNE_THROW(BCRSMatrixError,"Given size of row " << row
<< " (" << row_size
<< ") does not match number of passed entries (" << (col_end - col_begin) << ")");
}
template<typename It>
void setIndices(size_type row, It begin, It end)
{
size_type row_size = r[row].size();
size_type* col_begin = r[row].getindexptr();
size_type* col_end;
// consistency check between allocated row size and number of passed column indices
if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
DUNE_THROW(BCRSMatrixError,"Given size of row " << row
<< " (" << row_size
<< ") does not match number of passed entries (" << (col_end - col_begin) << ")");
std::sort(col_begin,col_end);
}
void endindices ()
{
if (build_mode!=random)
DUNE_THROW(BCRSMatrixError,"requires random build mode");
if (ready==built)
DUNE_THROW(BCRSMatrixError,"matrix already built up");
if (ready==building)
DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
if (ready==notAllocated)
DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
// check if there are undefined indices
RowIterator endi=end();
for (RowIterator i=begin(); i!=endi; ++i)
{
ColIterator endj = (*i).end();
for (ColIterator j=(*i).begin(); j!=endj; ++j) {
if (j.index() >= m) {
dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
<<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
nnz_ -= ((*i).end().offset() - j.offset());
r[i.index()].setsize(j.offset());
break;
}
}
}
allocateData();
setDataPointers();
// if not, set matrix to built
ready = built;
}
//===== implicit creation interface
B& entry(size_type row, size_type col)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (build_mode!=implicit)
DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
if (ready==built)
DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
if (ready==notAllocated)
DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
if (ready!=building)
DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
if (row >= n)
DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
if (col >= m)
DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
#endif
size_type* begin = r[row].getindexptr();
size_type* end = begin + r[row].getsize();
size_type* pos = std::find(begin, end, col);
//treat the case that there was a match in the array
if (pos != end)
if (*pos == col)
{
std::ptrdiff_t offset = pos - r[row].getindexptr();
B* aptr = r[row].getptr() + offset;
return *aptr;
}
//determine whether overflow has to be taken into account or not
if (r[row].getsize() == avg)
return overflow[std::make_pair(row,col)];
else
{
//modify index array
*end = col;
//do simultaneous operations on data array a
std::ptrdiff_t offset = end - r[row].getindexptr();
B* apos = r[row].getptr() + offset;
//increase rowsize
r[row].setsize(r[row].getsize()+1);
//return reference to the newly created entry
return *apos;
}
}
CompressionStatistics compress()
{
if (build_mode!=implicit)
DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
if (ready==built)
DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
if (ready==notAllocated)
DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
if (ready!=building)
DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
//calculate statistics
CompressionStatistics stats;
stats.overflow_total = overflow.size();
stats.maximum = 0;
//get insertion iterators pointing to one before start (for later use of ++it)
size_type* jiit = j_.get();
B* aiit = a;
//get iterator to the smallest overflow element
typename OverflowType::iterator oit = overflow.begin();
//store a copy of index pointers on which to perform sorting
std::vector<size_type*> perm;
//iterate over all rows and copy elements into their position in the compressed array
for (size_type i=0; i<n; i++)
{
//get old pointers into a and j and size without overflow changes
size_type* begin = r[i].getindexptr();
//B* apos = r[i].getptr();
size_type size = r[i].getsize();
perm.resize(size);
typename std::vector<size_type*>::iterator it = perm.begin();
for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
*it = iit;
//sort permutation array
std::sort(perm.begin(),perm.end(), [](const size_type* l, const size_type* r){ return *l < *r; });
//change row window pointer to their new positions
r[i].setindexptr(jiit);
r[i].setptr(aiit);
for (it = perm.begin(); it != perm.end(); ++it)
{
//check whether there are elements in the overflow area which take precedence
while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
{
//check whether there is enough memory to write to
if (jiit > begin)
DUNE_THROW(Dune::ImplicitModeCompressionBufferExhausted,
"Allocated memory for BCRSMatrix exhausted during compress()!"
"Please increase either the average number of entries per row or the compressionBufferSize value."
);
//copy an element from the overflow area to the insertion position in a and j
*jiit = oit->first.second;
++jiit;
*aiit = oit->second;
++aiit;
++oit;
r[i].setsize(r[i].getsize()+1);
}
//check whether there is enough memory to write to
if (jiit > begin)
DUNE_THROW(Dune::ImplicitModeCompressionBufferExhausted,
"Allocated memory for BCRSMatrix exhausted during compress()!"
"Please increase either the average number of entries per row or the compressionBufferSize value."
);
//copy element from array
*jiit = **it;
++jiit;
B* apos = *it - j_.get() + a;
*aiit = *apos;
++aiit;
}
//copy remaining elements from the overflow area
while ((oit!=overflow.end()) && (oit->first.first == i))
{
//check whether there is enough memory to write to
if (jiit > begin)
DUNE_THROW(Dune::ImplicitModeCompressionBufferExhausted,
"Allocated memory for BCRSMatrix exhausted during compress()!"
"Please increase either the average number of entries per row or the compressionBufferSize value."
);
//copy and element from the overflow area to the insertion position in a and j
*jiit = oit->first.second;
++jiit;
*aiit = oit->second;
++aiit;
++oit;
r[i].setsize(r[i].getsize()+1);
}
// update maximum row size
if (r[i].getsize()>stats.maximum)
stats.maximum = r[i].getsize();
}
// overflow area may be cleared
overflow.clear();
//determine average number of entries and memory usage
if ( n == 0)
{
stats.avg = 0;
stats.mem_ratio = 1;
}
else
{
std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
nnz_ = diff;
stats.avg = (double) (nnz_) / (double) n;
stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
}
//matrix is now built
ready = built;
return stats;
}
//===== vector space arithmetic
BCRSMatrix& operator*= (const field_type& k)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif
if (nnz_ > 0)
{
// process 1D array
for (size_type i=0; i<nnz_; i++)
a[i] *= k;
}
else
{
RowIterator endi=end();
for (RowIterator i=begin(); i!=endi; ++i)
{
ColIterator endj = (*i).end();
for (ColIterator j=(*i).begin(); j!=endj; ++j)
(*j) *= k;
}
}
return *this;
}
BCRSMatrix& operator/= (const field_type& k)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif
if (nnz_ > 0)
{
// process 1D array
for (size_type i=0; i<nnz_; i++)
a[i] /= k;
}
else
{
RowIterator endi=end();
for (RowIterator i=begin(); i!=endi; ++i)
{
ColIterator endj = (*i).end();
for (ColIterator j=(*i).begin(); j!=endj; ++j)
(*j) /= k;
}
}
return *this;
}
BCRSMatrix& operator+= (const BCRSMatrix& b)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built || b.ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if(N()!=b.N() || M() != b.M())
DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
RowIterator endi=end();
ConstRowIterator j=b.begin();
for (RowIterator i=begin(); i!=endi; ++i, ++j) {
i->operator+=(*j);
}
return *this;
}
BCRSMatrix& operator-= (const BCRSMatrix& b)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built || b.ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if(N()!=b.N() || M() != b.M())
DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
RowIterator endi=end();
ConstRowIterator j=b.begin();
for (RowIterator i=begin(); i!=endi; ++i, ++j) {
i->operator-=(*j);
}
return *this;
}
BCRSMatrix& axpy(field_type alpha, const BCRSMatrix& b)
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built || b.ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if(N()!=b.N() || M() != b.M())
DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
RowIterator endi=end();
ConstRowIterator j=b.begin();
for(RowIterator i=begin(); i!=endi; ++i, ++j)
i->axpy(alpha, *j);
return *this;
}
//===== linear maps
template<class X, class Y>
void mv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
"Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
"Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
y[i.index()]=0;
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xj = Impl::asVector(x[j.index()]);
auto&& yi = Impl::asVector(y[i.index()]);
Impl::asMatrix(*j).umv(xj, yi);
}
}
}
template<class X, class Y>
void umv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xj = Impl::asVector(x[j.index()]);
auto&& yi = Impl::asVector(y[i.index()]);
Impl::asMatrix(*j).umv(xj,yi);
}
}
}
template<class X, class Y>
void mmv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xj = Impl::asVector(x[j.index()]);
auto&& yi = Impl::asVector(y[i.index()]);
Impl::asMatrix(*j).mmv(xj,yi);
}
}
}
template<class X, class Y, class F>
void usmv (F&& alpha, const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xj = Impl::asVector(x[j.index()]);
auto&& yi = Impl::asVector(y[i.index()]);
Impl::asMatrix(*j).usmv(alpha,xj,yi);
}
}
}
template<class X, class Y>
void mtv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
for(size_type i=0; i<y.N(); ++i)
y[i]=0;
umtv(x,y);
}
template<class X, class Y>
void umtv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xi = Impl::asVector(x[i.index()]);
auto&& yj = Impl::asVector(y[j.index()]);
Impl::asMatrix(*j).umtv(xi,yj);
}
}
}
template<class X, class Y>
void mmtv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xi = Impl::asVector(x[i.index()]);
auto&& yj = Impl::asVector(y[j.index()]);
Impl::asMatrix(*j).mmtv(xi,yj);
}
}
}
template<class X, class Y>
void usmtv (const field_type& alpha, const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xi = Impl::asVector(x[i.index()]);
auto&& yj = Impl::asVector(y[j.index()]);
Impl::asMatrix(*j).usmtv(alpha,xi,yj);
}
}
}
template<class X, class Y>
void umhv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xi = Impl::asVector(x[i.index()]);
auto&& yj = Impl::asVector(y[j.index()]);
Impl::asMatrix(*j).umhv(xi,yj);
}
}
}
template<class X, class Y>
void mmhv (const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xi = Impl::asVector(x[i.index()]);
auto&& yj = Impl::asVector(y[j.index()]);
Impl::asMatrix(*j).mmhv(xi,yj);
}
}
}
template<class X, class Y>
void usmhv (const field_type& alpha, const X& x, Y& y) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
ConstRowIterator endi=end();
for (ConstRowIterator i=this->begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
{
auto&& xi = Impl::asVector(x[i.index()]);
auto&& yj = Impl::asVector(y[j.index()]);
Impl::asMatrix(*j).usmhv(alpha,xi,yj);
}
}
}
//===== norms
typename FieldTraits<field_type>::real_type frobenius_norm2 () const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif
typename FieldTraits<field_type>::real_type sum=0;
for (auto&& row : (*this))
for (auto&& entry : row)
sum += Impl::asMatrix(entry).frobenius_norm2();
return sum;
}
typename FieldTraits<field_type>::real_type frobenius_norm () const
{
return sqrt(frobenius_norm2());
}
template <typename ft = field_type,
typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
typename FieldTraits<ft>::real_type infinity_norm() const {
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
using real_type = typename FieldTraits<ft>::real_type;
using std::max;
real_type norm = 0;
for (auto const &x : *this) {
real_type sum = 0;
for (auto const &y : x)
sum += Impl::asMatrix(y).infinity_norm();
norm = max(sum, norm);
}
return norm;
}
template <typename ft = field_type,
typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
typename FieldTraits<ft>::real_type infinity_norm_real() const {
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
using real_type = typename FieldTraits<ft>::real_type;
using std::max;
real_type norm = 0;
for (auto const &x : *this) {
real_type sum = 0;
for (auto const &y : x)
sum += Impl::asMatrix(y).infinity_norm_real();
norm = max(sum, norm);
}
return norm;
}
template <typename ft = field_type,
typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
typename FieldTraits<ft>::real_type infinity_norm() const {
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
using real_type = typename FieldTraits<ft>::real_type;
using std::max;
real_type norm = 0;
real_type is_nan = 1;
for (auto const &x : *this) {
real_type sum = 0;
for (auto const &y : x)
sum += Impl::asMatrix(y).infinity_norm();
norm = max(sum, norm);
is_nan += sum;
}
return norm * (is_nan / is_nan);
}
template <typename ft = field_type,
typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
typename FieldTraits<ft>::real_type infinity_norm_real() const {
if (ready != built)
DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
using real_type = typename FieldTraits<ft>::real_type;
using std::max;
real_type norm = 0;
real_type is_nan = 1;
for (auto const &x : *this) {
real_type sum = 0;
for (auto const &y : x)
sum += Impl::asMatrix(y).infinity_norm_real();
norm = max(sum, norm);
is_nan += sum;
}
return norm * (is_nan / is_nan);
}
//===== sizes
size_type N () const
{
return n;
}
size_type M () const
{
return m;
}
size_type nonzeroes () const
{
// in case of row-wise allocation
if( nnz_ <= 0 )
nnz_ = std::accumulate( begin(), end(), size_type( 0 ), [] ( size_type s, const row_type &row ) { return s+row.getsize(); } );
return nnz_;
}
BuildStage buildStage() const
{
return ready;
}
BuildMode buildMode() const
{
return build_mode;
}
//===== query
bool exists (size_type i, size_type j) const
{
#ifdef DUNE_ISTL_WITH_CHECKING
if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
#endif
return (r[i].size() && r[i].find(j) != r[i].end());
}
protected:
// state information
BuildMode build_mode; // row wise or whole matrix
BuildStage ready; // indicate the stage the matrix building is in
// The allocator used for memory management
typename std::allocator_traits<A>::template rebind_alloc<B> allocator_;
typename std::allocator_traits<A>::template rebind_alloc<row_type> rowAllocator_;
typename std::allocator_traits<A>::template rebind_alloc<size_type> sizeAllocator_;
// size of the matrix
size_type n; // number of rows
size_type m; // number of columns
mutable size_type nnz_; // number of nonzeroes contained in the matrix
size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
// zero means that memory is allocated separately for each row.
// the rows are dynamically allocated
row_type* r; // [n] the individual rows having pointers into a,j arrays
// dynamically allocated memory
B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
// If a single array of column indices is used, it can be shared
// between different matrices with the same sparsity pattern
std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
// additional data is needed in implicit buildmode
size_type avg;
double compressionBufferSize_;
typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
OverflowType overflow;
void setWindowPointers(ConstRowIterator row)
{
row_type current_row(a,j_.get(),0); // Pointers to current row data
for (size_type i=0; i<n; i++, ++row) {
// set row i
size_type s = row->getsize();
if (s>0) {
// setup pointers and size
r[i].set(s,current_row.getptr(), current_row.getindexptr());
// update pointer for next row
current_row.setptr(current_row.getptr()+s);
current_row.setindexptr(current_row.getindexptr()+s);
} else{
// empty row
r[i].set(0,nullptr,nullptr);
}
}
}
void setColumnPointers(ConstRowIterator row)
{
size_type* jptr = j_.get();
for (size_type i=0; i<n; ++i, ++row) {
// set row i
size_type s = row->getsize();
if (s>0) {
// setup pointers and size
r[i].setsize(s);
r[i].setindexptr(jptr);
} else{
// empty row
r[i].set(0,nullptr,nullptr);
}
// advance position in global array
jptr += s;
}
}
void setDataPointers()
{
B* aptr = a;
for (size_type i=0; i<n; ++i) {
// set row i
if (r[i].getsize() > 0) {
// setup pointers and size
r[i].setptr(aptr);
} else{
// empty row
r[i].set(0,nullptr,nullptr);
}
// advance position in global array
aptr += r[i].getsize();
}
}
void copyWindowStructure(const BCRSMatrix& Mat)
{
setWindowPointers(Mat.begin());
// copy data
for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
// finish off
build_mode = row_wise; // dummy
ready = built;
}
void deallocate(bool deallocateRows=true)
{
if (notAllocated)
return;
if (allocationSize_>0)
{
// a,j_ have been allocated as one long vector
j_.reset();
if (a)
{
for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
std::allocator_traits<decltype(allocator_)>::destroy(allocator_, aiter);
allocator_.deallocate(a,allocationSize_);
a = nullptr;
}
}
else if (r)
{
// check if memory for rows have been allocated individually
for (size_type i=0; i<n; i++)
if (r[i].getsize()>0)
{
for (B *col=r[i].getptr()+(r[i].getsize()-1),
*colend = r[i].getptr()-1; col!=colend; --col) {
std::allocator_traits<decltype(allocator_)>::destroy(allocator_, col);
}
sizeAllocator_.deallocate(r[i].getindexptr(),1);
allocator_.deallocate(r[i].getptr(),1);
// clear out row data in case we don't want to deallocate the rows
// otherwise we might run into a double free problem here later
r[i].set(0,nullptr,nullptr);
}
}
// deallocate the rows
if (n>0 && deallocateRows && r) {
for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
rowAllocator_.deallocate(r,n);
r = nullptr;
}
// Mark matrix as not built at all.
ready=notAllocated;
}
void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
{
// Store size
n = rows;
m = columns;
nnz_ = allocationSize;
allocationSize_ = allocationSize;
// allocate rows
if(allocateRows) {
if (n>0) {
if (r)
DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
r = rowAllocator_.allocate(rows);
// initialize row entries
for(row_type* ri=r; ri!=r+rows; ++ri)
std::allocator_traits<decltype(rowAllocator_)>::construct(rowAllocator_, ri, row_type());
}else{
r = 0;
}
}
// allocate a and j_ array
if (allocate_data)
allocateData();
// allocate column indices only if not yet present (enable sharing)
if (allocationSize_>0) {
// we copy allocator and size to the deleter since _j may outlive this class
if (!j_.get())
j_.reset(sizeAllocator_.allocate(allocationSize_),
[alloc = sizeAllocator_, size = allocationSize_](auto ptr) mutable {
alloc.deallocate(ptr, size);
});
}else{
j_.reset();
}
// Mark the matrix as not built.
ready = building;
}
void allocateData()
{
if (a)
DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
if (allocationSize_>0) {
a = allocator_.allocate(allocationSize_);
// use placement new to call constructor that allocates
// additional memory.
new (a) B[allocationSize_];
} else {
a = nullptr;
}
}
void implicit_allocate(size_type _n, size_type _m)
{
if (build_mode != implicit)
DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
if (ready != notAllocated)
DUNE_THROW(InvalidStateException,"memory has already been allocated");
// check to make sure the user has actually set the parameters
if (compressionBufferSize_ < 0)
DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
//calculate size of overflow region, add buffer for row sort!
size_type osize = (size_type) (_n*avg)*compressionBufferSize_ + 4*avg;
allocationSize_ = _n*avg + osize;
allocate(_n, _m, allocationSize_,true,true);
//set row pointers correctly
size_type* jptr = j_.get() + osize;
B* aptr = a + osize;
for (size_type i=0; i<n; i++)
{
r[i].set(0,aptr,jptr);
jptr = jptr + avg;
aptr = aptr + avg;
}
ready = building;
}
};
template<class B, class A>
struct FieldTraits< BCRSMatrix<B, A> >
{
using field_type = typename BCRSMatrix<B, A>::field_type;
using real_type = typename FieldTraits<field_type>::real_type;
};
} // end namespace
#endif
matrixutils.hh
Some handy generic functions for ISTL matrices.
matrixindexset.hh
blocklevel.hh
Helper functions for determining the vector/matrix block level.
istlexception.hh
bvector.hh
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Dune
Definition allocator.hh:11
Dune::get
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:293
Dune::col
MatrixSparsityPatternGatherScatter< M, I >::ColIter MatrixSparsityPatternGatherScatter< M, I >::col
Definition matrixredistribute.hh:591
Dune::type
const Dtype_t BaseGetSuperLUType< T >::type
Definition supermatrix.hh:135
Dune::MatrixDimension
Definition matrixutils.hh:211
Dune::CompressionStatistics
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:89
Dune::CompressionStatistics::overflow_total
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition bcrsmatrix.hh:95
Dune::CompressionStatistics::maximum
size_type maximum
maximum number of non-zeroes per row.
Definition bcrsmatrix.hh:93
Dune::CompressionStatistics::avg
double avg
average number of non-zeroes per row.
Definition bcrsmatrix.hh:91
Dune::CompressionStatistics::mem_ratio
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition bcrsmatrix.hh:100
Dune::ImplicitMatrixBuilder::block_type
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition bcrsmatrix.hh:126
Dune::ImplicitMatrixBuilder::ImplicitMatrixBuilder
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition bcrsmatrix.hh:171
Dune::ImplicitMatrixBuilder::Matrix
M_ Matrix
The underlying matrix.
Definition bcrsmatrix.hh:123
Dune::ImplicitMatrixBuilder::ImplicitMatrixBuilder
ImplicitMatrixBuilder(Matrix &m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
Sets up matrix m for implicit construction using the given parameters and creates an ImplicitBmatrixu...
Definition bcrsmatrix.hh:195
Dune::ImplicitMatrixBuilder::M
size_type M() const
The number of columns in the matrix.
Definition bcrsmatrix.hh:218
Dune::ImplicitMatrixBuilder::size_type
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition bcrsmatrix.hh:129
Dune::ImplicitMatrixBuilder::operator[]
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition bcrsmatrix.hh:206
Dune::ImplicitMatrixBuilder::N
size_type N() const
The number of rows in the matrix.
Definition bcrsmatrix.hh:212
Dune::ImplicitMatrixBuilder::row_object
Proxy row object for entry access.
Definition bcrsmatrix.hh:138
Dune::ImplicitMatrixBuilder::row_object::operator[]
block_type & operator[](size_type j) const
Returns entry in column j.
Definition bcrsmatrix.hh:143
Dune::size_type
Dune::BCRSMatrix::field_type
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition bcrsmatrix.hh:489
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::rowAllocator_
std::allocator_traits< A >::template rebind_alloc< row_type > rowAllocator_
Definition bcrsmatrix.hh:2057
Dune::BCRSMatrix::exists
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition bcrsmatrix.hh:2039
Dune::BCRSMatrix::buildStage
BuildStage buildStage() const
The current build stage of the matrix.
Definition bcrsmatrix.hh:2025
Dune::BCRSMatrix::begin
Iterator begin()
Get iterator to first row.
Definition bcrsmatrix.hh:672
Dune::BCRSMatrix::CreateIterator
friend class CreateIterator
allow CreateIterator to access internal data
Definition bcrsmatrix.hh:1094
Dune::BCRSMatrix::copyWindowStructure
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition bcrsmatrix.hh:2154
Dune::BCRSMatrix::entry
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition bcrsmatrix.hh:1320
Dune::BCRSMatrix::mmhv
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition bcrsmatrix.hh:1844
Dune::BCRSMatrix::usmhv
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition bcrsmatrix.hh:1867
Dune::BCRSMatrix::umtv
void umtv(const X &x, Y &y) const
y += A^T x
Definition bcrsmatrix.hh:1754
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::compressionBufferSize_
double compressionBufferSize_
Definition bcrsmatrix.hh:2079
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::m
size_type m
Definition bcrsmatrix.hh:2063
Dune::BCRSMatrix::const_iterator
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition bcrsmatrix.hh:704
Dune::BCRSMatrix::allocate
void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
Allocate memory for the matrix structure.
Definition bcrsmatrix.hh:2238
Dune::BCRSMatrix::axpy
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition bcrsmatrix.hh:1624
Dune::BCRSMatrix::infinity_norm_real
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition bcrsmatrix.hh:1937
Dune::BCRSMatrix::~BCRSMatrix
~BCRSMatrix()
destructor
Definition bcrsmatrix.hh:822
Dune::BCRSMatrix::allocateData
void allocateData()
Definition bcrsmatrix.hh:2279
Dune::BCRSMatrix::deallocate
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition bcrsmatrix.hh:2171
Dune::BCRSMatrix::RowIterator
Iterator RowIterator
rename the iterators for easier access
Definition bcrsmatrix.hh:698
Dune::BCRSMatrix::operator[]
row_type & operator[](size_type i)
random access to the rows
Definition bcrsmatrix.hh:546
Dune::BCRSMatrix::BCRSMatrix
BCRSMatrix()
an empty matrix
Definition bcrsmatrix.hh:746
Dune::BCRSMatrix::endrowsizes
void endrowsizes()
indicate that size of all rows is defined
Definition bcrsmatrix.hh:1149
Dune::BCRSMatrix::incrementrowsize
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition bcrsmatrix.hh:1138
Dune::BCRSMatrix::mtv
void mtv(const X &x, Y &y) const
y = A^T x
Definition bcrsmatrix.hh:1739
Dune::BCRSMatrix::umhv
void umhv(const X &x, Y &y) const
y += A^H x
Definition bcrsmatrix.hh:1821
Dune::BCRSMatrix::nonzeroes
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero)
Definition bcrsmatrix.hh:2016
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::allocationSize_
size_type allocationSize_
Definition bcrsmatrix.hh:2065
Dune::BCRSMatrix::ConstRowIterator
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition bcrsmatrix.hh:735
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::ready
BuildStage ready
Definition bcrsmatrix.hh:2052
Dune::BCRSMatrix::setIndicesNoSort
void setIndicesNoSort(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition bcrsmatrix.hh:1235
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::build_mode
BuildMode build_mode
Definition bcrsmatrix.hh:2051
Dune::BCRSMatrix::setrowsize
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition bcrsmatrix.hh:1117
Dune::BCRSMatrix::Iterator
RealRowIterator< row_type > Iterator
Definition bcrsmatrix.hh:669
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::nnz_
size_type nnz_
Definition bcrsmatrix.hh:2064
Dune::BCRSMatrix::operator*=
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition bcrsmatrix.hh:1516
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::sizeAllocator_
std::allocator_traits< A >::template rebind_alloc< size_type > sizeAllocator_
Definition bcrsmatrix.hh:2059
Dune::BCRSMatrix::iterator
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition bcrsmatrix.hh:668
Dune::BCRSMatrix::usmtv
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition bcrsmatrix.hh:1798
Dune::BCRSMatrix::beforeBegin
ConstIterator beforeBegin() const
Definition bcrsmatrix.hh:729
Dune::BCRSMatrix::ConstIterator
RealRowIterator< const row_type > ConstIterator
Definition bcrsmatrix.hh:705
Dune::BCRSMatrix::beforeBegin
Iterator beforeBegin()
Definition bcrsmatrix.hh:692
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::a
FieldMatrix< T, n, m > * a
Definition bcrsmatrix.hh:2072
Dune::BCRSMatrix::BuildMode
BuildMode
we support two modes
Definition bcrsmatrix.hh:507
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::implicit
@ implicit
Definition bcrsmatrix.hh:536
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::unknown
@ unknown
Definition bcrsmatrix.hh:540
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::random
@ random
Definition bcrsmatrix.hh:527
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::row_wise
@ row_wise
Definition bcrsmatrix.hh:518
Dune::BCRSMatrix::BCRSMatrix
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition bcrsmatrix.hh:753
Dune::BCRSMatrix::row_type
Imp::CompressedBlockVectorWindow< B, size_type > row_type
implement row_type with compressed vector
Definition bcrsmatrix.hh:501
Dune::BCRSMatrix::CompressionStatistics
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress().
Definition bcrsmatrix.hh:504
Dune::BCRSMatrix::operator-=
BCRSMatrix & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition bcrsmatrix.hh:1599
Dune::BCRSMatrix::BCRSMatrix
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition bcrsmatrix.hh:803
Dune::BCRSMatrix::end
Iterator end()
Get iterator to one beyond last row.
Definition bcrsmatrix.hh:678
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::r
row_type * r
Definition bcrsmatrix.hh:2069
Dune::BCRSMatrix::setIndices
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition bcrsmatrix.hh:1258
Dune::BCRSMatrix::addindex
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition bcrsmatrix.hh:1191
Dune::BCRSMatrix::OverflowType
std::map< std::pair< size_type, size_type >, B > OverflowType
Definition bcrsmatrix.hh:2081
Dune::BCRSMatrix::ColIterator
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition bcrsmatrix.hh:701
Dune::BCRSMatrix::frobenius_norm
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition bcrsmatrix.hh:1909
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::size_type
A::size_type size_type
Definition bcrsmatrix.hh:498
Dune::BCRSMatrix::operator/=
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition bcrsmatrix.hh:1544
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::overflow
OverflowType overflow
Definition bcrsmatrix.hh:2082
Dune::BCRSMatrix::operator+=
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition bcrsmatrix.hh:1577
Dune::BCRSMatrix::BCRSMatrix
BCRSMatrix(size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
construct matrix with a known average number of entries per row
Definition bcrsmatrix.hh:782
Dune::BCRSMatrix::createend
CreateIterator createend()
get create iterator pointing to one after the last block
Definition bcrsmatrix.hh:1103
Dune::BCRSMatrix::frobenius_norm2
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition bcrsmatrix.hh:1892
Dune::BCRSMatrix::beforeEnd
Iterator beforeEnd()
Definition bcrsmatrix.hh:685
Dune::BCRSMatrix::ConstColIterator
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition bcrsmatrix.hh:738
Dune::BCRSMatrix::usmv
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition bcrsmatrix.hh:1716
Dune::BCRSMatrix::getrowsize
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition bcrsmatrix.hh:1128
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::M
size_type M() const
Definition bcrsmatrix.hh:2010
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::n
size_type n
Definition bcrsmatrix.hh:2062
Dune::BCRSMatrix::createbegin
CreateIterator createbegin()
get initial create iterator
Definition bcrsmatrix.hh:1097
Dune::BCRSMatrix::BuildStage
BuildStage
Definition bcrsmatrix.hh:470
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::rowSizesBuilt
@ rowSizesBuilt
Definition bcrsmatrix.hh:481
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::built
@ built
Definition bcrsmatrix.hh:483
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::notbuilt
@ notbuilt
Definition bcrsmatrix.hh:472
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::notAllocated
@ notAllocated
Definition bcrsmatrix.hh:474
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::building
@ building
Definition bcrsmatrix.hh:476
Dune::BCRSMatrix::buildMode
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition bcrsmatrix.hh:2031
Dune::BCRSMatrix::mmv
void mmv(const X &x, Y &y) const
y -= A x
Definition bcrsmatrix.hh:1693
Dune::BCRSMatrix::infinity_norm
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition bcrsmatrix.hh:1917
Dune::BCRSMatrix::mv
void mv(const X &x, Y &y) const
y = A x
Definition bcrsmatrix.hh:1644
Dune::BCRSMatrix::block_type
B block_type
export the type representing the components
Definition bcrsmatrix.hh:492
Dune::BCRSMatrix::mmtv
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition bcrsmatrix.hh:1777
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::avg
size_type avg
Definition bcrsmatrix.hh:2078
Dune::BCRSMatrix::umv
void umv(const X &x, Y &y) const
y += A x
Definition bcrsmatrix.hh:1670
Dune::BCRSMatrix::implicit_allocate
void implicit_allocate(size_type _n, size_type _m)
organizes allocation implicit mode calculates correct array size to be allocated and sets the the win...
Definition bcrsmatrix.hh:2298
Dune::BCRSMatrix::setImplicitBuildModeParameters
void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
Set parameters needed for creation in implicit build mode.
Definition bcrsmatrix.hh:887
Dune::BCRSMatrix::BCRSMatrix
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition bcrsmatrix.hh:762
Dune::BCRSMatrix::endindices
void endindices()
indicate that all indices are defined, check consistency
Definition bcrsmatrix.hh:1272
Dune::BCRSMatrix::compress
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition bcrsmatrix.hh:1384
Dune::BCRSMatrix::setDataPointers
void setDataPointers()
Set data pointers for all rows.
Definition bcrsmatrix.hh:2135
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::allocator_
std::allocator_traits< A >::template rebind_alloc< FieldMatrix< T, n, m > > allocator_
Definition bcrsmatrix.hh:2055
Dune::BCRSMatrix::N
size_type N() const
number of rows (counted in blocks)
Definition bcrsmatrix.hh:2004
Dune::BCRSMatrix::setBuildMode
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition bcrsmatrix.hh:831
Dune::BCRSMatrix::setSize
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition bcrsmatrix.hh:859
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A >::j_
std::shared_ptr< size_type > j_
Definition bcrsmatrix.hh:2075
Dune::BCRSMatrix::setWindowPointers
void setWindowPointers(ConstRowIterator row)
Definition bcrsmatrix.hh:2084
Dune::BCRSMatrix::operator=
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition bcrsmatrix.hh:909
Dune::BCRSMatrix::setColumnPointers
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition bcrsmatrix.hh:2109
Dune::BCRSMatrix::end
ConstIterator end() const
Get const iterator to one beyond last row.
Definition bcrsmatrix.hh:715
Dune::BCRSMatrix::begin
ConstIterator begin() const
Get const iterator to first row.
Definition bcrsmatrix.hh:709
Dune::BCRSMatrix::allocator_type
A allocator_type
export the allocator type
Definition bcrsmatrix.hh:495
Dune::BCRSMatrix::beforeEnd
ConstIterator beforeEnd() const
Definition bcrsmatrix.hh:722
Dune::BCRSMatrix::RealRowIterator
Iterator access to matrix rows
Definition bcrsmatrix.hh:576
Dune::BCRSMatrix::RealRowIterator::RealRowIterator
RealRowIterator()
empty constructor, use with care!
Definition bcrsmatrix.hh:593
Dune::BCRSMatrix::RealRowIterator::equals
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition bcrsmatrix.hh:621
Dune::BCRSMatrix::RealRowIterator< row_type >::ValueType
std::remove_const< row_type >::type ValueType
Definition bcrsmatrix.hh:580
Dune::BCRSMatrix::RealRowIterator::RealRowIterator
RealRowIterator(const RealRowIterator< ValueType > &it)
Definition bcrsmatrix.hh:597
Dune::BCRSMatrix::RealRowIterator::equals
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition bcrsmatrix.hh:628
Dune::BCRSMatrix::RealRowIterator::RealRowIterator
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition bcrsmatrix.hh:588
Dune::BCRSMatrix::RealRowIterator::distanceTo
std::ptrdiff_t distanceTo(const RealRowIterator< const ValueType > &other) const
Definition bcrsmatrix.hh:614
Dune::BCRSMatrix::RealRowIterator::index
size_type index() const
return index
Definition bcrsmatrix.hh:603
Dune::BCRSMatrix::RealRowIterator::distanceTo
std::ptrdiff_t distanceTo(const RealRowIterator< ValueType > &other) const
Definition bcrsmatrix.hh:608
Dune::BCRSMatrix::CreateIterator::operator==
bool operator==(const CreateIterator &it) const
equality
Definition bcrsmatrix.hh:1052
Dune::BCRSMatrix::CreateIterator::operator++
CreateIterator & operator++()
prefix increment
Definition bcrsmatrix.hh:975
Dune::BCRSMatrix::CreateIterator::index
size_type index() const
The number of the row that the iterator currently points to.
Definition bcrsmatrix.hh:1058
Dune::BCRSMatrix::CreateIterator::operator!=
bool operator!=(const CreateIterator &it) const
inequality
Definition bcrsmatrix.hh:1046
Dune::BCRSMatrix::CreateIterator::CreateIterator
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition bcrsmatrix.hh:958
Dune::BCRSMatrix::CreateIterator::insert
void insert(size_type j)
put column index in row
Definition bcrsmatrix.hh:1064
Dune::BCRSMatrix::CreateIterator::contains
bool contains(size_type j)
return true if column index is in row
Definition bcrsmatrix.hh:1070
Dune::BCRSMatrix::CreateIterator::size
size_type size() const
Get the current row size.
Definition bcrsmatrix.hh:1079
Dune::FieldTraits< BCRSMatrix< B, A > >::field_type
typename BCRSMatrix< B, A >::field_type field_type
Definition bcrsmatrix.hh:2332
Dune::FieldTraits< BCRSMatrix< B, A > >::real_type
typename FieldTraits< field_type >::real_type real_type
Definition bcrsmatrix.hh:2333
Dune::BCRSMatrixError
Error specific to BCRSMatrix.
Definition istlexception.hh:24
Dune::ImplicitModeCompressionBufferExhausted
Thrown when the compression buffer used by the implicit BCRSMatrix construction is exhausted.
Definition istlexception.hh:37
Dune::Matrix
A generic dynamic dense matrix.
Definition matrix.hh:561
Dune::Matrix::size_type
A::size_type size_type
Type for indices and sizes.
Definition matrix.hh:577
Dune::Matrix::block_type
T block_type
Export the type representing the components.
Definition matrix.hh:568