diff --git a/dune/istl/diagonalmatrix.hh b/dune/istl/diagonalmatrix.hh index 332d3b6e37c4a697c3252d49057eb4437159bd17..128a0e391485915370b5dd8eba110bd7c6f0450b 100644 --- a/dune/istl/diagonalmatrix.hh +++ b/dune/istl/diagonalmatrix.hh @@ -1060,8 +1060,8 @@ namespace Dune { - template<class K, int n> - void istl_assign_to_fmatrix(FieldMatrix<K,n,n>& fm, const DiagonalMatrix<K,n>& s) + 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) diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh index f57e16aac86ab3f6a51a87e1749b851a9983b3bc..425211b7329210f36f8a0ff6d66f1e520ea9962d 100644 --- a/dune/istl/paamg/amg.hh +++ b/dune/istl/paamg/amg.hh @@ -388,6 +388,35 @@ namespace Dune template<class M, class X, class S, class PI, class A> void AMG<M,X,S,PI,A>::pre(Domain& x, Range& b) { + // Detect Matrix rows where all offdiagonal entries are + // zero and set x such that A_dd*x_d=b_d + // Thus users can be more careless when setting up their linear + // systems. + typedef typename M::matrix_type Matrix; + typedef typename Matrix::ConstRowIterator RowIter; + 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; + for(ColIter col=row->begin(); col!=row->end(); ++col) { + if(row.index()==col.index()) { + diagonal = *col; + hasDiagonal = false; + }else{ + if(*col!=zero) + isDirichlet = false; + } + } + if(isDirichlet) + diagonal.solve(x[row.index()], b[row.index()]); + } + if(smoothers_.levels()>0) smoothers_.finest()->pre(x,b); else diff --git a/dune/istl/test/.gitignore b/dune/istl/test/.gitignore index a0787d050b89d48b2aee8c920d3c0dc199da1f1d..3614259ef262aeb760e2f7382af3b96af0a47d7f 100644 --- a/dune/istl/test/.gitignore +++ b/dune/istl/test/.gitignore @@ -4,6 +4,7 @@ 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 752ab52f6680c3791b46849c7c8ed73c1ed900cf..cc3d48b4ce79ef8f13eca80a7a17af1db7394404 100644 --- a/dune/istl/test/Makefile.am +++ b/dune/istl/test/Makefile.am @@ -17,8 +17,21 @@ if PARDISO endif # which tests where program to build and run are equal -NORMALTESTS = basearraytest matrixutilstest matrixtest mmtest bvectortest complexrhstest dotproducttest vbvectortest \ - bcrsbuildtest matrixiteratortest mv iotest scaledidmatrixtest seqmatrixmarkettest +NORMALTESTS = basearraytest \ + bcrsbuildtest \ + bvectortest \ + complexrhstest \ + diagonalmatrixtest \ + dotproducttest \ + iotest \ + matrixiteratortest \ + matrixtest \ + matrixutilstest \ + mmtest \ + mv \ + scaledidmatrixtest \ + seqmatrixmarkettest \ + vbvectortest # list of tests to run (indicestest is special case) TESTS = $(NORMALTESTS) $(MPITESTS) $(SUPERLUTESTS) $(PARDISOTEST) $(PARMETISTESTS) @@ -75,6 +88,8 @@ 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 new file mode 100644 index 0000000000000000000000000000000000000000..92247b28f899a7908835524c53eadf99adc184f3 --- /dev/null +++ b/dune/istl/test/diagonalmatrixtest.cc @@ -0,0 +1,64 @@ +// -*- 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/scaledidmatrixtest.cc b/dune/istl/test/scaledidmatrixtest.cc index c0f7403ae3eaed207e55153c70190e4cb400f3ae..dad86d8bba01a8edd3e9af624db68d976aa7a310 100644 --- a/dune/istl/test/scaledidmatrixtest.cc +++ b/dune/istl/test/scaledidmatrixtest.cc @@ -46,6 +46,9 @@ void test_matrix() std::cout << A << std::endl; // print vector std::cout << f << std::endl; + + // Construction of FieldMatrix from ScaledIdentityMatrix + FieldMatrix<K,n,n> AFM DUNE_UNUSED = FieldMatrix<K,n,n>(A); } int main()