Skip to content
Snippets Groups Projects
Commit ced5acb6 authored by Markus Blatt's avatar Markus Blatt
Browse files

Merged latest changes from the trunk

[[Imported from SVN: r1710]]
parents d894a599 432e2a6e
No related branches found
No related tags found
No related merge requests found
......@@ -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)
......
......@@ -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
......
......@@ -4,6 +4,7 @@ Makefile.in
.libs
semantic.cache
complexrhstest
diagonalmatrixtest
dotproducttest
bvectortest
matrixutilstest
......
......@@ -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
......
// -*- 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;
}
}
......@@ -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()
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment