Skip to content
Snippets Groups Projects
Commit bce95bfd authored by Uli Sack's avatar Uli Sack Committed by usack
Browse files

updated sumoperatortest to latest changes in SumOperator

[[Imported from SVN: r5484]]
parent 370d0b0d
Branches
Tags
No related merge requests found
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
#include <dune/istl/bvector.hh> #include <dune/istl/bvector.hh>
#include <dune/istl/matrix.hh> #include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/diagonalmatrix.hh> #include <dune/istl/diagonalmatrix.hh>
#include <dune/istl/scaledidmatrix.hh> #include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/io.hh> #include <dune/istl/io.hh>
...@@ -24,24 +25,28 @@ bool check() ...@@ -24,24 +25,28 @@ bool check()
{ {
bool passed = true; bool passed = true;
typedef LowRankOperator<LowRankFactorType> LowRankMatrixType;
typedef typename SparseMatrixType::block_type sp_block_type; typedef typename SparseMatrixType::block_type sp_block_type;
typedef typename LowRankFactorType::block_type lr_block_type; typedef typename LowRankFactorType::block_type lr_block_type;
size_t block_size = sp_block_type::rows; size_t block_size = sp_block_type::rows;
assert(block_size==lr_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::FieldVector<typename sp_block_type::field_type, sp_block_type::rows> Vblock_type;
typedef Dune::BlockVector<Vblock_type> Vector; 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) for (size_t i = 0; i<n; ++i)
// x[i] = 1.0; // x[i] = 1.0;
x[i] = (1.0*rand()) / RAND_MAX; 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); 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 // create sparse Matrix as tridiagonal matrix
typedef typename SparseMatrixType::CreateIterator Iter; typedef typename SparseMatrixType::CreateIterator Iter;
...@@ -57,15 +62,25 @@ bool check() ...@@ -57,15 +62,25 @@ bool check()
} }
// Now the sparsity pattern is fully set up and we can add values // 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) for (size_t i=0; i<n; ++i)
{ {
if (i>0) 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) 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 // set lowrankfactor to random values
typedef typename lr_block_type::RowIterator BlockRowIterator; typedef typename lr_block_type::RowIterator BlockRowIterator;
typedef typename lr_block_type::ColIterator BlockColIterator; typedef typename lr_block_type::ColIterator BlockColIterator;
...@@ -82,7 +97,9 @@ bool check() ...@@ -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) << " > ... "; std::cout << "Checking SumOperator<" << Dune::className(sparse_matrix) << ", " << Dune::className(lr_op) << " > ... ";
...@@ -93,8 +110,9 @@ bool check() ...@@ -93,8 +110,9 @@ bool check()
sparse_matrix.umv(x,ref_vec); sparse_matrix.umv(x,ref_vec);
sum_op.umv(x,b); sum_op.umv(x,b);
sum_op_default.umv(x,c);
for (size_t i=0; i<b.size(); ++i) 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; std::cout << "SumOperator::umv does not yield correct result." << std::endl;
passed = false; passed = false;
...@@ -104,19 +122,31 @@ bool check() ...@@ -104,19 +122,31 @@ bool check()
break; break;
} }
b=1.0; c=b=1.0;
sum_op.mv(x,b); sum_op.mv(x,b);
sum_op_default.mv(x,c);
for (size_t i=0; i<b.size(); ++i) 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; std::cout << "SumOperator::mv does not yield correct result." << std::endl;
passed = false; passed = false;
break; break;
} }
LowRankOperator<LowRankFactorType> lr_factor_reborn = sum_op.getLowRankMatrix(); sum_op.mmv(x,b);
SparseMatrixType sparse_matrix_reborn = sum_op.getSparseMatrix(); 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) if (passed)
std::cout << "passed"; std::cout << "passed";
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment