Skip to content
Snippets Groups Projects
Commit 96542556 authored by Oliver Sander's avatar Oliver Sander
Browse files

Move code into the generic template class SuperLU<M>

Previously, all code was in the specialization of that class
for BCRSMatrix<FieldMatrix>.  That was overly restrictive:
most of the code does not actually rely on the precise type
of the matrix.  As a first step towards allowing more general
matrix type we move all the code into the general class.
parent 883d8a2f
Branches
Tags
1 merge request!270Implement superlu for scalar valued matrices
......@@ -31,10 +31,6 @@ namespace Dune
* @author Markus Blatt
* @brief Classes for using SuperLU with ISTL matrices.
*/
template<class Matrix>
class SuperLU
{};
template<class M, class T, class TM, class TD, class TA>
class SeqOverlappingSchwarz;
......@@ -246,6 +242,26 @@ namespace Dune
};
#endif
namespace Impl
{
template<class M>
struct SuperLUVectorChooser
{};
template<typename T, typename A, int n, int m>
struct SuperLUVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
{
/** @brief The type of the domain of the solver */
using domain_type = BlockVector<
FieldVector<T,m>,
typename A::template rebind<FieldVector<T,m> >::other>;
/** @brief The type of the range of the solver */
using range_type = BlockVector<
FieldVector<T,n>,
typename A::template rebind<FieldVector<T,n> >::other>;
};
}
/**
* @brief SuperLu Solver
*
......@@ -259,30 +275,25 @@ namespace Dune
* 2: std::complex<float>, 3: std::complex<double>)
* if the numeric type should be different from double.
*/
template<typename T, typename A, int n, int m>
class SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A > >
template<typename M>
class SuperLU
: public InverseOperator<
BlockVector<FieldVector<T,m>,
typename A::template rebind<FieldVector<T,m> >::other>,
BlockVector<FieldVector<T,n>,
typename A::template rebind<FieldVector<T,n> >::other> >
typename Impl::SuperLUVectorChooser<M>::domain_type,
typename Impl::SuperLUVectorChooser<M>::range_type >
{
using T = typename M::field_type;
public:
/** @brief The matrix type. */
typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
using Matrix = M;
using matrix_type = M;
/** @brief The corresponding SuperLU Matrix type.*/
typedef Dune::SuperLUMatrix<Matrix> SuperLUMatrix;
/** @brief Type of an associated initializer class. */
typedef SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> > MatrixInitializer;
typedef SuperMatrixInitializer<Matrix> MatrixInitializer;
/** @brief The type of the domain of the solver. */
typedef Dune::BlockVector<
FieldVector<T,m>,
typename A::template rebind<FieldVector<T,m> >::other> domain_type;
using domain_type = typename Impl::SuperLUVectorChooser<M>::domain_type;
/** @brief The type of the range of the solver. */
typedef Dune::BlockVector<
FieldVector<T,n>,
typename A::template rebind<FieldVector<T,n> >::other> range_type;
using range_type = typename Impl::SuperLUVectorChooser<M>::range_type;
//! Category of the solver (see SolverCategory::Category)
virtual SolverCategory::Category category() const
......@@ -356,7 +367,7 @@ namespace Dune
const char* name() { return "SuperLU"; }
private:
template<class M,class X, class TM, class TD, class T1>
template<class Mat,class X, class TM, class TD, class T1>
friend class SeqOverlappingSchwarz;
friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
......@@ -377,16 +388,16 @@ namespace Dune
bool first, verbose, reusevector;
};
template<typename T, typename A, int n, int m>
SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
template<typename M>
SuperLU<M>
::~SuperLU()
{
if(mat.N()+mat.M()>0)
free();
}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
template<typename M>
void SuperLU<M>::free()
{
delete[] perm_c;
delete[] perm_r;
......@@ -405,8 +416,8 @@ namespace Dune
mat.free();
}
template<typename T, typename A, int n, int m>
SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
template<typename M>
SuperLU<M>
::SuperLU(const Matrix& mat_, bool verbose_, bool reusevector_)
: work(0), lwork(0), first(true), verbose(verbose_),
reusevector(reusevector_)
......@@ -414,19 +425,19 @@ namespace Dune
setMatrix(mat_);
}
template<typename T, typename A, int n, int m>
SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
template<typename M>
SuperLU<M>::SuperLU()
: work(0), lwork(0),verbose(false),
reusevector(false)
{}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(bool v)
template<typename M>
void SuperLU<M>::setVerbosity(bool v)
{
verbose=v;
}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(const Matrix& mat_)
template<typename M>
void SuperLU<M>::setMatrix(const Matrix& mat_)
{
if(mat.N()+mat.M()>0) {
free();
......@@ -438,9 +449,9 @@ namespace Dune
decompose();
}
template<typename T, typename A, int n, int m>
template<typename M>
template<class S>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(const Matrix& mat_,
void SuperLU<M>::setSubMatrix(const Matrix& mat_,
const S& mrs)
{
if(mat.N()+mat.M()>0) {
......@@ -453,8 +464,8 @@ namespace Dune
decompose();
}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
template<typename M>
void SuperLU<M>::decompose()
{
first = true;
......@@ -544,8 +555,8 @@ namespace Dune
options.Fact = FACTORED;
}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
template<typename M>
void SuperLU<M>
::apply(domain_type& x, range_type& b, InverseOperatorResult& res)
{
if (mat.N() != b.dim())
......@@ -637,8 +648,8 @@ namespace Dune
}
}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
template<typename M>
void SuperLU<M>
::apply(T* x, T* b)
{
if(mat.N()+mat.M()==0)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment