diff --git a/dune/istl/bvector.hh b/dune/istl/bvector.hh
index 02e88a1f4547f6c6f72b12eb038d57d84c8c4fdc..32b7e823b7464cbacdc661c0484ab6a7f4479eed 100644
--- a/dune/istl/bvector.hh
+++ b/dune/istl/bvector.hh
@@ -126,7 +126,7 @@ namespace Dune {
* @param y other (compatible) vector
* @return
*/
- template<class OtherB, class OtherA=std::allocator<OtherB> >
+ template<class OtherB, class OtherA>
typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
{
typedef typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType PromotedType;
@@ -147,7 +147,7 @@ namespace Dune {
* @param y other (compatible) vector
* @return
*/
- template<class OtherB, class OtherA=std::allocator<OtherB> >
+ template<class OtherB, class OtherA>
typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
{
typedef typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType PromotedType;
@@ -879,7 +879,7 @@ namespace Dune {
{
typename V::ConstIterator e=this->end();
for (size_type i=0; i<y.n; i++)
- if (find(y.j[i])==e)
+ if (this->find(y.j[i])==e)
return false;
return true;
}
diff --git a/dune/istl/diagonalmatrix.hh b/dune/istl/diagonalmatrix.hh
index 128a0e391485915370b5dd8eba110bd7c6f0450b..533a690293d28450f39ecb340f859d7c2da56b84 100644
--- a/dune/istl/diagonalmatrix.hh
+++ b/dune/istl/diagonalmatrix.hh
@@ -1,1072 +1,10 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_DIAGONAL_MATRIX_HH
-#define DUNE_DIAGONAL_MATRIX_HH
-
/*! \file
\brief This file implements a quadratic diagonal matrix of fixed size.
+ \deprecated Use the corresponding file from dune-common instead!
*/
-#include <cmath>
-#include <cstddef>
-#include <complex>
-#include <iostream>
-#include <memory>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/typetraits.hh>
-#include <dune/common/genericiterator.hh>
-
-
-
-namespace Dune {
-
- template< class K, int n > class DiagonalRowVectorConst;
- template< class K, int n > class DiagonalRowVector;
- template< class DiagonalMatrixType > class DiagonalMatrixWrapper;
- template< class C, class T, class R> class ContainerWrapperIterator;
-
- /**
- @addtogroup DenseMatVec
- @{
- */
-
- /**
- *@brief A diagonal matrix of static size.
- *
- * This is meant to be a replacement of FieldMatrix for the case that
- * it is a diagonal matrix.
- *
- * \tparam K Type used for scalars
- * \tparam n Matrix size
- */
- template<class K, int n>
- class DiagonalMatrix
- {
- typedef DiagonalMatrixWrapper< DiagonalMatrix<K,n> > WrapperType;
-
- public:
- //===== type definitions and constants
-
- //! export the type representing the field
- typedef K field_type;
-
- //! export the type representing the components
- typedef K block_type;
-
- //! The type used for the index access and size operations.
- typedef std::size_t size_type;
-
- //! We are at the leaf of the block recursion
- enum {
- //! The number of block levels we contain. This is 1.
- blocklevel = 1
- };
-
- //! Each row is implemented by a field vector
- typedef DiagonalRowVector<K,n> row_type;
- typedef row_type reference;
- typedef DiagonalRowVectorConst<K,n> const_row_type;
- typedef const_row_type const_reference;
-
- //! export size
- enum {
- //! The number of rows
- rows = n,
- //! The number of columns
- cols = n
- };
-
-
-
- //===== constructors
-
- //! Default constructor
- DiagonalMatrix () {}
-
- //! Constructor initializing the whole matrix with a scalar
- DiagonalMatrix (const K& k)
- : diag_(k)
- {}
-
- //! Constructor initializing the diagonal with a vector
- DiagonalMatrix (const FieldVector<K,n>& diag)
- : diag_(diag)
- {}
-
-
- /** \brief Assignment from a scalar */
- DiagonalMatrix& operator= (const K& k)
- {
- diag_ = k;
- return *this;
- }
-
- /** \brief Check if matrix is the same object as the other matrix */
- bool identical(const DiagonalMatrix<K,n>& other) const
- {
- return (this==&other);
- }
-
- //===== iterator interface to rows of the matrix
- //! Iterator class for sequential access
- typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator;
- //! typedef for stl compliant access
- typedef Iterator iterator;
- //! rename the iterators for easier access
- typedef Iterator RowIterator;
- //! rename the iterators for easier access
- typedef typename row_type::Iterator ColIterator;
-
- //! begin iterator
- Iterator begin ()
- {
- return Iterator(WrapperType(this),0);
- }
-
- //! end iterator
- Iterator end ()
- {
- return Iterator(WrapperType(this),n);
- }
-
- //! @returns an iterator that is positioned before
- //! the end iterator of the rows, i.e. at the last row.
- Iterator beforeEnd ()
- {
- return Iterator(WrapperType(this),n-1);
- }
-
- //! @returns an iterator that is positioned before
- //! the first row of the matrix.
- Iterator beforeBegin ()
- {
- return Iterator(WrapperType(this),-1);
- }
-
-
- //! Iterator class for sequential access
- typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator;
- //! typedef for stl compliant access
- typedef ConstIterator const_iterator;
- //! rename the iterators for easier access
- typedef ConstIterator ConstRowIterator;
- //! rename the iterators for easier access
- typedef typename const_row_type::ConstIterator ConstColIterator;
-
- //! begin iterator
- ConstIterator begin () const
- {
- return ConstIterator(WrapperType(this),0);
- }
-
- //! end iterator
- ConstIterator end () const
- {
- return ConstIterator(WrapperType(this),n);
- }
-
- //! @returns an iterator that is positioned before
- //! the end iterator of the rows. i.e. at the last row.
- ConstIterator beforeEnd() const
- {
- return ConstIterator(WrapperType(this),n-1);
- }
-
- //! @returns an iterator that is positioned before
- //! the first row of the matrix.
- ConstIterator beforeBegin () const
- {
- return ConstIterator(WrapperType(this),-1);
- }
-
-
-
- //===== vector space arithmetic
-
- //! vector space addition
- DiagonalMatrix& operator+= (const DiagonalMatrix& y)
- {
- diag_ += y.diag_;
- return *this;
- }
-
- //! vector space subtraction
- DiagonalMatrix& operator-= (const DiagonalMatrix& y)
- {
- diag_ -= y.diag_;
- return *this;
- }
-
- //! vector space multiplication with scalar
- DiagonalMatrix& operator+= (const K& k)
- {
- diag_ += k;
- return *this;
- }
-
- //! vector space division by scalar
- DiagonalMatrix& operator-= (const K& k)
- {
- diag_ -= k;
- return *this;
- }
-
- //! vector space multiplication with scalar
- DiagonalMatrix& operator*= (const K& k)
- {
- diag_ *= k;
- return *this;
- }
-
- //! vector space division by scalar
- DiagonalMatrix& operator/= (const K& k)
- {
- diag_ /= k;
- return *this;
- }
-
- //===== comparison ops
-
- //! comparison operator
- bool operator==(const DiagonalMatrix& other) const
- {
- return diag_==other.diagonal();
- }
-
- //! incomparison operator
- bool operator!=(const DiagonalMatrix& other) const
- {
- return diag_!=other.diagonal();
- }
-
-
- //===== linear maps
-
- //! y = A x
- template<class X, class Y>
- void mv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; ++i)
- y[i] = diag_[i] * x[i];
- }
-
- //! y = A^T x
- template<class X, class Y>
- void mtv (const X& x, Y& y) const
- {
- mv(x, y);
- }
-
- //! y += A x
- template<class X, class Y>
- void umv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; ++i)
- y[i] += diag_[i] * x[i];
- }
-
- //! y += A^T x
- template<class X, class Y>
- void umtv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; ++i)
- y[i] += diag_[i] * x[i];
- }
-
- //! y += A^H x
- template<class X, class Y>
- void umhv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; i++)
- y[i] += conjugateComplex(diag_[i])*x[i];
- }
-
- //! y -= A x
- template<class X, class Y>
- void mmv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; ++i)
- y[i] -= diag_[i] * x[i];
- }
-
- //! y -= A^T x
- template<class X, class Y>
- void mmtv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; ++i)
- y[i] -= diag_[i] * x[i];
- }
-
- //! y -= A^H x
- template<class X, class Y>
- void mmhv (const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; i++)
- y[i] -= conjugateComplex(diag_[i])*x[i];
- }
-
- //! y += alpha A x
- template<class X, class Y>
- void usmv (const K& alpha, const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; i++)
- y[i] += alpha * diag_[i] * x[i];
- }
-
- //! y += alpha A^T x
- template<class X, class Y>
- void usmtv (const K& alpha, const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; i++)
- y[i] += alpha * diag_[i] * x[i];
- }
-
- //! y += alpha A^H x
- template<class X, class Y>
- void usmhv (const K& alpha, const X& x, Y& y) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
- if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
-#endif
- for (size_type i=0; i<n; i++)
- y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
- }
-
- //===== norms
-
- //! frobenius norm: sqrt(sum over squared values of entries)
- double frobenius_norm () const
- {
- return diag_.two_norm();
- }
-
- //! square of frobenius norm, need for block recursion
- double frobenius_norm2 () const
- {
- return diag_.two_norm2();
- }
-
- //! infinity norm (row sum norm, how to generalize for blocks?)
- double infinity_norm () const
- {
- return diag_.infinity_norm();
- }
-
- //! simplified infinity norm (uses Manhattan norm for complex values)
- double infinity_norm_real () const
- {
- return diag_.infinity_norm_real();
- }
-
-
-
- //===== solve
-
- //! Solve system A x = b
- template<class V>
- void solve (V& x, const V& b) const
- {
- for (int i=0; i<n; i++)
- x[i] = b[i]/diag_[i];
- }
-
- //! Compute inverse
- void invert()
- {
- for (int i=0; i<n; i++)
- diag_[i] = 1/diag_[i];
- }
-
- //! calculates the determinant of this matrix
- K determinant () const
- {
- K det = diag_[0];
- for (int i=1; i<n; i++)
- det *= diag_[i];
- return det;
- }
-
-
-
- //===== sizes
-
- //! number of blocks in row direction
- size_type N () const
- {
- return n;
- }
-
- //! number of blocks in column direction
- size_type M () const
- {
- return n;
- }
-
-
-
- //===== query
-
- //! return true when (i,j) is in pattern
- bool exists (size_type i, size_type j) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range");
- if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range");
-#endif
- return i==j;
- }
-
-
-
- //! Sends the matrix to an output stream
- friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
- {
- for (size_type i=0; i<n; i++) {
- for (size_type j=0; j<n; j++)
- s << ((i==j) ? a.diag_[i] : 0) << " ";
- s << std::endl;
- }
- return s;
- }
-
- //! Return reference object as row replacement
- reference operator[](size_type i)
- {
- return reference(const_cast<K*>(&diag_[i]), i);
- }
-
- //! Return const_reference object as row replacement
- const_reference operator[](size_type i) const
- {
- return const_reference(const_cast<K*>(&diag_[i]), i);
- }
-
- //! Get const reference to diagonal entry
- const K& diagonal(size_type i) const
- {
- return diag_[i];
- }
-
- //! Get reference to diagonal entry
- K& diagonal(size_type i)
- {
- return diag_[i];
- }
-
- //! Get const reference to diagonal vector
- const FieldVector<K,n>& diagonal() const
- {
- return diag_;
- }
-
- //! Get reference to diagonal vector
- FieldVector<K,n>& diagonal()
- {
- return diag_;
- }
-
- private:
-
- // the data, a FieldVector storing the diagonal
- FieldVector<K,n> diag_;
- };
-
-#ifndef DOXYGEN // hide specialization
- /** \brief Special type for 1x1 matrices
- */
- template< class K >
- class DiagonalMatrix<K, 1> : public FieldMatrix<K, 1, 1>
- {
- typedef FieldMatrix<K,1,1> Base;
- public:
- //! The type used for index access and size operations
- typedef typename Base::size_type size_type;
-
- //! We are at the leaf of the block recursion
- enum {
- //! The number of block levels we contain.
- //! This is always one for this type.
- blocklevel = 1
- };
-
- typedef typename Base::row_type row_type;
-
- typedef typename Base::row_reference row_reference;
- typedef typename Base::const_row_reference const_row_reference;
-
- //! export size
- enum {
- //! \brief The number of rows.
- //! This is always one for this type.
- rows = 1,
- //! \brief The number of columns.
- //! This is always one for this type.
- cols = 1
- };
-
-
- //! Default Constructor
- DiagonalMatrix()
- {}
-
- //! Constructor initializing the whole matrix with a scalar
- DiagonalMatrix(const K& scalar)
- {
- (*this)[0][0] = scalar;
- }
-
- //! Constructor initializing the whole matrix with a scalar from convertible type
- template<typename T>
- DiagonalMatrix(const T& t)
- {
- DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
- }
-
- //! Get const reference to diagonal entry
- const K& diagonal(size_type) const
- {
- return (*this)[0][0];
- }
-
- //! Get reference to diagonal entry
- K& diagonal(size_type)
- {
- return (*this)[0][0];
- }
-
- //! Get const reference to diagonal vector
- const FieldVector<K,1>& diagonal() const
- {
- return (*this)[0];
- }
-
- //! Get reference to diagonal vector
- FieldVector<K,1>& diagonal()
- {
- return (*this)[0];
- }
-
- };
-#endif
-
-
- template<class DiagonalMatrixType>
- class DiagonalMatrixWrapper
- {
- typedef typename DiagonalMatrixType::reference reference;
- typedef typename DiagonalMatrixType::const_reference const_reference;
- typedef typename DiagonalMatrixType::field_type K;
- typedef DiagonalRowVector<K, DiagonalMatrixType::rows> row_type;
- typedef std::size_t size_type;
- typedef DiagonalMatrixWrapper< DiagonalMatrixType> MyType;
-
- friend class ContainerWrapperIterator<const MyType, reference, reference>;
- friend class ContainerWrapperIterator<const MyType, const_reference, const_reference>;
-
- public:
-
- DiagonalMatrixWrapper() :
- mat_(0)
- {}
-
- DiagonalMatrixWrapper(const DiagonalMatrixType* mat) :
- mat_(const_cast<DiagonalMatrixType*>(mat))
- {}
-
- size_type realIndex(int i) const
- {
- return i;
- }
-
- row_type* pointer(int i) const
- {
- row_ = row_type(&(mat_->diagonal(i)), i);
- return &row_;
- }
-
- bool identical(const DiagonalMatrixWrapper& other) const
- {
- return mat_==other.mat_;
- }
-
- private:
-
- mutable DiagonalMatrixType* mat_;
- mutable row_type row_;
- };
-
- /** \brief
- *
- */
- template< class K, int n >
- class DiagonalRowVectorConst
- {
- template<class DiagonalMatrixType>
- friend class DiagonalMatrixWrapper;
- friend class ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&>;
-
- public:
- // remember size of vector
- enum { dimension = n };
-
- // standard constructor and everything is sufficient ...
-
- //===== type definitions and constants
-
- //! export the type representing the field
- typedef K field_type;
-
- //! export the type representing the components
- typedef K block_type;
-
- //! The type used for the index access and size operation
- typedef std::size_t size_type;
-
- //! We are at the leaf of the block recursion
- enum {
- //! The number of block levels we contain
- blocklevel = 1
- };
-
- //! export size
- enum {
- //! The size of this vector.
- size = n
- };
-
- //! Constructor making uninitialized vector
- DiagonalRowVectorConst() :
- p_(0),
- row_(0)
- {}
-
- //! Constructor making vector with identical coordinates
- explicit DiagonalRowVectorConst (K* p, int col) :
- p_(p),
- row_(col)
- {}
-
- //===== access to components
-
- //! same for read only access
- const K& operator[] (size_type i) const
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (i!=row_)
- DUNE_THROW(FMatrixError,"index is not contained in pattern");
-#endif
- return *p_;
- }
-
- // check if row is identical to other row (not only identical values)
- // since this is a proxy class we need to check equality of the stored pointer
- bool identical(const DiagonalRowVectorConst<K,n>& other) const
- {
- return ((p_ == other.p_)and (row_ == other.row_));
- }
-
- //! ConstIterator class for sequential access
- typedef ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&> ConstIterator;
- //! typedef for stl compliant access
- typedef ConstIterator const_iterator;
-
- //! begin ConstIterator
- ConstIterator begin () const
- {
- return ConstIterator(*this,0);
- }
-
- //! end ConstIterator
- ConstIterator end () const
- {
- return ConstIterator(*this,1);
- }
-
- //! @returns an iterator that is positioned before
- //! the end iterator of the rows. i.e. at the row.
- ConstIterator beforeEnd() const
- {
- return ConstIterator(*this,0);
- }
-
- //! @returns an iterator that is positioned before
- //! the first row of the matrix.
- ConstIterator beforeBegin () const
- {
- return ConstIterator(*this,-1);
- }
-
- //! Binary vector comparison
- bool operator== (const DiagonalRowVectorConst& y) const
- {
- return ((p_==y.p_)and (row_==y.row_));
- }
-
- //===== sizes
-
- //! number of blocks in the vector (are of size 1 here)
- size_type N () const
- {
- return n;
- }
-
- //! dimension of the vector space
- size_type dim () const
- {
- return n;
- }
-
- //! index of this row in surrounding matrix
- size_type rowIndex() const
- {
- return row_;
- }
-
- //! the diagonal value
- const K& diagonal() const
- {
- return *p_;
- }
-
- protected:
-
- size_type realIndex(int i) const
- {
- return rowIndex();
- }
-
- K* pointer(size_type i) const
- {
- return const_cast<K*>(p_);
- }
-
- DiagonalRowVectorConst* operator&()
- {
- return this;
- }
-
- // the data, very simply a pointer to the diagonal value and the row number
- K* p_;
- size_type row_;
- };
-
- template< class K, int n >
- class DiagonalRowVector : public DiagonalRowVectorConst<K,n>
- {
- template<class DiagonalMatrixType>
- friend class DiagonalMatrixWrapper;
- friend class ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&>;
-
- public:
- // standard constructor and everything is sufficient ...
-
- //===== type definitions and constants
-
- //! export the type representing the field
- typedef K field_type;
-
- //! export the type representing the components
- typedef K block_type;
-
- //! The type used for the index access and size operation
- typedef std::size_t size_type;
-
- //! Constructor making uninitialized vector
- DiagonalRowVector() : DiagonalRowVectorConst<K,n>()
- {}
-
- //! Constructor making vector with identical coordinates
- explicit DiagonalRowVector (K* p, int col) : DiagonalRowVectorConst<K,n>(p, col)
- {}
-
- //===== assignment from scalar
- //! Assignment operator for scalar
- DiagonalRowVector& operator= (const K& k)
- {
- *p_ = k;
- return *this;
- }
-
- //===== access to components
-
- //! random access
- K& operator[] (size_type i)
- {
-#ifdef DUNE_FMatrix_WITH_CHECKING
- if (i!=row_)
- DUNE_THROW(FMatrixError,"index is contained in pattern");
-#endif
- return *p_;
- }
-
- //! Iterator class for sequential access
- typedef ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&> Iterator;
- //! typedef for stl compliant access
- typedef Iterator iterator;
-
- //! begin iterator
- Iterator begin ()
- {
- return Iterator(*this, 0);
- }
-
- //! end iterator
- Iterator end ()
- {
- return Iterator(*this, 1);
- }
-
- //! @returns an iterator that is positioned before
- //! the end iterator of the rows, i.e. at the last row.
- Iterator beforeEnd ()
- {
- return Iterator(*this, 0);
- }
-
- //! @returns an iterator that is positioned before
- //! the first row of the matrix.
- Iterator beforeBegin ()
- {
- return Iterator(*this, -1);
- }
-
- //! ConstIterator class for sequential access
- typedef ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&> ConstIterator;
- //! typedef for stl compliant access
- typedef ConstIterator const_iterator;
-
- using DiagonalRowVectorConst<K,n>::identical;
- using DiagonalRowVectorConst<K,n>::operator[];
- using DiagonalRowVectorConst<K,n>::operator==;
- using DiagonalRowVectorConst<K,n>::begin;
- using DiagonalRowVectorConst<K,n>::end;
- using DiagonalRowVectorConst<K,n>::beforeEnd;
- using DiagonalRowVectorConst<K,n>::beforeBegin;
- using DiagonalRowVectorConst<K,n>::N;
- using DiagonalRowVectorConst<K,n>::dim;
- using DiagonalRowVectorConst<K,n>::rowIndex;
- using DiagonalRowVectorConst<K,n>::diagonal;
-
- protected:
-
- DiagonalRowVector* operator&()
- {
- return this;
- }
-
- private:
-
- using DiagonalRowVectorConst<K,n>::p_;
- using DiagonalRowVectorConst<K,n>::row_;
- };
-
-
- // implement type traits
- template<class K, int n>
- struct const_reference< DiagonalRowVector<K,n> >
- {
- typedef DiagonalRowVectorConst<K,n> type;
- };
-
- template<class K, int n>
- struct const_reference< DiagonalRowVectorConst<K,n> >
- {
- typedef DiagonalRowVectorConst<K,n> type;
- };
-
- template<class K, int n>
- struct mutable_reference< DiagonalRowVector<K,n> >
- {
- typedef DiagonalRowVector<K,n> type;
- };
-
- template<class K, int n>
- struct mutable_reference< DiagonalRowVectorConst<K,n> >
- {
- typedef DiagonalRowVector<K,n> type;
- };
-
-
-
- /** \brief Iterator class for sparse vector-like containers
- *
- * This class provides an iterator for sparse vector like containers.
- * It contains a ContainerWrapper that must provide the translation
- * from the position in the underlying container to the index
- * in the sparse container.
- *
- * The ContainerWrapper must be default and copy-constructable.
- * Furthermore it must provide the methods:
- *
- * bool identical(other) - check if this is identical to other (same container, not only equal)
- * T* pointer(position) - get pointer to data at position in underlying container
- * size_t realIndex(position) - get index in sparse container for position in underlying container
- *
- * Notice that the iterator stores a ContainerWrapper.
- * This allows to use proxy classes as underlying container
- * and as returned reference type.
- *
- * \tparam CW The container wrapper class
- * \tparam T The contained type
- * \tparam R The reference type returned by dereference
- */
- template<class CW, class T, class R>
- class ContainerWrapperIterator : public BidirectionalIteratorFacade<ContainerWrapperIterator<CW,T,R>,T, R, int>
- {
- typedef typename remove_const<CW>::type NonConstCW;
-
- friend class ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type>;
- friend class ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type>;
-
- typedef ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type> MyType;
- typedef ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type> MyConstType;
-
- public:
-
- // Constructors needed by the facade iterators.
- ContainerWrapperIterator() :
- containerWrapper_(),
- position_(0)
- {}
-
- ContainerWrapperIterator(CW containerWrapper, int position) :
- containerWrapper_(containerWrapper),
- position_(position)
- {}
-
- template<class OtherContainerWrapperIteratorType>
- ContainerWrapperIterator(OtherContainerWrapperIteratorType& other) :
- containerWrapper_(other.containerWrapper_),
- position_(other.position_)
- {}
-
- ContainerWrapperIterator(const MyType& other) :
- containerWrapper_(other.containerWrapper_),
- position_(other.position_)
- {}
-
- ContainerWrapperIterator(const MyConstType& other) :
- containerWrapper_(other.containerWrapper_),
- position_(other.position_)
- {}
-
- template<class OtherContainerWrapperIteratorType>
- ContainerWrapperIterator& operator=(OtherContainerWrapperIteratorType& other)
- {
- containerWrapper_ = other.containerWrapper_;
- position_ = other.position_;
- }
-
- // This operator is needed since we can not get the address of the
- // temporary object returned by dereference
- T* operator->() const
- {
- return containerWrapper_.pointer(position_);
- }
-
- // Methods needed by the forward iterator
- bool equals(const MyType& other) const
- {
- return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
- }
-
- bool equals(const MyConstType& other) const
- {
- return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
- }
-
- R dereference() const
- {
- return *containerWrapper_.pointer(position_);
- }
-
- void increment()
- {
- ++position_;
- }
-
- // Additional function needed by BidirectionalIterator
- void decrement()
- {
- --position_;
- }
-
- // Additional function needed by RandomAccessIterator
- R elementAt(int i) const
- {
- return *containerWrapper_.pointer(position_+i);
- }
-
- void advance(int n)
- {
- position_=position_+n;
- }
-
- template<class OtherContainerWrapperIteratorType>
- std::ptrdiff_t distanceTo(OtherContainerWrapperIteratorType& other) const
- {
- assert(containerWrapper_.identical(other));
- return other.position_ - position_;
- }
-
- std::ptrdiff_t index() const
- {
- return containerWrapper_.realIndex(position_);
- }
-
- private:
- NonConstCW containerWrapper_;
- size_t position_;
- };
-
-
-
- template<class M, class K, int n>
- void istl_assign_to_fmatrix(DenseMatrix<M>& fm, const DiagonalMatrix<K,n>& s)
- {
- fm = K();
- for(int i=0; i<n; ++i)
- fm[i][i] = s.diagonal()[i];
- }
- /* @} */
-} // end namespace
-#endif
+#warning The file diagonalmatrix.hh has moved from dune-istl to dune-common,\
+ and will be removed from dune-istl after dune 2.3!
+#include <dune/common/diagonalmatrix.hh>
diff --git a/dune/istl/io.hh b/dune/istl/io.hh
index 85b6644a949826794aa8b090822e93a7a3a320b7..d44c8e914ef1171c4b737a54d4c371af926e60a9 100644
--- a/dune/istl/io.hh
+++ b/dune/istl/io.hh
@@ -15,9 +15,9 @@
#include "istlexception.hh"
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
+#include <dune/common/diagonalmatrix.hh>
#include <dune/istl/matrix.hh>
-#include <dune/istl/diagonalmatrix.hh>
#include <dune/istl/scaledidmatrix.hh>
#include "bcrsmatrix.hh"
diff --git a/dune/istl/matrixutils.hh b/dune/istl/matrixutils.hh
index a4e2fb7ff628cc909ec5ee3a936698ac27831242..28f29ae78b71e8b4772fe2fb3c39eb651ee1098f 100644
--- a/dune/istl/matrixutils.hh
+++ b/dune/istl/matrixutils.hh
@@ -9,7 +9,7 @@
#include <dune/common/typetraits.hh>
#include <dune/common/static_assert.hh>
#include <dune/common/fmatrix.hh>
-#include <dune/istl/diagonalmatrix.hh>
+#include <dune/common/diagonalmatrix.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include "istlexception.hh"
diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 425211b7329210f36f8a0ff6d66f1e520ea9962d..78e49e01d15ddc171dd665d28da02426542186e6 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -397,13 +397,13 @@ namespace Dune
typedef typename Matrix::ConstColIterator ColIter;
typedef typename Matrix::block_type Block;
Block zero;
- Block diagonal;
zero=typename Matrix::field_type();
const Matrix& mat=matrices_->matrices().finest()->getmat();
for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
bool isDirichlet = true;
bool hasDiagonal = false;
+ Block diagonal;
for(ColIter col=row->begin(); col!=row->end(); ++col) {
if(row.index()==col.index()) {
diagonal = *col;
@@ -413,7 +413,7 @@ namespace Dune
isDirichlet = false;
}
}
- if(isDirichlet)
+ if(isDirichlet && hasDiagonal)
diagonal.solve(x[row.index()], b[row.index()]);
}
diff --git a/dune/istl/scaledidmatrix.hh b/dune/istl/scaledidmatrix.hh
index 054f89540e28d2875b433ef3f254fb2a5550a8a5..162ef4777f8fbeb3052552e66deb15e79a855050 100644
--- a/dune/istl/scaledidmatrix.hh
+++ b/dune/istl/scaledidmatrix.hh
@@ -15,7 +15,7 @@
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
-#include <dune/istl/diagonalmatrix.hh>
+#include <dune/common/diagonalmatrix.hh>
namespace Dune {
diff --git a/dune/istl/test/.gitignore b/dune/istl/test/.gitignore
index 3614259ef262aeb760e2f7382af3b96af0a47d7f..a0787d050b89d48b2aee8c920d3c0dc199da1f1d 100644
--- a/dune/istl/test/.gitignore
+++ b/dune/istl/test/.gitignore
@@ -4,7 +4,6 @@ Makefile.in
.libs
semantic.cache
complexrhstest
-diagonalmatrixtest
dotproducttest
bvectortest
matrixutilstest
diff --git a/dune/istl/test/Makefile.am b/dune/istl/test/Makefile.am
index cc3d48b4ce79ef8f13eca80a7a17af1db7394404..bad932b37ae77056579dbd8f8b47be374f67c1c2 100644
--- a/dune/istl/test/Makefile.am
+++ b/dune/istl/test/Makefile.am
@@ -21,7 +21,6 @@ NORMALTESTS = basearraytest \
bcrsbuildtest \
bvectortest \
complexrhstest \
- diagonalmatrixtest \
dotproducttest \
iotest \
matrixiteratortest \
@@ -88,8 +87,6 @@ complexrhstest_LDADD= $(SUPERLU_LIBS)
complexrhstest_LDFLAGS= $(AM_LDFLAGS) $(SUPERLU_LDFLAGS)
complexrhstest_CPPFLAGS=$(AM_CPPFLAGS) $(SUPERLU_CPPFLAGS) -DSUPERLU_NTYPE=3
-diagonalmatrixtest_SOURCES = diagonalmatrixtest.cc
-
dotproducttest_SOURCES = dotproducttest.cc
vbvectortest_SOURCES = vbvectortest.cc
diff --git a/dune/istl/test/diagonalmatrixtest.cc b/dune/istl/test/diagonalmatrixtest.cc
deleted file mode 100644
index 92247b28f899a7908835524c53eadf99adc184f3..0000000000000000000000000000000000000000
--- a/dune/istl/test/diagonalmatrixtest.cc
+++ /dev/null
@@ -1,64 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#include "config.h"
-
-#include <dune/istl/diagonalmatrix.hh>
-
-#include <iostream>
-#include <algorithm>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/exceptions.hh>
-
-using namespace Dune;
-
-
-template<class K, int n>
-void test_matrix()
-{
- typedef typename DiagonalMatrix<K,n>::size_type size_type;
-
- DiagonalMatrix<K,n> A(1);
- FieldVector<K,n> f;
- FieldVector<K,n> v;
-
- // assign matrix
- A=2;
-
- // assign vector
- f = 1;
- v = 2;
-
- // matrix vector product
- A.umv(v,f);
-
-
- // test norms
- A.frobenius_norm();
- A.frobenius_norm2();
- A.infinity_norm();
- A.infinity_norm_real();
-
- std::sort(v.begin(), v.end());
-
- // print matrix
- std::cout << A << std::endl;
- // print vector
- std::cout << f << std::endl;
-
- // assign to FieldMatrix
- FieldMatrix<K,n,n> AFM = FieldMatrix<K,n,n>(A);
-}
-
-int main()
-{
- try {
- test_matrix<float, 1>();
- test_matrix<double, 1>();
- test_matrix<double, 5>();
- }
- catch (Dune::Exception & e)
- {
- std::cerr << "Exception: " << e << std::endl;
- }
-}
diff --git a/dune/istl/test/iotest.cc b/dune/istl/test/iotest.cc
index 7a7de78f8686a91c68c9888478c31aa83ebb544c..d2d63601399ad227a2f40456c7175dc75d4a6179 100644
--- a/dune/istl/test/iotest.cc
+++ b/dune/istl/test/iotest.cc
@@ -2,7 +2,7 @@
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
#include <dune/common/fmatrix.hh>
-#include <dune/istl/diagonalmatrix.hh>
+#include <dune/common/diagonalmatrix.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/io.hh>
diff --git a/dune/istl/test/matrixtest.cc b/dune/istl/test/matrixtest.cc
index e71743cf7960cb0788537ed09351700fd8a67b95..3a06aaf9f0326c2f518cc0f5dc7bf8dc0c301621 100644
--- a/dune/istl/test/matrixtest.cc
+++ b/dune/istl/test/matrixtest.cc
@@ -8,12 +8,12 @@
#include <fenv.h>
#include <dune/common/fmatrix.hh>
+#include <dune/common/diagonalmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bdmatrix.hh>
#include <dune/istl/btdmatrix.hh>
#include <dune/istl/scaledidmatrix.hh>
-#include <dune/istl/diagonalmatrix.hh>
using namespace Dune;
@@ -197,7 +197,7 @@ void testMatrix(MatrixType& matrix, X& x, Y& y)
matrix = 0;
// The copy constructor
- MatrixType thirdMatrix(matrix);
+ DUNE_UNUSED MatrixType thirdMatrix(matrix);
// ///////////////////////////////////////////////////////
// Test component-wise operations