diff --git a/dune/solvers/test/sumoperatortest.cc b/dune/solvers/test/sumoperatortest.cc index 4ec4c848506929765c3aa74088f0cfe403395269..ffd2e21e43699c3127c17cd249687dfa163802fb 100644 --- a/dune/solvers/test/sumoperatortest.cc +++ b/dune/solvers/test/sumoperatortest.cc @@ -9,6 +9,7 @@ #include <dune/istl/bvector.hh> #include <dune/istl/matrix.hh> +#include <dune/istl/matrixindexset.hh> #include <dune/istl/diagonalmatrix.hh> #include <dune/istl/scaledidmatrix.hh> #include <dune/istl/io.hh> @@ -24,24 +25,28 @@ bool check() { bool passed = true; + typedef LowRankOperator<LowRankFactorType> LowRankMatrixType; 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; + size_t m=3, n=10; 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); + Vector x(n), b(n), c(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; + b = c = ref_vec = 0.0; + + SumOperator<SparseMatrixType,LowRankMatrixType> sum_op_default; LowRankFactorType lr_factor(m,n); - SparseMatrixType sparse_matrix(n,n,4,SparseMatrixType::row_wise); + SparseMatrixType sparse_matrix(n,n,3*n,SparseMatrixType::row_wise); + // create sparse Matrix as tridiagonal matrix typedef typename SparseMatrixType::CreateIterator Iter; @@ -57,15 +62,25 @@ bool check() } // Now the sparsity pattern is fully set up and we can add values + double entry=0.0; 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; + { + entry = (1.0*rand())/RAND_MAX; + sparse_matrix[i][i-1] = entry; + } + entry = (1.0*rand())/RAND_MAX; + sparse_matrix[i][i]= entry; if (i<n-1) - sparse_matrix[i][i+1] = (1.0*rand())/RAND_MAX; + { + entry = (1.0*rand())/RAND_MAX; + sparse_matrix[i][i+1] = entry; + } } + sum_op_default.sparseMatrix() = sparse_matrix; + // set lowrankfactor to random values typedef typename lr_block_type::RowIterator BlockRowIterator; typedef typename lr_block_type::ColIterator BlockColIterator; @@ -82,7 +97,9 @@ bool check() } } - LowRankOperator<LowRankFactorType> lr_op(lr_factor); + sum_op_default.lowRankMatrix().lowRankFactor() = lr_factor; + + LowRankMatrixType lr_op(lr_factor); std::cout << "Checking SumOperator<" << Dune::className(sparse_matrix) << ", " << Dune::className(lr_op) << " > ... "; @@ -93,8 +110,9 @@ bool check() sparse_matrix.umv(x,ref_vec); sum_op.umv(x,b); + sum_op_default.umv(x,c); for (size_t i=0; i<b.size(); ++i) - if ((b[i] - ref_vec[i]).two_norm()>1e-15) + if ((b[i] - ref_vec[i]).two_norm()>1e-12 or (c[i] - ref_vec[i]).two_norm()>1e-12) { std::cout << "SumOperator::umv does not yield correct result." << std::endl; passed = false; @@ -104,19 +122,31 @@ bool check() break; } - b=1.0; + c=b=1.0; sum_op.mv(x,b); + sum_op_default.mv(x,c); for (size_t i=0; i<b.size(); ++i) - if ((b[i] - ref_vec[i]).two_norm()>1e-15) + if ((b[i] - ref_vec[i]).two_norm()>1e-12 or (c[i] - ref_vec[i]).two_norm()>1e-12) { 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(); + sum_op.mmv(x,b); + sum_op_default.mmv(x,c); + for (size_t i=0; i<b.size(); ++i) + if (b[i].two_norm()>1e-12 or c[i].two_norm()>1e-12) + { + std::cout << "SumOperator::mmv does not yield correct result." << std::endl; + passed = false; + break; + } + + + LowRankOperator<LowRankFactorType> lr_factor_reborn = sum_op.lowRankMatrix(); + SparseMatrixType sparse_matrix_reborn = sum_op.sparseMatrix(); if (passed) std::cout << "passed";