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

now, here are the classes and tests

[[Imported from SVN: r5238]]
parent a3946d5c
Branches
No related tags found
No related merge requests found
SUBDIRS =
operatorsdir = $(includedir)/dune/solvers/operators
operators_HEADERS = lowrankoperator.hh \
nulloperator.hh \
sumoperator.hh
include $(top_srcdir)/am/global-rules
#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
#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
#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
# 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
#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;
}
#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;
}
#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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment