From 38d22f814b19f91b0cf696a874c2b3b8229d7022 Mon Sep 17 00:00:00 2001
From: Uli Sack <usack@math.fu-berlin.de>
Date: Tue, 10 Jan 2012 15:36:49 +0000
Subject: [PATCH] now, here are the classes and tests
[[Imported from SVN: r5238]]
---
dune/solvers/operators/Makefile.am | 8 +
dune/solvers/operators/lowrankoperator.hh | 57 +++++++
dune/solvers/operators/nulloperator.hh | 87 +++++++++++
dune/solvers/operators/sumoperator.hh | 62 ++++++++
dune/solvers/test/Makefile.am | 38 +++++
dune/solvers/test/lowrankoperatortest.cc | 172 ++++++++++++++++++++++
dune/solvers/test/nulloperatortest.cc | 104 +++++++++++++
dune/solvers/test/sumoperatortest.cc | 157 ++++++++++++++++++++
8 files changed, 685 insertions(+)
create mode 100644 dune/solvers/operators/Makefile.am
create mode 100644 dune/solvers/operators/lowrankoperator.hh
create mode 100644 dune/solvers/operators/nulloperator.hh
create mode 100644 dune/solvers/operators/sumoperator.hh
create mode 100644 dune/solvers/test/Makefile.am
create mode 100644 dune/solvers/test/lowrankoperatortest.cc
create mode 100644 dune/solvers/test/nulloperatortest.cc
create mode 100644 dune/solvers/test/sumoperatortest.cc
diff --git a/dune/solvers/operators/Makefile.am b/dune/solvers/operators/Makefile.am
new file mode 100644
index 00000000..e1c13db9
--- /dev/null
+++ b/dune/solvers/operators/Makefile.am
@@ -0,0 +1,8 @@
+SUBDIRS =
+
+operatorsdir = $(includedir)/dune/solvers/operators
+operators_HEADERS = lowrankoperator.hh \
+ nulloperator.hh \
+ sumoperator.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/solvers/operators/lowrankoperator.hh b/dune/solvers/operators/lowrankoperator.hh
new file mode 100644
index 00000000..916112d3
--- /dev/null
+++ b/dune/solvers/operators/lowrankoperator.hh
@@ -0,0 +1,57 @@
+#ifndef LOWRANKOPERATOR_HH
+#define LOWRANKOPERATOR_HH
+
+
+/** \brief Class that behaves like a filled in low-rank matrix defined via a factor \f$ M\in\mathbb R^{k\times n}, k<<n \f$
+ *
+ * The full matrix is given by \f$ A=M^TM\f$
+ *
+ * \tparam LRFactorType The (Matrix-) type of the defining factor M
+ */
+template <class LRFactorType>
+class LowRankOperator
+{
+ public:
+ //! export the type of the defining factor
+ typedef LRFactorType LowRankFactorType;
+
+ //! constructor taking the defining factor as argument
+ LowRankOperator(LowRankFactorType& lrFactor):
+ lowRankFactor_(lrFactor)
+ {}
+
+ //! b += Ax
+ template <class LVectorType, class RVectorType>
+ void umv(const LVectorType& x, RVectorType& b) const
+ {
+ LVectorType temp(lowRankFactor_.N());
+ lowRankFactor_.mv(x,temp);
+
+ for (size_t i=0; i<b.size(); ++i)
+ for (size_t j=0; j<lowRankFactor_.N(); ++j)
+ if (lowRankFactor_.exists(j,i))
+ lowRankFactor_[j][i].umv(temp[j],b[i]);
+ }
+
+ //! b = Ax
+ template <class LVectorType, class RVectorType>
+ void mv(const LVectorType& x, RVectorType& b) const
+ {
+ b = 0.0;
+ umv(x,b);
+ }
+
+ //! returns the defining factor
+ const LowRankFactorType& getLowRankFactor() const
+ {
+ return lowRankFactor_;
+ }
+
+ private:
+ //! the defining factor
+ LowRankFactorType& lowRankFactor_;
+
+};
+
+#endif
+
diff --git a/dune/solvers/operators/nulloperator.hh b/dune/solvers/operators/nulloperator.hh
new file mode 100644
index 00000000..f9b2888b
--- /dev/null
+++ b/dune/solvers/operators/nulloperator.hh
@@ -0,0 +1,87 @@
+#ifndef NULLOPERATOR_HH
+#define NULLOPERATOR_HH
+
+/** \brief Represents the null operator in the needed space
+ *
+ */
+template <class BlockType>
+class NullOperator
+{
+ private:
+ /** \brief A dummy class for the rows
+ *
+ * Contains a zero block to refer to
+ */
+ class RowDummy
+ {
+ private:
+ BlockType zero_;
+
+ public:
+ RowDummy(): zero_(0){};
+
+ const BlockType& operator[](size_t i) const
+ {
+ return zero_;
+ }
+ };
+
+ //! a dummy row
+ RowDummy rowDummy_;
+
+ public:
+ //! export the block type
+ typedef BlockType block_type;
+
+ //! Default constructor
+ NullOperator(){}
+
+ /** \brief constructor taking anything
+ *
+ * This is here to allow for NullOperator as block_type (this is needed in the constructor of RowDummy)
+ */
+ template <class T>
+ NullOperator(const T& t){}
+
+ /** \brief Matrix-Vector multiplication
+ *
+ * Implements b += Nx and hence does nothing (N=0 !)
+ */
+ template <class LVectorType, class RVectorType>
+ void umv(const LVectorType& x, RVectorType& b) const
+ {}
+
+ /** \brief Matrix-Vector multiplication with scalar multiplication
+ *
+ * Implements b += a*Nx and hence does nothing (N=0 !)
+ */
+ template <class LVectorType, class RVectorType>
+ void usmv(const double a, const LVectorType& x, RVectorType& b) const
+ {}
+
+ /** \brief Matrix-Vector multiplication
+ *
+ * Implements b = Nx and hence does nothing but set b=0 (N=0 !)
+ */
+ template <class LVectorType, class RVectorType>
+ void mv(const LVectorType& x, RVectorType& b) const
+ {
+ b = 0.0;
+ }
+
+ //! random access operator
+ const RowDummy& operator[](size_t i) const
+ {
+ return rowDummy_;
+ }
+
+ //! return j-th diagonal entry
+ const block_type& diagonal(size_t i) const
+ {
+ return rowDummy_[0];
+ }
+
+};
+
+#endif
+
diff --git a/dune/solvers/operators/sumoperator.hh b/dune/solvers/operators/sumoperator.hh
new file mode 100644
index 00000000..7869d376
--- /dev/null
+++ b/dune/solvers/operators/sumoperator.hh
@@ -0,0 +1,62 @@
+#ifndef SUMOPERATOR_HH
+#define SUMOPERATOR_HH
+
+/** \brief represents the sum of two linear operators
+ *
+ * The summands may be of arbitrary type. Their names are for historical illustrative reasons where one is a sparse and the other a filled in lowrank matrix
+ * \tparam SparseMatrix the type of one of the summands
+ * \tparam LowRankMatrix the type of the other summand
+ */
+template <class SparseMatrix, class LowRankMatrix>
+class SumOperator
+{
+ public:
+ //! export the summand type
+ typedef SparseMatrix SparseMatrixType;
+ //! export the summand type
+ typedef LowRankMatrix LowRankMatrixType;
+
+ //! construct from summands
+ SumOperator(SparseMatrixType& A, LowRankMatrixType& M):
+ sparse_matrix_(A),
+ lowrank_matrix_(M)
+ {}
+
+ //! b += (A+M)x
+ template <class LVectorType, class RVectorType>
+ void umv(const LVectorType& x, RVectorType& b) const
+ {
+ sparse_matrix_.umv(x,b);
+ lowrank_matrix_.umv(x,b);
+ }
+
+ //! b = (A+M)x
+ template <class LVectorType, class RVectorType>
+ void mv(const LVectorType& x, RVectorType& b)
+ {
+ sparse_matrix_.mv(x,b);
+ lowrank_matrix_.umv(x,b);
+ }
+
+ //! return one summand
+ SparseMatrixType& getSparseMatrix()
+ {
+ return sparse_matrix_;
+ }
+
+ //! return "the other" summand
+ LowRankMatrixType& getLowRankMatrix()
+ {
+ return lowrank_matrix_;
+ }
+
+ private:
+ //! one summand
+ SparseMatrixType& sparse_matrix_;
+ //! the other summand
+ LowRankMatrixType& lowrank_matrix_;
+
+};
+
+#endif
+
diff --git a/dune/solvers/test/Makefile.am b/dune/solvers/test/Makefile.am
new file mode 100644
index 00000000..df9e99f2
--- /dev/null
+++ b/dune/solvers/test/Makefile.am
@@ -0,0 +1,38 @@
+
+# list of tests to run
+TESTS = lowrankoperatortest \
+ nulloperatortest \
+ sumoperatortest
+
+# programs just to build when "make check" is used
+check_PROGRAMS = $(TESTS)
+
+noinst_HEADERS =
+
+# collect some common flags
+COMMON_CPPFLAGS = $(AM_CPPFLAGS) $(DUNEMPICPPFLAGS)
+COMMON_LDADD = $(DUNE_LDFLAGS) $(DUNE_LIBS) $(DUNEMPILIBS) $(LDADD)
+COMMON_LDFLAGS = $(AM_LDFLAGS) $(DUNEMPILDFLAGS) $(DUNE_LDFLAGS)
+
+GRID_CPPFLAGS = $(UG_CPPFLAGS) $(ALUGRID_CPPFLAGS)
+GRID_LDADD = $(UG_LDFLAGS) $(UG_LIBS) $(ALUGRID_LDFLAGS) $(ALUGRID_LIBS)
+GRID_LDFLAGS = $(UG_LDFLAGS) $(ALUGRID_LDFLAGS)
+
+
+# define the programs
+lowrankoperatortest_SOURCES = lowrankoperatortest.cc
+lowrankoperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
+lowrankoperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
+lowrankoperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
+
+nulloperatortest_SOURCES = nulloperatortest.cc
+nulloperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
+nulloperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
+nulloperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
+
+sumoperatortest_SOURCES = sumoperatortest.cc
+sumoperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
+sumoperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
+sumoperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/solvers/test/lowrankoperatortest.cc b/dune/solvers/test/lowrankoperatortest.cc
new file mode 100644
index 00000000..8a456e29
--- /dev/null
+++ b/dune/solvers/test/lowrankoperatortest.cc
@@ -0,0 +1,172 @@
+#include <config.h>
+
+#include <stdio.h>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/istl/diagonalmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/operators/lowrankoperator.hh>
+
+#define BLOCKSIZE 2
+
+// This tests the LowRankOperator's methods for consistency with corresponding operations of the full matrix
+template <class StaticMatrixType>
+class MatrixMultiplier
+{
+ public:
+ static StaticMatrixType matXmat(const StaticMatrixType& A, const StaticMatrixType& B)
+ {
+ DUNE_THROW(Dune::NotImplemented,"MatXmat not implemented for this MatrixType.");
+ }
+};
+
+template <class field_type, int N>
+class MatrixMultiplier<Dune::ScaledIdentityMatrix<field_type,N> >
+{
+ typedef Dune::ScaledIdentityMatrix<field_type,N> StaticMatrixType;
+ public:
+ static StaticMatrixType matXmat(const StaticMatrixType& A, const StaticMatrixType& B)
+ {
+ return StaticMatrixType(A.scalar()*B.scalar());
+ }
+ static StaticMatrixType matTXmat(const StaticMatrixType& A, const StaticMatrixType& B)
+ {
+ return matXmat(A,B);
+ }
+};
+
+template <class field_type, int N>
+class MatrixMultiplier<Dune::DiagonalMatrix<field_type,N> >
+{
+ typedef Dune::DiagonalMatrix<field_type,N> StaticMatrixType;
+ public:
+ static StaticMatrixType matXmat(const StaticMatrixType& A, const StaticMatrixType& B)
+ {
+ StaticMatrixType r(A);
+ for (size_t i=0; i<N; ++i)
+ r.diagonal(i) *= B.diagonal(i);
+ return r;
+ }
+ static StaticMatrixType matTXmat(const StaticMatrixType& A, const StaticMatrixType& B)
+ {
+ return matXmat(A,B);
+ }
+};
+
+
+template <class LowRankFactorType>
+bool check()
+{
+
+ bool passed = true;
+ typedef typename LowRankFactorType::block_type block_type;
+ size_t block_size = block_type::rows;
+
+ size_t m=3, n=1000;
+ typedef Dune::FieldVector<double, block_type::rows> Vblock_type;
+ typedef Dune::BlockVector<Vblock_type> Vector;
+
+ Vector x(n), b(n), ref_vec(n);
+ for (size_t i = 0; i<n; ++i)
+// x[i] = 1.0;
+ x[i] = (1.0*rand()) / RAND_MAX;
+ b = ref_vec = 0.0;
+
+ LowRankFactorType lr_factor(m,n);
+ std::cout << "Checking LowRankOperator<" << Dune::className(lr_factor) << " > ... ";
+
+ // set lowrankfactor to random values
+ typedef typename block_type::RowIterator BlockRowIterator;
+ typedef typename block_type::ColIterator BlockColIterator;
+ for (size_t row=0; row<m; ++row)
+ for (size_t col=0; col<n; ++col)
+ {
+ BlockRowIterator row_end = lr_factor[row][col].end();
+ for (BlockRowIterator row_it = lr_factor[row][col].begin(); row_it!=row_end; ++ row_it)
+ {
+ BlockColIterator col_end = row_it->end();
+ for (BlockColIterator col_it = row_it->begin(); col_it!=col_end; ++col_it)
+// *col_it = col+1;
+ *col_it = (1.0*rand())/RAND_MAX;
+ }
+ }
+
+ LowRankOperator<LowRankFactorType> lr_op(lr_factor);
+
+ LowRankFactorType full_op(n,n);
+
+ for (size_t i=0; i<n; ++i)
+ for (size_t j=0; j<n; ++j)
+ {
+ if (j>=i)
+ {
+ full_op[i][j] = 0.0;
+ for (size_t k=0; k<m; ++k)
+ {
+ full_op[i][j] += MatrixMultiplier<block_type>::matTXmat(lr_factor[k][i], lr_factor[k][j]);
+ }
+ }
+ else
+ full_op[i][j] = full_op[j][i];
+ }
+
+ full_op.mv(x,ref_vec);
+
+ lr_op.umv(x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if ((b[i] - ref_vec[i]).two_norm()/b[i].two_norm()>1e-12)
+ {
+ std::cout << "LowRankOperator::umv does not yield correct result. Difference = " << (b[i] - ref_vec[i]).two_norm() << std::endl;
+ passed = false;
+ for (size_t j=0; j<b.size(); ++j)
+ std::cout << b[j] << " "<<ref_vec[j] << std::endl;
+ std::cout << std::endl;
+ break;
+ }
+
+ b=1.0;
+ lr_op.mv(x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if ((b[i] - ref_vec[i]).two_norm()/b[i].two_norm()>1e-12)
+ {
+ std::cout << "LowRankOperator::mv does not yield correct result." << std::endl;
+ passed = false;
+ break;
+ }
+
+ LowRankFactorType lr_factor_reborn = lr_op.getLowRankFactor();
+
+ if (passed)
+ std::cout << "passed";
+ std::cout << std::endl;
+ return passed;
+}
+
+int main(int argc, char** argv) try
+{
+ static const int block_size = BLOCKSIZE;
+
+ bool passed(true);
+
+ passed = check<Dune::Matrix<Dune::ScaledIdentityMatrix<double, block_size> > >();
+ passed = passed and check<Dune::Matrix<Dune::DiagonalMatrix<double, block_size> > >();
+
+ return passed ? 0 : 1;
+}
+
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+
+ }
+
+
+
diff --git a/dune/solvers/test/nulloperatortest.cc b/dune/solvers/test/nulloperatortest.cc
new file mode 100644
index 00000000..c7f345cd
--- /dev/null
+++ b/dune/solvers/test/nulloperatortest.cc
@@ -0,0 +1,104 @@
+#include <config.h>
+
+#include <stdio.h>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/diagonalmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+#include <dune/solvers/operators/nulloperator.hh>
+
+#define BLOCKSIZE 2
+
+// this test tests the NullOperator class for several square block_types.
+// The test is mainly there to check if it compiles. It uses all implemented methods of NullOperator and
+// checks whether they come up with the expected behaviour.
+
+template <class Mblock_type>
+bool check()
+{
+ size_t block_size = Mblock_type::rows;
+ size_t n=1000;
+ typedef Dune::FieldVector<double, Mblock_type::rows> Vblock_type;
+ typedef Dune::BlockVector<Vblock_type> Vector;
+
+ Vector x(n), b(n), null_vec(n);
+ for (size_t i = 0; i<n; ++i)
+ x[i] = (1.0*rand()) / RAND_MAX;
+ b = null_vec = 0.0;
+
+ NullOperator<Mblock_type> null_op;
+
+ null_op.umv(x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if (b[i] != null_vec[i])
+ {
+ std::cout << "NullOperator::umv does not leave second argument unchanged." << std::endl;
+ return false;
+ }
+
+ null_op.usmv(1.0,x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if (b[i] != null_vec[i])
+ {
+ std::cout << "NullOperator::usmv does not leave second argument unchanged." << std::endl;
+ return false;
+ }
+
+ b=1.0;
+ null_op.mv(x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if (b[i] != null_vec[i])
+ {
+ std::cout << "NullOperator::mv does not set second argument to zero." << std::endl;
+ return false;
+ }
+
+ Mblock_type zero_block(0);
+
+ {
+ const Mblock_type block(null_op[n/2][n/4]);
+ if (block != zero_block)
+ {
+ std::cout << "NullOperator[][] does not yield correct zero block" << std::endl;
+ return false;
+ }
+ }
+
+ {
+ const Mblock_type block(null_op.diagonal(n/4));
+ if (block != zero_block)
+ {
+ std::cout << "NullOperator::diagonal does not yield correct zero block" << std::endl;
+ return false;
+ }
+ }
+
+ return true;
+}
+
+int main(int argc, char** argv) try
+{
+ static const int block_size = BLOCKSIZE;
+
+ bool passed;
+
+ passed = check<Dune::ScaledIdentityMatrix<double, block_size> >();
+ passed = passed and check<Dune::DiagonalMatrix<double, block_size> >();
+ passed = passed and check<Dune::FieldMatrix<double, block_size, block_size> >();
+
+ return passed ? 0 : 1;
+}
+
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+
+ }
+
+
diff --git a/dune/solvers/test/sumoperatortest.cc b/dune/solvers/test/sumoperatortest.cc
new file mode 100644
index 00000000..4ec4c848
--- /dev/null
+++ b/dune/solvers/test/sumoperatortest.cc
@@ -0,0 +1,157 @@
+#include <config.h>
+
+#include <stdio.h>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/istl/diagonalmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/operators/sumoperator.hh>
+#include <dune/solvers/operators/lowrankoperator.hh>
+
+#define BLOCKSIZE 2
+
+// This tests the SumOperator's methods for consistency with successive execution of corresponding operations
+template <class SparseMatrixType, class LowRankFactorType>
+bool check()
+{
+
+ bool passed = true;
+ typedef typename SparseMatrixType::block_type sp_block_type;
+ typedef typename LowRankFactorType::block_type lr_block_type;
+ size_t block_size = sp_block_type::rows;
+
+ assert(block_size==lr_block_type::rows);
+
+ size_t m=1, n=2;
+ typedef Dune::FieldVector<typename sp_block_type::field_type, sp_block_type::rows> Vblock_type;
+ typedef Dune::BlockVector<Vblock_type> Vector;
+
+ Vector x(n), b(n), ref_vec(n);
+ for (size_t i = 0; i<n; ++i)
+// x[i] = 1.0;
+ x[i] = (1.0*rand()) / RAND_MAX;
+ b = ref_vec = 0.0;
+
+ LowRankFactorType lr_factor(m,n);
+ SparseMatrixType sparse_matrix(n,n,4,SparseMatrixType::row_wise);
+
+ // create sparse Matrix as tridiagonal matrix
+ typedef typename SparseMatrixType::CreateIterator Iter;
+ Iter c_end = sparse_matrix.createend();
+ for(Iter row=sparse_matrix.createbegin(); row!=c_end; ++row)
+ {
+ // Add nonzeros for left neighbour, diagonal and right neighbour
+ if(row.index()>0)
+ row.insert(row.index()-1);
+ row.insert(row.index());
+ if(row.index()<sparse_matrix.N()-1)
+ row.insert(row.index()+1);
+ }
+
+ // Now the sparsity pattern is fully set up and we can add values
+ for (size_t i=0; i<n; ++i)
+ {
+ if (i>0)
+ sparse_matrix[i][i-1] = (1.0*rand())/RAND_MAX;
+ sparse_matrix[i][i]= (1.0*rand())/RAND_MAX;
+ if (i<n-1)
+ sparse_matrix[i][i+1] = (1.0*rand())/RAND_MAX;
+ }
+
+ // set lowrankfactor to random values
+ typedef typename lr_block_type::RowIterator BlockRowIterator;
+ typedef typename lr_block_type::ColIterator BlockColIterator;
+ for (size_t row=0; row<m; ++row)
+ for (size_t col=0; col<n; ++col)
+ {
+ BlockRowIterator row_end = lr_factor[row][col].end();
+ for (BlockRowIterator row_it = lr_factor[row][col].begin(); row_it!=row_end; ++ row_it)
+ {
+ BlockColIterator col_end = row_it->end();
+ for (BlockColIterator col_it = row_it->begin(); col_it!=col_end; ++col_it)
+// *col_it = col+1;
+ *col_it = (1.0*rand())/RAND_MAX;
+ }
+ }
+
+ LowRankOperator<LowRankFactorType> lr_op(lr_factor);
+
+ std::cout << "Checking SumOperator<" << Dune::className(sparse_matrix) << ", " << Dune::className(lr_op) << " > ... ";
+
+ SumOperator<SparseMatrixType,LowRankOperator<LowRankFactorType> > sum_op(sparse_matrix, lr_op);
+
+
+ lr_op.umv(x,ref_vec);
+ sparse_matrix.umv(x,ref_vec);
+
+ sum_op.umv(x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if ((b[i] - ref_vec[i]).two_norm()>1e-15)
+ {
+ std::cout << "SumOperator::umv does not yield correct result." << std::endl;
+ passed = false;
+ for (size_t j=0; j<b.size(); ++j)
+ std::cout << b[j] << " "<<ref_vec[j] << std::endl;
+ std::cout << std::endl;
+ break;
+ }
+
+ b=1.0;
+
+ sum_op.mv(x,b);
+ for (size_t i=0; i<b.size(); ++i)
+ if ((b[i] - ref_vec[i]).two_norm()>1e-15)
+ {
+ std::cout << "SumOperator::mv does not yield correct result." << std::endl;
+ passed = false;
+ break;
+ }
+
+ LowRankOperator<LowRankFactorType> lr_factor_reborn = sum_op.getLowRankMatrix();
+ SparseMatrixType sparse_matrix_reborn = sum_op.getSparseMatrix();
+
+ if (passed)
+ std::cout << "passed";
+ std::cout << std::endl;
+ return passed;
+}
+
+int main(int argc, char** argv) try
+{
+ static const int block_size = BLOCKSIZE;
+
+ bool passed(true);
+
+ typedef Dune::ScaledIdentityMatrix<double, block_size> Block1;
+ typedef Dune::DiagonalMatrix<double, block_size> Block2;
+ typedef Dune::FieldMatrix<double, block_size, block_size> Block3;
+
+ passed = check<Dune::BCRSMatrix<Block1>, Dune::Matrix<Block1> >();
+ passed = passed and check<Dune::BCRSMatrix<Block1>, Dune::Matrix<Block2> >();
+ passed = passed and check<Dune::BCRSMatrix<Block1>, Dune::Matrix<Block3> >();
+
+ passed = passed and check<Dune::BCRSMatrix<Block2>, Dune::Matrix<Block1> >();
+ passed = passed and check<Dune::BCRSMatrix<Block2>, Dune::Matrix<Block2> >();
+ passed = passed and check<Dune::BCRSMatrix<Block2>, Dune::Matrix<Block3> >();
+
+ passed = passed and check<Dune::BCRSMatrix<Block3>, Dune::Matrix<Block1> >();
+ passed = passed and check<Dune::BCRSMatrix<Block3>, Dune::Matrix<Block2> >();
+ passed = passed and check<Dune::BCRSMatrix<Block3>, Dune::Matrix<Block3> >();
+
+ return passed ? 0 : 1;
+}
+
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+
+ }
+
--
GitLab