diff --git a/AMDiS/src/AMDiS.cc b/AMDiS/src/AMDiS.cc
index 9ea6d73f083da4d98de20226545a14811d5dbb70..a5b3e4440152387dfe84bfe5566c53df33e3c90f 100644
--- a/AMDiS/src/AMDiS.cc
+++ b/AMDiS/src/AMDiS.cc
@@ -29,6 +29,10 @@
#endif
#include "boost/program_options.hpp"
+#if defined HAVE_PETSC || defined HAVE_SEQ_PETSC || defined HAVE_PARALLEL_PETSC
+#include <petsc.h>
+#endif
+
namespace AMDiS {
using namespace std;
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index f27c34ab5bf960bbc7790114e6c34920cc51d161..762b35de9cd50b1b3113a15e727f1942717c06d5 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -154,9 +154,7 @@
#include "parallel/StdMpi.h"
#include "parallel/ParallelProblemStat.h"
-#if HAVE_PARALLEL_MTL4
-// #include "parallel/PMTL4Solver.h"
-#else
+#if !HAVE_PARALLEL_MTL4
#include "parallel/PetscSolver.h"
#include "parallel/PetscSolverNavierStokes.h"
#endif
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index 8ea6e49f88d3556b3926bc6af5576ff1a31f8996..229f7616bb889209759a59f03cde371a1cc5eb6b 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -35,6 +35,7 @@ namespace AMDiS {
class AdaptStationary;
class Assembler;
class BasisFunction;
+ struct BlockMTLMatrix;
class BoundaryManager;
// class CGSolver;
class CoarseningManager;
diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc
index 36c001a8c01c40d1a375f2830be0ff5748df194f..5a88adb91654713a267e15eea6676c9a66d2cce0 100644
--- a/AMDiS/src/CreatorMap.cc
+++ b/AMDiS/src/CreatorMap.cc
@@ -25,6 +25,7 @@
#include "MTL4Types.h"
#include "solver/LinearSolver.h"
#include "solver/ITL_Solver.h"
+#include "solver/BITL_Solver.h"
#include "solver/ITL_Preconditioner.h"
#include "solver/UmfPackSolver.h"
#include "solver/KrylovPreconditioner.h"
@@ -61,6 +62,9 @@ namespace AMDiS {
{
LinearSolverCreator *creator;
+ // creators for MTL4 krylov solvers
+ // _________________________________________________________________________
+
creator = new CGSolver::Creator;
addCreator("mtl_cg", creator);
@@ -115,11 +119,64 @@ namespace AMDiS {
addCreator("mtl_hypre", creator);
#endif
+ // creators for block krylov solvers
+ // _________________________________________________________________________
+
+ creator = new B_CGSolver::Creator;
+ addCreator("bmtl_cg", creator);
+
+ creator = new B_CGSSolver::Creator;
+ addCreator("bmtl_cgs", creator);
+
+// creator = new B_BiCGSolver::Creator;
+// addCreator("mtl_bicg", creator);
+
+ creator = new B_BiCGStabSolver::Creator;
+ addCreator("bmtl_bicgstab", creator);
+
+// creator = new B_BiCGStab2Solver::Creator;
+// addCreator("bmtl_bicgstab2", creator); // MTL Error
+
+// creator = new B_BiCGStabEllSolver::Creator;
+// addCreator("bmtl_bicgstab_ell", creator); // MTL Error
+
+// creator = new QMRSolver::Creator;
+// addCreator("mtl_qmr", creator);
+
+ creator = new B_TFQMRSolver::Creator;
+ addCreator("bmtl_tfqmr", creator);
+
+ creator = new B_GMResSolver::Creator;
+ addCreator("bmtl_gmres", creator);
+
+ creator = new B_GcrSolver::Creator;
+ addCreator("bmtl_gcr", creator);
+
+ creator = new B_FGMResSolver::Creator;
+ addCreator("bmtl_fgmres", creator);
+
+// creator = new B_IDRsSolver::Creator;
+// addCreator("bmtl_idr_s", creator); // MTL Error
+
+ creator = new B_MinResSolver::Creator;
+ addCreator("bmtl_minres", creator);
+
+ creator = new B_PreOnly::Creator;
+ addCreator("bmtl_preonly", creator);
+ addCreator("bmtl_richardson", creator);
+
+ // creators for PETSc solvers
+ // _________________________________________________________________________
+
#if defined HAVE_SEQ_PETSC
// sequential PETSc-Solver
- creator = new PetscSolver<>::Creator;
+ creator = new PetscSolver::Creator;
addCreator("petsc_petsc", creator); // standard creator for petsc solver
addCreator("petsc", creator);
+
+ LinearSolverCreator* creator2 = new PetscSolverNested::Creator;
+ addCreator("bpetsc_petsc", creator2); // standard creator for petsc solver
+ addCreator("bpetsc", creator2);
std::map<std::string,std::string>::iterator it;
PetscParameters params;
@@ -127,7 +184,8 @@ namespace AMDiS {
it!= params.solverMap.end();
it++) {
CreatorMap< LinearSolver >::addCreator("petsc_" + it->first, creator);
- }
+ CreatorMap< LinearSolver >::addCreator("bpetsc_" + it->first, creator2);
+ }
#endif
}
@@ -159,6 +217,22 @@ namespace AMDiS {
}
+ template<>
+ void CreatorMap<ITL_BasePreconditioner<BlockMTLMatrix, MTLTypes::MTLVector> >::addDefaultCreators()
+ {
+ addCreator("no", new BlockIdentityPreconditioner::Creator);
+ addCreator("diag", new BlockDiagonalPreconditioner::Creator);
+ }
+
+
+#if defined HAVE_SEQ_PETSC
+ template<>
+ void CreatorMap<PetscPreconditioner>::addDefaultCreators() { }
+
+ template<>
+ void CreatorMap<PetscPreconditionerNested>::addDefaultCreators() { }
+#endif
+
template<>
void CreatorMap<Estimator>::addDefaultCreators()
{
diff --git a/AMDiS/src/CreatorMap.h b/AMDiS/src/CreatorMap.h
index 212fd6727c0ae0fabaae392cfbd1942bf62872cb..8da0c1fe7d647a3dcfc755a8454485bbeb946419 100644
--- a/AMDiS/src/CreatorMap.h
+++ b/AMDiS/src/CreatorMap.h
@@ -43,13 +43,16 @@ namespace AMDiS {
template<typename BaseClass>
class CreatorMap
{
+ public:
+ typedef std::map< std::string, CreatorInterface<BaseClass>* > CreatorMapType;
+
public:
/// Adds a new creator together with the given key to the map.
static void addCreator(std::string key, CreatorInterface<BaseClass>* creator)
{
FUNCNAME("CreatorMap::addCreator()");
init();
- TEST_EXIT(creatorMap[key] == nullptr)
+ TEST_EXIT(creatorMap[key] == NULL)
("there is already a creator for key %s\n",key.c_str());
creatorMap[key] = creator;
}
@@ -66,8 +69,10 @@ namespace AMDiS {
FUNCNAME("CreatorMap::getCreator()");
init();
CreatorInterface<BaseClass> *creator = creatorMap[key];
- TEST_EXIT(creator)("No creator for key \"%s\" defined in init file for parameter \"%s\"\n",
- key.c_str(), initFileStr.c_str());
+ TEST_EXIT(creator)
+ ("No creator for key \"%s\" defined in init file for parameter \"%s\"\n",
+ key.c_str(), initFileStr.c_str());
+
return creator;
}
@@ -89,13 +94,13 @@ namespace AMDiS {
protected:
/// STL map containing the pairs of keys and creators.
- static std::map< std::string, CreatorInterface<BaseClass>* > creatorMap;
+ static CreatorMapType creatorMap;
static bool initialized;
};
template<typename BaseClass>
- std::map< std::string, CreatorInterface<BaseClass>* > CreatorMap<BaseClass>::creatorMap;
+ typename CreatorMap<BaseClass>::CreatorMapType CreatorMap<BaseClass>::creatorMap;
template<typename BaseClass>
bool CreatorMap<BaseClass>::initialized = false;
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index ea2aae87bc72470d2edca30644acfcf51f542380..a4960615119cc1b7ee0e939075a9c6bc56e9c55f 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -67,10 +67,6 @@ namespace AMDiS {
if (!colFeSpace)
colFeSpace = rowFeSpace;
-#if 0
- if (rowFeSpace && rowFeSpace->getAdmin())
- (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->addDOFIndexed(this);
-#endif
boundaryManager = new BoundaryManager(rowFeSpace);
nRow = rowFeSpace->getBasisFcts()->getNumber();
@@ -87,10 +83,6 @@ namespace AMDiS {
FUNCNAME("DOFMatrix::DOFMatrix()");
*this = rhs;
-#if 0
- if (rowFeSpace && rowFeSpace->getAdmin())
- (const_cast<DOFAdmin*>( rowFeSpace->getAdmin()))->addDOFIndexed(this);
-#endif
TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
inserter = 0;
}
@@ -98,10 +90,6 @@ namespace AMDiS {
DOFMatrix::~DOFMatrix()
{
-#if 0
- if (rowFeSpace && rowFeSpace->getAdmin())
- (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->removeDOFIndexed(this);
-#endif
if (boundaryManager)
delete boundaryManager;
if (inserter)
@@ -116,43 +104,6 @@ namespace AMDiS {
}
- bool DOFMatrix::symmetric()
- {
- double tol = 1e-5;
-
- using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
- namespace traits= mtl::traits;
- typedef base_matrix_type Matrix;
-
- traits::row<Matrix>::type row(matrix);
- traits::col<Matrix>::type col(matrix);
- traits::const_value<Matrix>::type value(matrix);
-
- typedef traits::range_generator<major, Matrix>::type cursor_type;
- typedef traits::range_generator<nz, cursor_type>::type icursor_type;
-
- for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
- for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
- // Compare each non-zero entry with its transposed
- if (abs(value(*icursor) - matrix[col(*icursor)][row(*icursor)]) > tol)
- return false;
- return true;
- }
-
-
- void DOFMatrix::test()
- {
- FUNCNAME("DOFMatrix::test()");
-
- int non_symmetric = !symmetric();
-
- if (non_symmetric)
- MSG("Matrix `%s' not symmetric.\n", name.data());
- else
- MSG("Matrix `%s' is symmetric.\n", name.data());
- }
-
-
DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
{
rowFeSpace = rhs.rowFeSpace;
@@ -290,14 +241,11 @@ namespace AMDiS {
}
+#if 0
double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
{
return matrix[a][b];
}
-
-#if 0
- void DOFMatrix::freeDOFContent(int index)
- {}
#endif
void DOFMatrix::assemble(double factor,
@@ -418,7 +366,7 @@ namespace AMDiS {
void DOFMatrix::finishAssembling()
{
- // call the operatos cleanup procedures
+ // call the operators cleanup procedures
for (std::vector<Operator*>::iterator it = operators.begin();
it != operators.end(); ++it)
(*it)->finishAssembling();
@@ -436,7 +384,7 @@ namespace AMDiS {
return fillFlag;
}
-
+#if 0
void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
{
matrix += a * x.matrix + y.matrix;
@@ -447,7 +395,7 @@ namespace AMDiS {
{
matrix *= b;
}
-
+#endif
void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor)
{
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 7f44bb925cdfa296ff5725a9ed3e5402f19aeb46..f104cf23ff326579dfbb95eba7f682100f51a2ff 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -69,6 +69,9 @@ namespace AMDiS {
/// Symbolic constant for an unused entry which is not followed by any
/// used entry in this row
static const int NO_MORE_ENTRIES = -2;
+
+ /// Minimal slot-size used for the inserter. See \ref startInsertion for details.
+ static const int MIN_NNZ_PER_ROW = 20;
public:
DOFMatrix();
@@ -100,56 +103,7 @@ namespace AMDiS {
{
return matrix;
}
-#if 0
- // Only to get rid of the abstract functions, I hope they are not used
- std::vector<bool>::iterator begin()
- {
- ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v;
- return v.begin();
- }
-
- std::vector<bool>::iterator end()
- {
- ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v;
- return v.end();
- }
-
- // Only to get rid of the abstract functions, I hope they are not used
- std::vector<bool>::const_iterator begin() const
- {
- ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v;
- return v.begin();
- }
-
- std::vector<bool>::const_iterator end() const
- {
- ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v;
- return v.end();
- }
-
- bool dummy; // Must be deleted later
- reference operator[](DegreeOfFreedom i)
- {
- ERROR_EXIT("Shouldn't be used, only fake.");
- std::vector<bool> d; d.push_back(dummy);
- return d[0];
- }
-
- const_reference operator[](DegreeOfFreedom i) const
- {
- ERROR_EXIT("Shouldn't be used, only fake.");
- std::vector<bool> d; d.push_back(dummy);
- return d[0];
- }
-
- /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
- void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF)
- {}
-
- /// Implements DOFIndexedBase::freeDOFContent()
- virtual void freeDOFContent(int index);
-#endif
/// Returns \ref coupleMatrix.
inline bool isCoupleMatrix()
{
@@ -162,11 +116,13 @@ namespace AMDiS {
coupleMatrix = c;
}
+#if 0
/// a * x + y
void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
/// Multiplication with a scalar.
void scal(double s);
+#endif
/// Adds an operator to the DOFMatrix. A factor, that is multipled to the
/// operator, and a multilier factor for the estimator may be also given.
@@ -267,7 +223,7 @@ namespace AMDiS {
/// entries per row (per column for CCS matrices). Choosing this parameter
/// too small can induce perceivable overhead for compressed matrices. Thus,
/// it's better to choose a bit too large than too small.
- void startInsertion(int nnz_per_row = 10);
+ void startInsertion(int nnz_per_row = MIN_NNZ_PER_ROW);
/// Finishes insertion. For compressed matrix types, this is where the
/// compression happens.
@@ -281,15 +237,6 @@ namespace AMDiS {
inserter= 0;
}
-#if 0
- /// Returns whether restriction should be performed after coarsening
- /// (false by default)
- virtual bool coarseRestrict()
- {
- return false;
- }
-#endif
-
/// Returns const \ref rowFeSpace
const FiniteElemSpace* getRowFeSpace() const
{
@@ -332,12 +279,6 @@ namespace AMDiS {
return name;
}
#if 0
- /// Resizes \ref matrix to n rows
- inline void resize(int n)
- {
- TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
- }
-#endif
/// Returns value at logical indices a,b
double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
@@ -349,6 +290,7 @@ namespace AMDiS {
void addSparseDOFEntry(double sign,
int irow, int jcol, double entry,
bool add = true);
+#endif
void clearDirichletRows();
@@ -361,11 +303,6 @@ namespace AMDiS {
set_to_zero(matrix);
}
- /// Test whether \ref matrix is symmetric. Exits if not.
- void test();
-
- bool symmetric();
-
inline std::vector<Operator*>& getOperators()
{
return operators;
@@ -398,8 +335,8 @@ namespace AMDiS {
if (num_rows(matrix) != 0)
nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2);
- if (nnzPerRow < 20)
- nnzPerRow = 20;
+ if (nnzPerRow < MIN_NNZ_PER_ROW)
+ nnzPerRow = MIN_NNZ_PER_ROW;
}
/// Returns \ref nnzPerRow.
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 7d46bc825b076db07bbb67a67ee70f8631261eee..f396a8144e3d09dab7d29b46c42dd963d562fd39 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -110,7 +110,7 @@ namespace AMDiS {
delete oldElInfo;
} else
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
-
+
if (oldElInfo)
oldElInfo = elInfo;
@@ -122,7 +122,7 @@ namespace AMDiS {
value = basFcts->evalUh(lambda, uh);
} else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- return 0.0;
+ value = 0.0;
#else
ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
#endif
@@ -132,6 +132,9 @@ namespace AMDiS {
if (oldElInfo == nullptr)
delete elInfo;
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ Parallel::mpi::globalAdd(value);
+#endif
return value;
}
diff --git a/AMDiS/src/GenericOperatorTerm.h b/AMDiS/src/GenericOperatorTerm.h
index d703d1e6b39a5113a553ba754c97980aac8ac0ab..6114c21ab79eb7926162e46921765d2bdaa54680 100644
--- a/AMDiS/src/GenericOperatorTerm.h
+++ b/AMDiS/src/GenericOperatorTerm.h
@@ -315,7 +315,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
l1lt(grdLambda, LALt[iq], term(iq));
}
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -336,7 +336,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
}
}
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result)
{
@@ -521,25 +521,6 @@ inline void addZOT(Operator& op, double term)
addZOT(op, constant(term));
}
-
-// .............................................................................
-// Rosenbrock method
-
-// template <typename Term>
-// void addZOT_(Operator* op, const Term& term)
-// {
-// if (op->getUhOld()) {
-// addZOT(op, jacobian(term,0) * ...
-//
-// }
-
-// template <typename Term>
-// void addZOT(Operator& op, const Term& term)
-// {
-// op.addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
-// }
-
-
// _____________________________________________________________________________
// first order term using FirstOrderTerm::l1
@@ -775,7 +756,7 @@ inline expressions::Function1<expressions::Wrapper<TOut,WorldVector<double> >, e
eval(AbstractFunction<TOut, WorldVector<double> >* fct) { return function_(wrap(fct), X()); }
/// Integrate an expression over a domain.
-/** IF the expression does not contain any DOFVector, a mesh must be given as second argument */
+/** If the expression does not contain any DOFVector, a mesh must be given as second argument */
template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
integrate(Term term, Mesh* mesh_ = NULL);
diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc
index da9725d28ce1dd612512d0e05d34fd7fd95a07d7..beae250a82f4080d4e8410b99f08b2dadb9e4574 100644
--- a/AMDiS/src/Global.cc
+++ b/AMDiS/src/Global.cc
@@ -188,7 +188,11 @@ namespace AMDiS {
PRINT_LINE((*error), buff);
va_end(arg);
#if defined HAVE_PARALLEL_DOMAIN_AMDIS && !defined HAVE_PARALLEL_MTL4 && (DEBUG == 0 || defined NDEBUG)
+ #if (PETSC_VERSION_MINOR >= 5)
+ PetscError(MPI_COMM_WORLD, __LINE__, "Msg::print_error_exit", "Global.cc", 1, PETSC_ERROR_INITIAL, buff);
+ #else
PetscError(MPI_COMM_WORLD, __LINE__, "Msg::print_error_exit", "Global.cc", "src/", 1, PETSC_ERROR_INITIAL, buff);
+ #endif
#else
throw std::runtime_error(buff);
#endif
diff --git a/AMDiS/src/LinearAlgebra.h b/AMDiS/src/LinearAlgebra.h
new file mode 100644
index 0000000000000000000000000000000000000000..7b29384f524729b1ec4957db0f3246f69ca3b144
--- /dev/null
+++ b/AMDiS/src/LinearAlgebra.h
@@ -0,0 +1,239 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+ /// This file is work in progress and may replace the assembling now part of
+ /// DOFMatrix directly. The idea is to assemble matrices to a data-structure
+ /// provided by a linear algebra backend. You have to specify how an element-
+ /// matrix can be added to the global matrix, how a Dirichlet-row is inserted
+ /// and how the matrix is initialized/finalized.
+
+
+/** \file LinearAlgebra.h */
+
+#ifndef AMDIS_LINEAR_ALGEBRA_H
+#define AMDIS_LINEAR_ALGEBRA_H
+
+#include "utility/calculate_nnz.hpp"
+
+namespace AMDiS {
+
+ namespace tag
+ {
+ struct distributed {};
+ struct non_distributed {};
+
+ // backends
+ struct mtl {};
+ struct pmtl {};
+ struct petsc {};
+ struct p_petsc {};
+ }
+
+ namespace traits
+ {
+ template< class T >
+ struct distributed_tag : boost::mpl::if_<
+ mtl::traits::is_distributed< T >,
+ tag::distributed,
+ tag::non_distributed
+ > {};
+ }
+
+ namespace detail
+ {
+ template< class Dummy = void >
+ struct DefaultParameters {};
+ }
+
+ /** fill for each backend a linear algebra class that
+ * contains the following typedefs and functions
+ **/
+ template< class Backend = tag::mtl, class Value = double, class Parameters = detail::DefaultParameters<> >
+ struct LinearAlgebra
+ {
+ public: // typedefs
+
+ typedef Value ValueType;
+ typedef unsigned int SizeType;
+
+ typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, SizeType> MatrixPara;
+ typedef mtl::matrix::compressed2D<ValueType, MatrixPara> MatrixBase;
+ typedef mtl::vector::dense_vector<ValueType> VectorBase;
+
+ typedef mtl::matrix::inserter<MatrixBase, mtl::operations::update_plus<ValueType> > Inserter;
+ typedef MPI::Intracomm Communicator;
+
+ static const int minSlotSize = 20;
+
+ public: // methods
+
+ /// constructor
+ LinearAlgebra(Communicator& comm)
+ : comm(comm), ins(NULL)
+ { }
+
+
+ /// prepare matrix for insertion
+ template< class M >
+ void init_matrix(MatrixBase& matrix, MapperBase<M>& mapper, bool add_to_matrix = false)
+ {
+ int slot_size = estimate_slot_size(matrix);
+ prepare_matrix(matrix, mapper, traits::distributed_tag<MatrixBase>::type() );
+
+ if (!add_to_matrix)
+ set_to_zero(matrix);
+
+ // prepare the inserter
+ if (ins) {
+ delete ins;
+ ins = NULL;
+ }
+ ins = new Inserter(matrix, slot_size);
+ }
+
+
+ /// finish insertion
+ template< class M >
+ void exit_matrix(MatrixBase& matrix, MapperBase<M>& mapper)
+ {
+ // destroy inserter, to start insertion
+ assert( ins );
+ delete ins;
+ inserter = NULL;
+ }
+
+
+ /// insert values of element_matrix into matrix
+ template< class ElementMatrix, class VectorType >
+ void add_element_matrix(MatrixBase& matrix, const ElementMatrix& elMatrix,
+ const VectorType& rowIndices, const VectorType& colIndices)
+ {
+ assert( ins );
+ *ins << mtl::matrix::element_matrix(elMatrix, rowIndices, colIndices);
+ }
+
+ protected: // methods
+
+ /// resize matrix to appropriate size
+ template< class M >
+ void prepare_matrix(MatrixBase& matrix, MapperBase<M>& mapper, tag::non_distributed)
+ {
+ matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
+ }
+
+
+ /// estimate the number of nonzeros per row
+ int estimate_slot_size(const MatrixBase& matrix)
+ {
+ int nnzPerRow = 0;
+ if (num_rows(matrix) != 0)
+ nnzPerRow = static_cast<int>((matrix.nnz() * 1.2) / num_rows(matrix));
+ return std::max( nnzPerRow, minSlotSize );
+ }
+
+ private: // member variables
+ Communicator& comm;
+ Inserter* ins;
+ };
+
+#ifdef HAVE_PARALLEL_MTL4
+ template< class Value, class Parameters >
+ struct LinearAlgebra< tag::pmtl, Value, Parameters > : public LinearAlgebra< tag::mtl, Value, Parameters >
+ {
+ typedef LinearAlgebra< tag::mtl, Value, Parameters > super;
+
+ /// constructor
+ LinearAlgebra(Communicator& comm)
+ : super(comm)
+ { }
+
+ protected:
+ // eventuell als freie funktion schreiben
+ template< class M >
+ void prepare_matrix(MatrixBase& matrix, MapperBase<M>& mapper, tag::distributed)
+ {
+ mtl::par::block_distribution dist(mapper.getNumRows());
+ dist.setup_from_local_size(mapper.getMap().getLocalDofs());
+ matrix.change_dim(0, 0);
+ matrix.init_distribution(dist, dist, mapper.getNumRows(), mapper.getNumRows());
+ }
+ };
+#endif
+
+
+#ifdef HAVE_SEQ_PETSC
+ template< class Value, class Parameters >
+ struct LinearAlgebra< tag::petsc, Value, Parameters >
+ {
+ public: // typedefs
+
+ typedef PetscValue ValueType;
+ typedef PetscInt SizeType;
+
+ typedef PetscMatrix MatrixBase;
+ typedef PetscVector VectorBase;
+
+ typedef MPI::Intracomm Communicator;
+
+ public: // methods
+
+ /// constructor
+ LinearAlgebra(Communicator& comm)
+ : comm(comm)
+ { }
+
+
+ /// prepare matrix for insertion
+ template< class M >
+ void init_matrix(MatrixBase& mat, MapperBase<M>& mapper, bool add_to_matrix = false)
+ {
+ std::vector<size_t> nnz_per_row(mapper.getNumRows());
+ calculate_nnz(dof_matrix /* ?? */, nnz_per_row); // siehe auch MatrixNnzStructure fuer parallele NNZ
+
+ // eventuell nur max_nnz angeben
+ MatCreateSeqAIJ(comm, mapper.getNumRows(), mapper.getNumCols(), 0, &(nnz_per_row[0]), &mat.matrix);
+ }
+
+
+ /// finish insertion
+ template< class M >
+ void exit_matrix(MatrixBase& mat, MapperBase<M>& mapper)
+ {
+ MatAssemblyBegin(mat.matrix, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat.matrix, MAT_FINAL_ASSEMBLY);
+ }
+
+
+ /// insert values of element_matrix into matrix
+ template< class ElementMatrix, class VectorType >
+ void add_element_matrix(MatrixBase& mat, const ElementMatrix& elMatrix,
+ const VectorType& rowIndices, const VectorType& colIndices)
+ {
+ // delete dirichlet-dofs by changing the sign of these indices
+ MatSetValues(mat.matrix, rowIndices.size(), &rowIndices[0], colIndices.size(), &colIndices[0], elMatrix.address_data(), ADD_VALUES);
+ }
+
+ private: // member variables
+ Communicator& comm;
+ };
+#endif
+}
+
+#endif // AMDIS_DOFMATRIX_H
diff --git a/AMDiS/src/ProblemInstat.cc b/AMDiS/src/ProblemInstat.cc
index 8b2f9df4bbb88aba879c8c49aad372bc0563b66a..43e51a434139c376ef8f2e2d840b7ec1191d5de5 100644
--- a/AMDiS/src/ProblemInstat.cc
+++ b/AMDiS/src/ProblemInstat.cc
@@ -89,13 +89,16 @@ namespace AMDiS {
ProblemStatBase &initialProb)
: ProblemInstatBase(sname, &initialProb),
problemStat(&prob),
- oldSolution(nullptr)
+ oldSolution(NULL)
{}
ProblemInstat::~ProblemInstat()
{
- delete oldSolution;
+ if (oldSolution) {
+ delete oldSolution;
+ oldSolution = NULL;
+ }
}
@@ -153,7 +156,8 @@ namespace AMDiS {
lastTimepoint = MPI::Wtime();
#endif
- oldSolution->copy(*(problemStat->getSolution()));
+ if (oldSolution)
+ oldSolution->copy(*(problemStat->getSolution()));
}
}
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index e17ef11c2e698514892e6766dc80f6aeae8964c4..1e767bcccf8c7c2c1369ba3019b6d9a22d2fe0cd 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -676,7 +676,7 @@ namespace AMDiS {
std::string componentName = "";
Parameters::get(numberedName + "->name", componentName);
std::vector<std::string> names;
- if (name != "")
+ if (componentName != "")
names.push_back(componentName);
std::vector<int> comp;
Parameters::get(numberedName + "->components", comp);
diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h
index 3b285a05a8419dbfebd70ae7e19b72f5c1dc4548..0ba96dc117f7e859b827cb5c60c59fc06ba1422a 100644
--- a/AMDiS/src/SecondOrderTerm.h
+++ b/AMDiS/src/SecondOrderTerm.h
@@ -140,7 +140,7 @@ namespace AMDiS {
l1lt(grdLambda, LALt[iq], (*factor));
}
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -161,7 +161,7 @@ namespace AMDiS {
}
}
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result)
{
@@ -200,7 +200,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -208,7 +208,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -257,7 +257,7 @@ namespace AMDiS {
lalt(Lambda, matrix, LALt[iq], symmetric, 1.0);
}
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -265,7 +265,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -317,7 +317,7 @@ namespace AMDiS {
lalt_kl(grdLambda, xi, xj, LALt[iq], (*factor));
}
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -330,7 +330,7 @@ namespace AMDiS {
result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac;
}
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result)
{
@@ -376,7 +376,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -384,7 +384,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -425,7 +425,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -433,7 +433,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -475,7 +475,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -483,7 +483,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -527,7 +527,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -535,7 +535,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -584,7 +584,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -592,7 +592,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -634,7 +634,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -642,7 +642,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -690,7 +690,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -698,7 +698,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -791,7 +791,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -799,7 +799,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -839,7 +839,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -847,7 +847,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -889,7 +889,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -897,7 +897,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
@@ -947,7 +947,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
- /// Implenetation of SecondOrderTerm::eval().
+ /// Implemetation of SecondOrderTerm::eval().
void eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -955,7 +955,7 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double factor);
- /// Implenetation of SecondOrderTerm::weakEval().
+ /// Implemetation of SecondOrderTerm::weakEval().
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result);
diff --git a/AMDiS/src/expressions/cmath_expr.hpp b/AMDiS/src/expressions/cmath_expr.hpp
index f99b978e6c566d1be64fbce4943ad8e1581cd970..ccfa46eea201f1d1186a9c96222eefcdca5b5c95 100644
--- a/AMDiS/src/expressions/cmath_expr.hpp
+++ b/AMDiS/src/expressions/cmath_expr.hpp
@@ -65,7 +65,7 @@ namespace AMDiS
int getDegree() const
{
- return 0;
+ return 2*super::term.getDegree();
}
inline value_type operator()(const int& iq) const { return (super::term(iq) > 0.0 ? 1.0 : -1.0); }
@@ -83,7 +83,7 @@ namespace AMDiS
Ceil(const Term& term_) : super(term_) {}
- int getDegree() const { return 0; }
+ int getDegree() const { return super::term.getDegree(); }
inline value_type operator()(const int& iq) const { return std::ceil(super::term(iq)); }
@@ -100,7 +100,7 @@ namespace AMDiS
Floor(const Term& term_) : super(term_) {}
- int getDegree() const { return 0; }
+ int getDegree() const { return super::term.getDegree(); }
inline value_type operator()(const int& iq) const { return std::floor(super::term(iq)); }
diff --git a/AMDiS/src/expressions/expressions.h b/AMDiS/src/expressions/expressions.h
index 78e27b953c332210f1178d832cd5b09d63cb30cb..7f27ea1ee68ad469e629ff7aa62c7ab0abb5a04e 100644
--- a/AMDiS/src/expressions/expressions.h
+++ b/AMDiS/src/expressions/expressions.h
@@ -30,6 +30,7 @@
#include "coords_expr.hpp" // coordinates / normal vectors
#include "valueOf.hpp" // value of DOFVector at QPs
#include "gradientOf.hpp" // gradient of DOFVector at QPs
+#include "hessianOf.hpp" // second derivatives of DOFVector at QPs
#include "add_expr.hpp" // add two expressions
#include "mult_expr.hpp" // multiply two expressions
diff --git a/AMDiS/src/expressions/hessianOf.hpp b/AMDiS/src/expressions/hessianOf.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1bb19bcca07caf8ae6f907605ab48be4151a97b2
--- /dev/null
+++ b/AMDiS/src/expressions/hessianOf.hpp
@@ -0,0 +1,234 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+
+/** \file hessianOf.hpp */
+
+#ifndef AMDIS_HESSIAN_OF_HPP
+#define AMDIS_HESSIAN_OF_HPP
+
+#include "AMDiS_fwd.h"
+#include "LazyOperatorTerm.h"
+#include "DOFVector.h"
+#include "traits/category.hpp"
+
+namespace AMDiS
+{
+ namespace expressions
+ {
+ /// Expressions that extracts the gradient of a DOFVector at QPs
+ template<typename Vector, typename Name>
+ struct HessianOf : public LazyOperatorTermBase
+ {
+ typedef typename traits::category<Vector>::value_type T;
+ typedef typename D2Type<T>::type value_type;
+ typedef Name id;
+
+ DOFVector<T>* vecDV;
+ mutable mtl::dense_vector<value_type> vec;
+ mutable mtl::dense_vector<T> coeff;
+
+ HessianOf(Vector& vector) : vecDV(&vector) {}
+ HessianOf(Vector* vector) : vecDV(vector) {}
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const
+ {
+ feSpaces.insert(vecDV->getFeSpace());
+ }
+
+ int getDegree() const
+ {
+ return vecDV->getFeSpace()->getBasisFcts()->getDegree() /* -1 */;
+ }
+
+ template<typename OT>
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
+ { FUNCNAME("HessianOf::initElement");
+
+ if (ot && subAssembler) {
+ ERROR_EXIT("Hessian expression not yet implemented for operator-terms!\n");
+ } else if (quad)
+ vecDV->getD2AtQPs(elInfo, quad, NULL, vec);
+ else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ for (size_t i = 0; i < nBasisFct; i++)
+ localBasisFct->evalD2Uh(*basisFct->getCoords(i), grdLambda, coeff, &vec[i]);
+ }
+ }
+
+
+ template<typename OT>
+ inline void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
+ { FUNCNAME("HessianOf::initElement");
+
+ ERROR_EXIT("Hessian expression not yet implemented for Dual-mesh!\n");
+ }
+
+ value_type operator()(const int& iq) const { return vec[iq]; }
+
+ std::string str() const { return std::string("hessian(") + vecDV->getName() + ")"; }
+ };
+
+
+ /// Expressions that extracts the partial derivative of a DOFVector at QPs
+ template<typename Vector, typename Name>
+ struct LaplacianOf : public LazyOperatorTermBase
+ {
+ typedef typename traits::category<Vector>::value_type T;
+ typedef T value_type;
+ typedef Name id;
+
+ DOFVector<T>* vecDV;
+ mutable mtl::dense_vector<typename D2Type<T>::type> vec_tmp;
+ mutable mtl::dense_vector<value_type> vec;
+ mutable mtl::dense_vector<T> coeff;
+ int comp;
+
+ LaplacianOf(Vector& vector) : vecDV(&vector) { }
+ LaplacianOf(Vector* vector) : vecDV(vector) { }
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const
+ {
+ feSpaces.insert(vecDV->getFeSpace());
+ }
+
+ int getDegree() const
+ {
+ return vecDV->getFeSpace()->getBasisFcts()->getDegree() /* -1 */;
+ }
+
+ template<typename OT>
+ void initElement(OT* ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
+ { FUNCNAME("LaplacianOf::initElement");
+
+ if (ot && subAssembler) {
+ ERROR_EXIT("Laplacian expression not yet implemented for operator-terms!\n");
+ } else if (quad) {
+ vecDV->getD2AtQPs(elInfo, quad, NULL, vec_tmp);
+ for (size_t i = 0; i < size(vec_tmp); i++) {
+ vec[i] = vec_tmp[i][0][0];
+ for (int j = 1; j < Global::getGeo(WORLD); j++)
+ vec[i] += vec_tmp[i][j][j];
+ }
+ } else if (basisFct) {
+ const BasisFunction *localBasisFct = vecDV->getFeSpace()->getBasisFcts();
+
+ // get coefficients of DOFVector
+ coeff.change_dim(localBasisFct->getNumber());
+ vecDV->getLocalVector(elInfo->getElement(), coeff);
+
+ // eval basisfunctions of DOFVector at coords of given basisFct
+ size_t nBasisFct = basisFct->getNumber();
+ vec.change_dim(nBasisFct);
+ WorldMatrix<double> hessian;
+
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ for (size_t i = 0; i < nBasisFct; i++) {
+ localBasisFct->evalD2Uh(*basisFct->getCoords(i), grdLambda, coeff, &hessian);
+ vec[i] = hessian[0][0];
+ for (int j = 1; j < Global::getGeo(WORLD); j++)
+ vec[i] += hessian[j][j];
+ }
+ }
+ }
+
+ template<typename OT>
+ void initElement(OT* ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad,
+ const BasisFunction *basisFct = NULL)
+ { FUNCNAME("LaplacianOf::initElement");
+
+ ERROR_EXIT("Laplacian expression not yet implemented for Dual-mesh!\n");
+ }
+
+ value_type operator()(const int& iq) const { return vec[iq]; }
+
+ std::string str() const { return std::string("laplace(") + vecDV->getName() + ")"; }
+ };
+
+ } // end namespace expressions
+
+
+ // gradient of a DOFVector
+ // _____________________________________________________________________________
+
+ // with Name
+ template<typename Name, typename T>
+ expressions::HessianOf<DOFVector<T>, Name > hessianOf(DOFVector<T>& vector)
+ { return expressions::HessianOf<DOFVector<T>, Name >(vector); }
+
+ template<typename Name, typename T>
+ expressions::HessianOf<DOFVector<T>, Name > hessianOf(DOFVector<T>* vector)
+ { return expressions::HessianOf<DOFVector<T>, Name >(vector); }
+
+ // without Name
+ template<typename T>
+ expressions::HessianOf<DOFVector<T>, _unknown > hessianOf(DOFVector<T>& vector)
+ { return expressions::HessianOf<DOFVector<T>, _unknown >(vector); }
+
+ template<typename T>
+ expressions::HessianOf<DOFVector<T>, _unknown > hessianOf(DOFVector<T>* vector)
+ { return expressions::HessianOf<DOFVector<T>, _unknown >(vector); }
+
+
+ // Partial derivative of a DOFVector
+ // _____________________________________________________________________________
+
+ // with Name
+ template<typename Name, typename T>
+ expressions::LaplacianOf<DOFVector<T>, Name > laplacianOf(DOFVector<T>& vector)
+ { return expressions::LaplacianOf<DOFVector<T>, Name >(vector); }
+
+ template<typename Name, typename T>
+ expressions::LaplacianOf<DOFVector<T>, Name > laplacianOf(DOFVector<T>* vector)
+ { return expressions::LaplacianOf<DOFVector<T>, Name >(vector); }
+
+ // without Name
+ template<typename T>
+ expressions::LaplacianOf<DOFVector<T>, _unknown > laplacianOf(DOFVector<T>& vector)
+ { return expressions::LaplacianOf<DOFVector<T>, _unknown >(vector); }
+
+ template<typename T>
+ expressions::LaplacianOf<DOFVector<T>, _unknown > laplacianOf(DOFVector<T>* vector)
+ { return expressions::LaplacianOf<DOFVector<T>, _unknown >(vector); }
+
+} // end namespace AMDiS
+
+
+#endif // AMDIS_HESSIAN_OF_HPP
diff --git a/AMDiS/src/parallel/MatrixNnzStructure.cc b/AMDiS/src/parallel/MatrixNnzStructure.cc
index dbfb1d45d79fe71d99e1a850a26889a17992658c..0c514c7bf3f77348d9ca1c589a8dc4f30878340e 100644
--- a/AMDiS/src/parallel/MatrixNnzStructure.cc
+++ b/AMDiS/src/parallel/MatrixNnzStructure.cc
@@ -58,11 +58,11 @@ namespace AMDiS { namespace Parallel {
{
FUNCNAME("MatrixNnzStructure::create()");
- int nRankRows = rowDofMap.getRankDofs();
- int rankStartRowIndex = rowDofMap.getStartDofs();
+ int nRankRows = rowDofMap.getRankDofs(); // Number of DOFs owned by rank.
+ int rankStartRowIndex = rowDofMap.getStartDofs(); // Smallest global index of a DOF owned by the rank.
int nRankCols = colDofMap.getRankDofs();
- int nOverallCols = colDofMap.getOverallDofs();
+ int nOverallCols = colDofMap.getOverallDofs(); // Number of global DOFs (this value is thus the same on all ranks).
int rankStartColIndex = colDofMap.getStartDofs();
create(nRankRows, (!localMatrix ? nRankRows : -1));
@@ -139,7 +139,7 @@ namespace AMDiS { namespace Parallel {
// === Prepare MTL4 iterators for the current DOF matrix. ===
- Matrix baseMat = dofMat->getBaseMatrix();
+ Matrix baseMat = dofMat->getBaseMatrix(); // TBD: warum hier keine Referenz?
traits::col<Matrix>::type getCol(baseMat);
traits::row<Matrix>::type getRow(baseMat);
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index dfbb9dc838a8e621ac1eeb653b2909652c16bd29..d57739321f31c65dcf5bf605b428ecdd2bfe2b0f 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -1294,7 +1294,7 @@ namespace AMDiS { namespace Parallel {
getImbalanceFactor(imbalanceFactor, minDofs, maxDofs, sumDofs);
if (mpiRank == 0)
- MSG("Imbalancing factor: %.2f [ minDofs = %d, maxDofs = %d, sumDofs = %d ]\n",
+ MSG("Imbalancing factor: %.2f %% [ minDofs = %d, maxDofs = %d, sumDofs = %d ]\n",
imbalanceFactor * 100.0, minDofs, maxDofs, sumDofs);
}
@@ -1447,13 +1447,11 @@ namespace AMDiS { namespace Parallel {
int repartitioning = 0;
double imbalanceFactor = getImbalanceFactor();
+ double imbalanceRepartitionBound = 0.2;
+ Parameters::get("parallel->repartitioning->imbalance",
+ imbalanceRepartitionBound);
if (mpiRank == 0) {
- double imbalanceRepartitionBound = 0.2;
- Parameters::get("parallel->repartitioning->imbalance",
- imbalanceRepartitionBound);
-
-
if (imbalanceFactor > imbalanceRepartitionBound)
repartitioning = 1;
@@ -1461,6 +1459,11 @@ namespace AMDiS { namespace Parallel {
} else {
mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
}
+#if (DEBUG != 0)
+ if (repartitioning == 0) {
+ MSG("imbalanceFactor = %f < %f = imbalanceRepartitionBound\n", imbalanceFactor, imbalanceRepartitionBound);
+ }
+#endif
return (repartitioning != 0);
}
diff --git a/AMDiS/src/parallel/PITL_Solver.h b/AMDiS/src/parallel/PITL_Solver.h
index 905f9acfb0c01fd06f6de175301b0e2a7332d1cc..0836c788f2b504f4483b6ebc030297601e37d670 100644
--- a/AMDiS/src/parallel/PITL_Solver.h
+++ b/AMDiS/src/parallel/PITL_Solver.h
@@ -37,10 +37,10 @@ namespace AMDiS
using namespace MTLTypes;
template< typename SolverType >
- struct PITL_Solver : MTL4Solver< PMTLMatrix, PMTLVector, ITL_Runner< SolverType, PMTLMatrix, PMTLVector >, true >
+ struct PITL_Solver : MTL4Solver< PMTLMatrix, PMTLVector, ITL_Runner< SolverType, PMTLMatrix, PMTLVector > >
{
PITL_Solver(std::string name)
- : MTL4Solver< PMTLMatrix, PMTLVector, ITL_Runner< SolverType, PMTLMatrix, PMTLVector >, true >(name) {}
+ : MTL4Solver< PMTLMatrix, PMTLVector, ITL_Runner< SolverType, PMTLMatrix, PMTLVector > >(name) {}
};
typedef PITL_Solver< cg_solver_type > P_CGSolver;
diff --git a/AMDiS/src/parallel/PMTL4Solver.h b/AMDiS/src/parallel/PMTL4Solver.h
deleted file mode 100644
index dfbfb47338149c870725562e566d11f800530a97..0000000000000000000000000000000000000000
--- a/AMDiS/src/parallel/PMTL4Solver.h
+++ /dev/null
@@ -1,114 +0,0 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-
-/** \file PMtl4Solver.h */
-
-#ifndef AMDIS_PMTL4_SOLVER_H
-#define AMDIS_PMTL4_SOLVER_H
-
-#ifdef HAVE_PARALLEL_MTL4
-
-#include "solver/MTL4SolverBase.h"
-#include "parallel/ParallelSolver.h"
-
-namespace AMDiS {
-
- namespace Parallel {
-
-
- /** \ingroup Solver
- *
- * \brief
- * Wrapper class for different PMTL4 solvers. These solvers
- * are parametrized by Matrix- and VectorType. The different
- * algorithms, like krylov subspace methods or other external
- * solvers where PMTL4 provides an interface, can be assigned
- * by different Runner objects.
- **/
- template< typename MatrixType, typename VectorType, typename Runner >
- class PMTL4Solver : public ParallelSolver, protected MTL4SolverBase<MatrixType, VectorType, Runner>
- {
- public:
- typedef PMTL4Solver<MatrixType, VectorType, Runner> self;
- typedef MTL4SolverBase<MatrixType, VectorType, Runner> details;
-
- using details::solve;
-
- /// Creator class used in the LinearSolverMap.
- class Creator : public LinearSolverCreator
- {
- public:
- virtual ~Creator() {}
-
- /// Returns a new MTL4Solver object.
- LinearSolver* create()
- {
- return new self(this->name);
- }
- };
-
- /// Constructor
- PMTL4Solver(std::string n)
- : ParallelSolver(n, false),
- details(this)
- {}
-
-
- /// Implementation of \ref LinearSolver::getRunner()
- virtual OEMRunner* getRunner()
- {
- return details::getRunner();
- }
-
-
- /// Implementation of \ref LinearSolver::getLeftPrecon()
- virtual OEMPreconditioner* getLeftPrecon()
- {
- return details::getLeftPrecon();
- }
-
-
- /// Implementation of \ref LinearSolver::getRightPrecon()
- virtual OEMPreconditioner* getRightPrecon()
- {
- return details::getRightPrecon();
- }
-
- protected:
- /// Implementation of \ref LinearSolver::solveLinearSystem()
- virtual int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
- SystemVector& x,
- SystemVector& b,
- bool createMatrixData,
- bool storeMatrixData)
- {
- ParallelMapper mapper(*ParallelSolver::getDofMapping());
- return details::solve(A, x, b, mapper, createMatrixData, storeMatrixData);
- }
- };
-
- } // end namespace Parallel
-} // end namespace AMDiS
-
-#endif // HAVE_PARALLEL_MTL4
-
-#endif // AMDIS_PMTL4_SOLVER_H
-
diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc b/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc
index 88cb0def56bc98903f707de65338d8632955d26a..e372bd8484e5b9ba1bf7903aeadb34d731395a45 100644
--- a/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc
+++ b/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc
@@ -84,6 +84,7 @@ namespace AMDiS { namespace Parallel {
}
+ /// initializes vecSol, vecRhs, mat und componentIthCoarseMap
void ParallelCoarseSpaceSolver::prepare()
{
FUNCNAME("ParallelCoarseSpaceSolver:prepare()");
diff --git a/AMDiS/src/parallel/PetscHelper.cc b/AMDiS/src/parallel/PetscHelper.cc
index f49ba98981dbc5c61df82a772914cf0f8da5b151..d053a5096bc6063b92484750e12b6d85acf222e3 100644
--- a/AMDiS/src/parallel/PetscHelper.cc
+++ b/AMDiS/src/parallel/PetscHelper.cc
@@ -327,6 +327,25 @@ namespace AMDiS
rtol, atol, maxIt);
}
+
+ void createSolver(MPI::Intracomm comm, KSP &ksp, Mat m, std::string kspPrefix, int info)
+ {
+ KSPCreate(comm, &ksp);
+ #if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, m, m);
+ #else
+ KSPSetOperators(ksp, m, m, SAME_NONZERO_PATTERN);
+ #endif
+ KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
+ KSPSetType(ksp, KSPBCGS);
+ KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
+
+ if (info >= 10)
+ KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
+ else if (info >= 20)
+ KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
+ }
+
} // end namespace petsc_helper
} // end namespace Parallel
} // end namespace AMDiS
diff --git a/AMDiS/src/parallel/PetscHelper.h b/AMDiS/src/parallel/PetscHelper.h
index b9f8a29a210f1d15ba8c4f2bc87b37ff828f490a..d102f463abca8f0c3c29344be5b87dc87030a829 100644
--- a/AMDiS/src/parallel/PetscHelper.h
+++ b/AMDiS/src/parallel/PetscHelper.h
@@ -115,6 +115,8 @@ namespace AMDiS
PetscReal rtol = PETSC_DEFAULT,
PetscReal atol = PETSC_DEFAULT,
PetscInt maxIt = PETSC_DEFAULT);
+
+ void createSolver(MPI::Intracomm comm, KSP &ksp, Mat m, std::string kspPrefix = "", int info = 0);
} // end namespace petsc_helper
} // end namespace Parallel
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index f3e666df3443bb95849d1c85eef1b0e8dbf09540..3f00e4e4503b260503647943796b863c85cf3373 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -74,12 +74,13 @@ namespace AMDiS { namespace Parallel {
super::init(fe0, fe1, createGlobalMapping);
}
- int PetscSolver::solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+ /// see seq::MTL4Solver::solveLinearSystem(...)
+ int PetscSolver::solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
bool createMatrixData,
bool storeMatrixData)
- {
+ {
FUNCNAME("PetscSolver::solveLinearSystem()");
TEST_EXIT(meshDistributor)("No meshDistributor provided. Should not happen!\n");
@@ -103,7 +104,7 @@ namespace AMDiS { namespace Parallel {
}
return 0;
- }
+ }
void PetscSolver::fillPetscMatrix(DOFMatrix* mat)
{
@@ -113,6 +114,7 @@ namespace AMDiS { namespace Parallel {
}
+ /// see seq::PetscSolver::solve(const MatrixType& A, VectorType& x, const VectorType& b)
void PetscSolver::solve(Vec &rhs, Vec &sol)
{
PetscErrorCode solverError = KSPSolve(kspInterior, rhs, sol);
@@ -135,6 +137,7 @@ namespace AMDiS { namespace Parallel {
}
+ /// was soll das denn genau machen?
void PetscSolver::solveGlobal(Vec &rhs, Vec &sol)
{
FUNCNAME("PetscSolver::solveGlobal()");
@@ -149,6 +152,7 @@ namespace AMDiS { namespace Parallel {
}
+ /// -> Helper-Funktion, benötigt nur domainComm
void PetscSolver::copyVec(Vec& originVec, Vec& destVec,
std::vector<int>& originIndex, std::vector<int>& destIndex)
{
@@ -178,6 +182,7 @@ namespace AMDiS { namespace Parallel {
}
+ /// -> Helper-Funktion
bool PetscSolver::testMatrixSymmetric(Mat mat, bool advancedTest)
{
FUNCNAME("PetscSolver::testMatrixSymmetric()");
diff --git a/AMDiS/src/parallel/PetscSolverCahnHilliard.cc b/AMDiS/src/parallel/PetscSolverCahnHilliard.cc
index 76fb678d743b2c4f14edff295d85aac69ec5224b..ae6352b2b9a9b05df8268fc371f2e114d4d1bca4 100644
--- a/AMDiS/src/parallel/PetscSolverCahnHilliard.cc
+++ b/AMDiS/src/parallel/PetscSolverCahnHilliard.cc
@@ -85,7 +85,11 @@ namespace AMDiS { namespace Parallel {
{
// Create FGMRES based outer solver
KSPCreate(meshDistributor->getMpiComm(0), &ksp);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior());
+#else
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
if (getInfo() >= 10)
KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
else if (getInfo() >= 20)
diff --git a/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc b/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc
index e71ac593c08a4e214c8011362de6b26da994d694..f901be28b1fd4eadf6b74181114345463df7898f 100644
--- a/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc
+++ b/AMDiS/src/parallel/PetscSolverCahnHilliard2.cc
@@ -67,7 +67,11 @@ namespace AMDiS { namespace Parallel {
/// create new solver for S
KSP kspS;
KSPCreate(*data->mpiCommGlobal, &kspS);
- KSPSetOperators(kspS, S, S, SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(kspS, S, S);
+#else
+ KSPSetOperators(kspS, S, S, SAME_NONZERO_PATTERN);
+#endif
petsc_helper::setSolver(kspS, "S_", KSPFGMRES, PCSHELL, 1e-6, 1e-8, 1);
{
PC pc;
@@ -112,8 +116,12 @@ namespace AMDiS { namespace Parallel {
// Create FGMRES based outer solver
MSG("CREATE POS 1: %p\n", &ksp);
- KSPCreate(domainComm, &ksp);
+ KSPCreate(domainComm, &ksp);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior());
+#else
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
if (getInfo() >= 10)
KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
else if (getInfo() >= 20)
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index e1e81b3abd532c20224c8524d1b57fa76870818f..ef54c0621e9e36079aea14265bda8da76211613f 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -1050,7 +1050,11 @@ namespace AMDiS { namespace Parallel {
}
KSPCreate(domainComm, &ksp_schur_primal);
- KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal);
+#else
+ KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
KSPSetType(ksp_schur_primal, KSPGMRES);
KSPSetFromOptions(ksp_schur_primal);
@@ -1070,8 +1074,11 @@ namespace AMDiS { namespace Parallel {
// === Create KSP solver object and set appropriate solver options. ===
KSPCreate(domainComm, &ksp_schur_primal);
- KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal);
+#else
+ KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
KSPSetType(ksp_schur_primal, KSPPREONLY);
PC pc_schur_primal;
@@ -1357,7 +1364,11 @@ namespace AMDiS { namespace Parallel {
}
KSPCreate(domainComm, &ksp_feti);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp_feti, mat_feti, mat_feti);
+#else
KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(ksp_feti, "feti_");
KSPSetType(ksp_feti, KSPGMRES);
KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
@@ -1431,7 +1442,11 @@ namespace AMDiS { namespace Parallel {
(void(*)(void))petscMultMatFetiInexact);
KSPCreate(domainComm, &ksp_feti);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp_feti, mat_feti, mat_feti);
+#else
KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(ksp_feti, "feti_");
KSPSetType(ksp_feti, KSPGMRES);
KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
@@ -1448,9 +1463,11 @@ namespace AMDiS { namespace Parallel {
createVec(localDofMap, fetiInexactPreconData.tmp_vec_b0);
KSPCreate(domainComm, &(fetiInexactPreconData.ksp_pc_feti));
- KSPSetOperators(fetiInexactPreconData.ksp_pc_feti,
- mat_lagrange, mat_lagrange,
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(fetiInexactPreconData.ksp_pc_feti, mat_lagrange, mat_lagrange);
+#else
+ KSPSetOperators(fetiInexactPreconData.ksp_pc_feti, mat_lagrange, mat_lagrange, SAME_NONZERO_PATTERN);
+#endif
KSPGetPC(fetiInexactPreconData.ksp_pc_feti,
&(fetiInexactPreconData.pc_feti));
createFetiPreconLumped(fetiInexactPreconData.pc_feti);
@@ -1582,8 +1599,11 @@ namespace AMDiS { namespace Parallel {
("Stokes mode does not yet support the Dirichlet precondition!\n");
KSPCreate(PETSC_COMM_SELF, &ksp_interior);
- KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior,
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior);
+#else
+ KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
KSPSetType(ksp_interior, KSPPREONLY);
PC pc_interior;
diff --git a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
index 858f7969168f15bee374b82d26d9160fa2cd5301..cbd5d6ccf9e8d48c0be216d18d54e04af5a06cee 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
@@ -150,8 +150,11 @@ namespace AMDiS { namespace Parallel {
FUNCNAME("PetscSolverGlobalBlockMatrix::initSolver()");
KSPCreate(domainComm, &ksp);
- KSPSetOperators(ksp, getMatInterior(), getMatInterior(),
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior());
+#else
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
KSPSetFromOptions(ksp);
}
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 5f72e8dbd75edbe47bb4a006ab680a963e70c961..6cb8273f6795b81607f4fc125e5a262e231cc89c 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -328,8 +328,11 @@ namespace AMDiS { namespace Parallel {
// === Create solver for the non primal (thus local) variables. ===
KSPCreate(domainComm, &kspInterior);
- KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(kspInterior, getMatInterior(), getMatInterior());
+#else
+ KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
KSPSetOptionsPrefix(kspInterior, "interior_");
KSPSetType(kspInterior, KSPPREONLY);
KSPGetPC(kspInterior, &pcInterior);
@@ -380,6 +383,11 @@ namespace AMDiS { namespace Parallel {
}
+ /// 1.) set startsolution
+ /// 2.) create null-space
+ /// 3.) solve Ax=b
+ /// 4.) destroy null-space
+ /// 5.) transfer solution back to DOFVector
void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec,
AdaptInfo *adaptInfo)
{
@@ -447,8 +455,11 @@ namespace AMDiS { namespace Parallel {
TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
MSG("Remove nullspace from rhs vector.\n");
-
+#if (PETSC_VERSION_MINOR >= 5)
+ MatNullSpaceRemove(matNullspace, getVecRhsInterior());
+#else
MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
+#endif
}
} else {
TEST_EXIT(removeRhsNullspace == false)
@@ -601,8 +612,11 @@ namespace AMDiS { namespace Parallel {
void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
{
KSPCreate(domainComm, &ksp);
- KSPSetOperators(ksp, getMatInterior(), getMatInterior(),
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior());
+#else
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetType(ksp, KSPBCGS);
KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
@@ -632,6 +646,7 @@ namespace AMDiS { namespace Parallel {
{ }
+ /// called by \ref fillPetscMatrix
void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
int rowComp, int colComp)
{
@@ -810,6 +825,7 @@ namespace AMDiS { namespace Parallel {
}
+ /// called by \ref fillPetscRhs
void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior,
Vec vecCoarse,
DOFVector<double>* vec,
diff --git a/AMDiS/src/parallel/PetscSolverNSCH.cc b/AMDiS/src/parallel/PetscSolverNSCH.cc
index cbeac6299597696270afba983959053843bc337a..2e6a142edaead8ba565f31363c8b2515e72edf25 100644
--- a/AMDiS/src/parallel/PetscSolverNSCH.cc
+++ b/AMDiS/src/parallel/PetscSolverNSCH.cc
@@ -184,7 +184,11 @@ namespace AMDiS { namespace Parallel {
// Create FGMRES based outer solver
MSG("CREATE POS 1: %p\n", &ksp);
KSPCreate(domainComm, &ksp);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior());
+#else
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCSHELL, getRelative(), getTolerance(), getMaxIterations());
setConstantNullSpace(ksp, componentSpaces[0]->getMesh()->getDim() , true);
@@ -333,7 +337,12 @@ namespace AMDiS { namespace Parallel {
///erstelle kspVelocity
KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspVelocity));
- KSPSetOperators(matShellContext.kspVelocity, matShellContext.velocityMat, matShellContext.velocityMat, SAME_NONZERO_PATTERN);
+
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(matShellContext.kspVelocity, matShellContext.velocityMat, matShellContext.velocityMat);
+#else
+ KSPSetOperators(matShellContext.kspVelocity, matShellContext.velocityMat, matShellContext.velocityMat, SAME_NONZERO_PATTERN);
+#endif
///regularisiere LaplaceMatrix
if (regularizeLaplace)
@@ -342,7 +351,11 @@ namespace AMDiS { namespace Parallel {
rows[0]=0;
MatZeroRows(laplaceMatrixSolver->getMatInterior(), 1, rows, 0, PETSC_NULL, PETSC_NULL);
KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspLaplace));
- KSPSetOperators(matShellContext.kspLaplace, laplaceMatrixSolver->getMatInterior(), laplaceMatrixSolver->getMatInterior(), SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(matShellContext.kspLaplace, laplaceMatrixSolver->getMatInterior(), laplaceMatrixSolver->getMatInterior());
+#else
+ KSPSetOperators(matShellContext.kspLaplace, laplaceMatrixSolver->getMatInterior(), laplaceMatrixSolver->getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
}
else
{ matShellContext.kspLaplace=laplaceMatrixSolver->getSolver();
diff --git a/AMDiS/src/parallel/PetscSolverNavierStokes.cc b/AMDiS/src/parallel/PetscSolverNavierStokes.cc
index 6509b5f9dfe0bb62bb73b851cf5c970dda24783c..4cb7f602c7747d110f5600fce06fd3bcd9d52957 100644
--- a/AMDiS/src/parallel/PetscSolverNavierStokes.cc
+++ b/AMDiS/src/parallel/PetscSolverNavierStokes.cc
@@ -117,7 +117,11 @@ namespace AMDiS { namespace Parallel {
MSG("CREATE POS 1: %p\n", &ksp);
KSPCreate(domainComm, &ksp);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp, getMatInterior(), getMatInterior());
+#else
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
if (getInfo() >= 10)
KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
else if (getInfo() >= 20)
diff --git a/AMDiS/src/parallel/PetscSolverSchur.cc b/AMDiS/src/parallel/PetscSolverSchur.cc
index ea6ff64c96e178cc15e116b6d51daecd6e10e506..0713cc8ffcfd0e906498bb3de68a53795e3c945b 100644
--- a/AMDiS/src/parallel/PetscSolverSchur.cc
+++ b/AMDiS/src/parallel/PetscSolverSchur.cc
@@ -292,8 +292,11 @@ namespace AMDiS { namespace Parallel {
KSPCreate(domainComm, &kspInterior);
- KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
- SAME_NONZERO_PATTERN);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(kspInterior, getMatInterior(), getMatInterior());
+#else
+ KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
+#endif
KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetFromOptions(kspInterior);
diff --git a/AMDiS/src/solver/BITL_Solver.h b/AMDiS/src/solver/BITL_Solver.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbffa4b9fd4028db9049f7bbe567e66ccd7eb83e
--- /dev/null
+++ b/AMDiS/src/solver/BITL_Solver.h
@@ -0,0 +1,101 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+/** \file BITL_Solver.h */
+
+#ifndef AMDIS_BITL_SOLVER_H
+#define AMDIS_BITL_SOLVER_H
+
+#include "solver/BlockMTLMatrix.h"
+
+#include "solver/MTL4Solver.h"
+#include "solver/ITL_Runner.h"
+#include "solver/ITL_Solver.h"
+#include "MTL4Types.h"
+
+#include "solver/itl/block_diagonal.hpp"
+
+namespace AMDiS {
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Wrapper for MTL4 itl-solvers applied to block matrices.
+ *
+ */
+ template< typename SolverType >
+ struct BITL_Solver : MTL4Solver< BlockMTLMatrix, MTLTypes::MTLVector, ITL_Runner< SolverType, BlockMTLMatrix, MTLTypes::MTLVector > >
+ {
+ BITL_Solver(std::string name)
+ : MTL4Solver< BlockMTLMatrix, MTLTypes::MTLVector, ITL_Runner< SolverType, BlockMTLMatrix, MTLTypes::MTLVector > >(name) {}
+ };
+
+ // ===================================================================================
+
+ typedef BITL_Solver< cg_solver_type > B_CGSolver;
+ typedef BITL_Solver< cgs_solver_type > B_CGSSolver;
+// typedef BITL_Solver< bicg_solver_type > B_BiCGSolver; // uses adjoint(A)
+ typedef BITL_Solver< bicgstab_type > B_BiCGStabSolver;
+// typedef BITL_Solver< bicgstab2_type > B_BiCGStab2Solver; // ERROR
+// typedef BITL_Solver< qmr_solver_type > B_QMRSolver; // uses trans(A)
+ typedef BITL_Solver< tfqmr_solver_type > B_TFQMRSolver;
+// typedef BITL_Solver< bicgstab_ell_type > B_BiCGStabEllSolver; // ERROR
+ typedef BITL_Solver< gmres_type > B_GMResSolver;
+// typedef BITL_Solver< idr_s_type > B_IDRsSolver; // ERROR
+ typedef BITL_Solver< minres_solver_type> B_MinResSolver;
+ typedef BITL_Solver< gcr_type > B_GcrSolver;
+ typedef BITL_Solver< fgmres_type > B_FGMResSolver;
+ typedef BITL_Solver< preonly_type > B_PreOnly;
+
+ // ===================================================================================
+
+ typedef ITL_Preconditioner<itl::pc::diagonal<BlockMTLMatrix>, BlockMTLMatrix, MTLTypes::MTLVector > BlockDiagonalPreconditioner;
+ typedef ITL_Preconditioner<itl::pc::identity<BlockMTLMatrix>, BlockMTLMatrix, MTLTypes::MTLVector > BlockIdentityPreconditioner;
+
+ template<typename MatrixType>
+ void initCreatorMap(std::string backend)
+ {
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<cg_solver_type, MatrixType, MTLTypes::MTLVector> > MatCgSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<cgs_solver_type, MatrixType, MTLTypes::MTLVector> > MatCgsSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<bicgstab_type, MatrixType, MTLTypes::MTLVector> > MatBicgstabSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<tfqmr_solver_type, MatrixType, MTLTypes::MTLVector> > MatTfqmrSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<gmres_type, MatrixType, MTLTypes::MTLVector> > MatGmresSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<minres_solver_type, MatrixType, MTLTypes::MTLVector> > MatMinresSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<gcr_type, MatrixType, MTLTypes::MTLVector> > MatGcrSolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<preonly_type, MatrixType, MTLTypes::MTLVector> > MatPreonlySolver;
+ typedef MTL4Solver<MatrixType, MTLTypes::MTLVector, ITL_Runner<fgmres_type, MatrixType, MTLTypes::MTLVector> > MatFgmresSolver;
+
+ CreatorMap<LinearSolver>::addCreator(backend + "_cg", new typename MatCgSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_cgs", new typename MatCgsSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_bicgstab", new typename MatBicgstabSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_tfqmr", new typename MatTfqmrSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_gmres", new typename MatGmresSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_minres", new typename MatMinresSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_gcr", new typename MatGcrSolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_preconly", new typename MatPreonlySolver::Creator);
+ CreatorMap<LinearSolver>::addCreator(backend + "_fgmres", new typename MatFgmresSolver::Creator);
+ }
+} // namespace AMDiS
+
+#endif // AMDIS_BITL_SOLVER
+
diff --git a/AMDiS/src/solver/BlockMTLMatrix.h b/AMDiS/src/solver/BlockMTLMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..78742d5aa4a3e0bf53f6cff578f8fc2204925945
--- /dev/null
+++ b/AMDiS/src/solver/BlockMTLMatrix.h
@@ -0,0 +1,180 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+/** \file BlockMTLMatrix.h */
+
+#ifndef AMDIS_BLOCK_MTL_MATRIX_H
+#define AMDIS_BLOCK_MTL_MATRIX_H
+
+#include "solver/SolverMatrix.h"
+#include "solver/Mapper.h"
+#include "MTL4Types.h"
+
+namespace AMDiS
+{
+
+ /// A wrapper for AMDiS::SolverMatrix to be used in MTL/ITL solvers
+ struct BlockMTLMatrix
+ {
+ typedef typename mtl::Collection<MTLTypes::MTLMatrix>::size_type size_type;
+
+ size_type n; // overall number of rows
+ size_type m; // overall number of columns
+
+ size_t n_rows; // number of row blocks
+ size_t n_cols; // number of column blocks
+
+ BlockMTLMatrix()
+ : A(NULL), initialized(false)
+ { }
+
+ BlockMTLMatrix(const SolverMatrix<Matrix<DOFMatrix*> >& A_)
+ : A(&A_), initialized(false)
+ {
+ init();
+ }
+
+ void setMatrix(const SolverMatrix<Matrix<DOFMatrix*> >& A_)
+ {
+ A = &A_;
+ init();
+ }
+
+ void init()
+ {
+ RectangularMapper mapper(*A);
+ n = mapper.getNumRows();
+ m = mapper.getNumCols();
+ n_rows = mapper.getNumRowComponents();
+ n_cols = mapper.getNumColComponents();
+
+ r_rows.resize(n_rows);
+ r_cols.resize(n_cols);
+ int start = 0;
+ mapper.setCol(0);
+ for (size_t i = 0; i < n_rows; i++) {
+ mapper.setRow(i+1);
+ int finish = mapper.row(0);
+ r_rows[i].set(start, finish);
+ start = finish;
+ }
+ start = 0;
+ mapper.setRow(0);
+ for (size_t i = 0; i < n_cols; i++) {
+ mapper.setCol(i+1);
+ int finish = mapper.col(0);
+ r_cols[i].set(start, finish);
+ start = finish;
+ }
+ initialized = true;
+ }
+
+ inline const mtl::irange& getRowRange(size_t i) const
+ {
+ return r_rows[i];
+ }
+
+ inline const mtl::irange& getColRange(size_t i) const
+ {
+ return r_cols[i];
+ }
+
+ const DOFMatrix::base_matrix_type& getSubMatrix(size_t i, size_t j) const
+ {
+ return A->getSubMatrix(i, j);
+ }
+
+ /// perform blockwise multiplication A*b -> x
+ template <typename VectorIn, typename VectorOut, typename Assign>
+ void mult(const VectorIn& b, VectorOut& x, Assign) const
+ {
+ TEST_EXIT_DBG(initialized)("Block matrix not initialized. Assign a SolverMatrix first!\n");
+
+ for (size_t i = 0; i < n_rows; i++) {
+ VectorOut x_i(x[r_rows[i]]);
+ bool first = true;
+ for (size_t j = 0; j < n_cols; j++) {
+ if ((*A->getOriginalMat())[i][j]) {
+ const VectorIn b_j(b[r_cols[j]]);
+ if (first) {
+ Assign::first_update(x_i, A->getSubMatrix(i, j) * b_j);
+ first = false;
+ } else {
+ Assign::update(x_i, A->getSubMatrix(i, j) * b_j);
+ }
+ }
+ }
+ }
+ }
+
+ template <typename VectorIn>
+ mtl::vector::mat_cvec_multiplier<BlockMTLMatrix, VectorIn> operator*(const VectorIn& v) const
+ { return mtl::vector::mat_cvec_multiplier<BlockMTLMatrix, VectorIn>(*this, v); }
+
+ protected:
+ bool initialized;
+
+ private:
+ const SolverMatrix<Matrix<DOFMatrix*> >* A;
+ std::vector<mtl::irange> r_rows, r_cols;
+ };
+
+ namespace dispatch
+ {
+ template< typename M >
+ void initMatrix(BlockMTLMatrix& m, MapperBase<M>& mapper)
+ { }
+
+ template< typename M >
+ void fillMatrix(BlockMTLMatrix& m, const SolverMatrix<Matrix<DOFMatrix*> >& source, MapperBase<M>& mapper)
+ {
+ m.setMatrix(source);
+ }
+
+ } // end namespace dispatch
+
+} // end namespace AMDiS
+
+
+inline AMDiS::BlockMTLMatrix::size_type size(const AMDiS::BlockMTLMatrix& A) { return (A.n) * (A.m); }
+inline AMDiS::BlockMTLMatrix::size_type num_rows(const AMDiS::BlockMTLMatrix& A) { return A.n; }
+inline AMDiS::BlockMTLMatrix::size_type num_cols(const AMDiS::BlockMTLMatrix& A) { return A.m; }
+
+namespace mtl
+{
+ template <>
+ struct Collection<AMDiS::BlockMTLMatrix>
+ {
+ typedef double value_type;
+ typedef int size_type;
+ };
+
+ namespace ashape
+ {
+ template <> struct ashape_aux<AMDiS::BlockMTLMatrix>
+ {
+ typedef nonscal type;
+ };
+ } // end namespace ashape
+
+} // end namespace mtl
+
+#endif // AMDIS_BLOCK_MTL_MATRIX_H
\ No newline at end of file
diff --git a/AMDiS/src/solver/BlockPreconditioner.h b/AMDiS/src/solver/BlockPreconditioner.h
new file mode 100644
index 0000000000000000000000000000000000000000..ec00bdda149b42bfde38e444afab42375bb056e7
--- /dev/null
+++ b/AMDiS/src/solver/BlockPreconditioner.h
@@ -0,0 +1,134 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+/** \file BlockPreconditioner.h */
+
+#ifndef AMDIS_BLOCK_PRECONDITIONER_H
+#define AMDIS_BLOCK_PRECONDITIONER_H
+
+#include "solver/ITL_Preconditioner.h"
+#include "solver/SolverMatrix.h"
+
+namespace AMDiS {
+
+ /// Basis preconditioner structure for block-preconditioners
+ template<typename MatrixType>
+ struct BlockPreconditioner : ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector>
+ {
+ typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
+ typedef ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector> super;
+ typedef BlockPreconditioner<MatrixType> self;
+
+ BlockPreconditioner()
+ : A(NULL), fullMatrix(NULL)
+ { }
+
+ virtual ~BlockPreconditioner() {}
+
+ /// extract iranges from BlockMatrix to be used to extract sub-vectors and sub-matrices.
+ void init(const SolverMatrix<Matrix<DOFMatrix*> >& A_, const MatrixType& fullMatrix_) override
+ {
+ A = &A_;
+ fullMatrix = &fullMatrix_;
+
+ BlockMapper mapper(A_);
+ rows.resize(mapper.getNumComponents());
+ int start = 0;
+ for (int i = 0; i < mapper.getNumComponents(); i++) {
+ mapper.setRow(i+1);
+ int finish = mapper.row(0);
+ rows[i].set(start, finish);
+ start = finish;
+ }
+ }
+
+ virtual void solve(const MTLTypes::MTLVector& b, MTLTypes::MTLVector& x) const
+ { FUNCNAME("BlockPreconditioner::solve()");
+ TEST_EXIT(false)("Must be implemented in derived classes!\n");
+ }
+
+ virtual void adjoint_solve(const MTLTypes::MTLVector& x, MTLTypes::MTLVector& y) const
+ { FUNCNAME("BlockPreconditioner::adjoint_solve()");
+ TEST_EXIT(false)("Must be implemented in derived classes!\n");
+ }
+
+ inline const mtl::irange& getRowRange(size_t i) const
+ {
+ return rows[i];
+ }
+
+ inline const mtl::irange& getColRange(size_t i) const
+ {
+ return rows[i];
+ }
+
+ const DOFMatrix::base_matrix_type& getSubMatrix(size_t i, size_t j) const
+ {
+ return A->getSubMatrix(i, j);
+ }
+
+ template<typename SolverType, typename RunnerType>
+ static void createSubSolver(std::string param, SolverType*& solver, RunnerType*& runner, std::string solverType = "0")
+ {
+ // definition of standard-backends
+#if defined HAVE_PARALLEL_PETSC
+ std::string backend("p_petsc");
+#elif defined HAVE_PARALLEL_MTL
+ std::string backend("p_mtl");
+#elif defined HAVE_PETSC
+ std::string backend("petsc");
+#else
+ std::string backend("mtl");
+#endif
+
+ // === read backend-name ===
+ std::string initFileStr = param + "->solver";
+ Parameters::get(initFileStr + "->backend", backend);
+
+ // === read solver-name ===
+ Parameters::get(initFileStr, solverType);
+
+ if (backend != "0" && backend != "no" && backend != "")
+ solverType = backend + "_" + solverType;
+
+ LinearSolverCreator *solverCreator =
+ dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolver>::getCreator(solverType, initFileStr));
+ TEST_EXIT(solverCreator)
+ ("No valid solver type found in parameter \"%s\"\n", initFileStr.c_str());
+ solverCreator->setName(initFileStr);
+ solver = dynamic_cast<SolverType*>(solverCreator->create());
+ assert(solver != NULL);
+
+ runner = dynamic_cast<RunnerType*>(solver->getRunner());
+ assert(runner != NULL);
+ }
+
+
+ protected:
+ const SolverMatrix<Matrix<DOFMatrix*> >* A;
+ const MatrixType* fullMatrix;
+
+ std::vector<mtl::irange> rows;
+ };
+
+} // end namespace AMDiS
+
+#endif // AMDIS_BLOCK_PRECONDITIONER_H
diff --git a/AMDiS/src/solver/CombinedPreconditioner.h b/AMDiS/src/solver/CombinedPreconditioner.h
new file mode 100644
index 0000000000000000000000000000000000000000..43e0ef7ba544f70f3dbfdb93b37a70c41eb34787
--- /dev/null
+++ b/AMDiS/src/solver/CombinedPreconditioner.h
@@ -0,0 +1,153 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+/** \file CombinedPreconditioner.h */
+
+#ifndef AMDIS_COMBINED_PRECONDITIONER_H
+#define AMDIS_COMBINED_PRECONDITIONER_H
+
+#include "solver/BlockPreconditioner.h"
+
+namespace AMDiS {
+
+ /// preconditioner structure that combines various block-preconditioners
+ template<typename MatrixType>
+ struct CombinedPreconditioner : BlockPreconditioner<MatrixType>
+ {
+ typedef BlockPreconditioner<MatrixType> super;
+ typedef ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector> precon_base;
+ typedef CombinedPreconditioner<MatrixType> self;
+
+ class Creator : public CreatorInterfaceName<precon_base>
+ {
+ public:
+ Creator(int l0)
+ { l.push_back(l0); };
+
+ Creator(int l0, int l1)
+ { l.push_back(l0); l.push_back(l1); };
+
+ Creator(int l0, int l1, int l2)
+ { l.push_back(l0); l.push_back(l1); l.push_back(l2); };
+
+ Creator(const std::vector<int>& l) : l(l) {}
+
+ virtual ~Creator() {}
+
+ precon_base* create() {
+ return new self(l);
+ }
+ private:
+ std::vector<int> l;
+ };
+
+ /// Constructor
+ CombinedPreconditioner(const std::vector<int>& parts_)
+ : super(), parts(parts_)
+ {
+ precon.resize(parts.size(), &identity);
+ subA.resize(parts.size());
+ }
+
+ virtual ~CombinedPreconditioner() {}
+
+ /// Extract iranges from SolverMatrix to be used to extract sub-vectors and sub-matrices.
+ void init(const SolverMatrix<Matrix<DOFMatrix*> >& A_, const MatrixType& fullMatrix_) override
+ {
+ super::A = &A_;
+ super::fullMatrix = &fullMatrix_;
+
+ BlockMapper mapper(A_);
+ rows.resize(mapper.getNumComponents());
+ int start = 0;
+ int sum = 0;
+ for (int p = 0; p < parts.size(); p++) {
+ sum += parts[p];
+ if (sum > mapper.getNumComponents())
+ throw std::runtime_error("range out of bound!");
+ mapper.setRow(sum);
+ int finish = mapper.row(0);
+ if (finish <= start)
+ throw std::runtime_error("ranges overlap!");
+ rows[p].set(start, finish);
+ start = finish;
+ }
+
+ const Matrix<DOFMatrix*>& mat = *A_.getOriginalMat();
+ sum = 0;
+ for (int p = 0; p < parts.size(); p++) {
+ Matrix<DOFMatrix*>* subMat = new Matrix<DOFMatrix*>(parts[p],parts[p]);
+ for (int i = 0; i < parts[p]; i++)
+ for (int j = 0; j < parts[p]; j++)
+ (*subMat)[i][j] = mat[i+sum][j+sum];
+ subA[p].setMatrix(*subMat);
+ precon[p]->init(subA[p], fullMatrix_);
+ }
+ }
+
+ void exit() override
+ {
+ for (int p = 0; p < parts.size(); p++) {
+ precon[p]->exit();
+ delete subA[p].getOriginalMat();
+ }
+ }
+
+
+ /// Apply the preconditioners block-wise
+ /**
+ * solve Px = b, with
+ * P = diag(P1, P2, ...)
+ **/
+ void solve(const MTLTypes::MTLVector& b, MTLTypes::MTLVector& x) const override
+ { FUNCNAME("CombinedPreconditioner::solve()");
+
+ x.change_dim(num_rows(b));
+
+ for (size_t i = 0; i < precon.size(); i++) {
+ const MTLTypes::MTLVector b_i(b[rows[i]]);
+ MTLTypes::MTLVector x_i(x[rows[i]]);
+ precon[i]->solve(b_i, x_i);
+ }
+ }
+
+ /// Sets a preconditioner for part i. If non is set, an identity preconditioner is used.
+ void setPreconditioner(size_t i, precon_base& p)
+ {
+ precon[i] = &p;
+ }
+
+ protected:
+ std::vector<mtl::irange> rows;
+
+ // length of blocks assigned to a seperate preconditioners
+ std::vector<int> parts;
+
+ std::vector<precon_base*> precon;
+ std::vector<SolverMatrix<Matrix<DOFMatrix*> > > subA;
+
+ DOFMatrix::base_matrix_type dummy;
+ ITL_Preconditioner<itl::pc::identity<MatrixType>, MatrixType, MTLTypes::MTLVector > identity;
+ };
+
+} // end namespace AMDiS
+
+#endif // AMDIS_COMBINED_PRECONDITIONER_H
diff --git a/AMDiS/src/solver/ITL_Preconditioner.h b/AMDiS/src/solver/ITL_Preconditioner.h
index ca048c12a116d27bb2989259bf039591e7df487a..8106043d1241acbae4f07e0cd0ebccb2731ca8f6 100644
--- a/AMDiS/src/solver/ITL_Preconditioner.h
+++ b/AMDiS/src/solver/ITL_Preconditioner.h
@@ -51,8 +51,7 @@ namespace AMDiS {
virtual ~ITL_BasePreconditioner() {}
- typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
- virtual void init(const BlockMatrix& A, const MatrixType& fullMatrix) = 0;
+ virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) = 0;
virtual void solve(const VectorType& x, VectorType& y) const = 0;
virtual void adjoint_solve(const VectorType& x, VectorType& y) const = 0;
};
diff --git a/AMDiS/src/solver/ITL_Runner.h b/AMDiS/src/solver/ITL_Runner.h
index aa1216da9542b93f32f5de56d2578d07f57452de..036198a39683be007ab6af4215cf4710c0eaeee4 100644
--- a/AMDiS/src/solver/ITL_Runner.h
+++ b/AMDiS/src/solver/ITL_Runner.h
@@ -82,7 +82,7 @@ namespace AMDiS {
/// Implementation of \ref MTL4Runner::init()
- void init(const typename super::BlockMatrix& A, const MatrixType& fullMatrix) override
+ void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
{
preconPair.l->init(A, fullMatrix);
preconPair.r->init(A, fullMatrix);
@@ -209,7 +209,7 @@ namespace AMDiS {
leftCreator = dynamic_cast<CreatorInterfaceName< ITL_BasePreconditioner<MatrixType, VectorType> >*>(
CreatorMap<ITL_BasePreconditioner<MatrixType, VectorType> >::getCreator(preconType, initFileStr) );
TEST_EXIT(leftCreator != nullptr)
- ("there is no creator for the given left preconditioner\n");
+ ("There is no creator for the given left preconditioner '%s'\n", preconType.c_str());
leftCreator->setName(initFileStr);
preconType = "no";
@@ -218,7 +218,7 @@ namespace AMDiS {
rightCreator = dynamic_cast<CreatorInterfaceName< ITL_BasePreconditioner<MatrixType, VectorType> >*>(
CreatorMap<ITL_BasePreconditioner<MatrixType, VectorType> >::getCreator(preconType, initFileStr) );
TEST_EXIT(rightCreator != nullptr)
- ("there is no creator for the given right preconditioner\n");
+ ("There is no creator for the given right preconditioner '%s'\n", preconType.c_str());
rightCreator->setName(initFileStr);
pair.l = leftCreator->create();
diff --git a/AMDiS/src/solver/KrylovPreconditioner.h b/AMDiS/src/solver/KrylovPreconditioner.h
index 8c573d2f87f13e90eb28de18935f2d2d7a7c704f..53bec1dd52ae0b3dffac545cc82db1c8735289bb 100644
--- a/AMDiS/src/solver/KrylovPreconditioner.h
+++ b/AMDiS/src/solver/KrylovPreconditioner.h
@@ -105,7 +105,7 @@ namespace AMDiS {
}
/// Implementation of \ref ITL_BasePreconditioner::init()
- void init(const typename super::BlockMatrix& A, const MatrixType& fullMatrix_) override
+ void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix_) override
{
fullMatrix = &fullMatrix_;
runner->init(A, fullMatrix_);
diff --git a/AMDiS/src/solver/LinearSolver.h b/AMDiS/src/solver/LinearSolver.h
index bfd1b1cae095095a11a5e9fd848fb4923760793e..e3f3e8924a884983ae5fafe7c52fe3d4937d54c1 100644
--- a/AMDiS/src/solver/LinearSolver.h
+++ b/AMDiS/src/solver/LinearSolver.h
@@ -123,9 +123,6 @@ namespace AMDiS {
rel_residual = -1.0;
int error_code = solveLinearSystem(A, x, b, createMatrixData, storeMatrixData);
- TEST_EXIT(error_code == 0)
- ("Error [%d] in LinearSolver. residual = %e, rel_residual = %e\n", error_code, residual, rel_residual);
-
// calculate and print resiual
if (info > 0) {
if (residual >= 0.0 && rel_residual >= 0.0) {
diff --git a/AMDiS/src/solver/MTL4Solver.h b/AMDiS/src/solver/MTL4Solver.h
index 0917b36cd016d39bda407ae46676b4b515f31944..518604b59a1cd02846f2b83888e55400e505d346 100644
--- a/AMDiS/src/solver/MTL4Solver.h
+++ b/AMDiS/src/solver/MTL4Solver.h
@@ -35,39 +35,80 @@ namespace AMDiS {
template< typename MatrixType, typename VectorType >
struct MTL4Runner : public OEMRunner
{
- typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
+// typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
+
+ virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) = 0;
- virtual void init(const BlockMatrix& A, const MatrixType& fullMatrix) = 0;
virtual int solve(const MatrixType& A, VectorType& x, const VectorType& b) = 0;
- virtual int adjoint_solve(const MatrixType& A, VectorType& x, const VectorType& b) {
+
+ virtual int adjoint_solve(const MatrixType& A, VectorType& x, const VectorType& b)
+ {
FUNCNAME("MTL4Runner::adjoint_solve()");
ERROR_EXIT("Must be implemented in derived class!\n");
return 0;
- };
+ }
};
-#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
- namespace Parallel {
- typedef LinearSolver ParallelSolver;
- }
- typedef BlockMapper ParallelMapper;
-#endif
-
- using Parallel::ParallelSolver;
-
/// Wrapper for template-argument dependent constructors
- template <bool parallel = false>
- struct MTL4SolverBase : LinearSolver
+ template < typename MatrixType, typename Enable = void>
+ struct MTL4SolverBase : public LinearSolver
{
- MTL4SolverBase(std::string name) : LinearSolver(name) {}
+ typedef BlockMapper Mapper;
+
+ MTL4SolverBase(std::string name)
+ : LinearSolver(name) {}
+
+ MatrixType& getMatrix()
+ {
+ return matrix;
+ }
+
+ /// create a sequential BlockMapper
+ void initMapper(const SolverMatrix<Matrix<DOFMatrix*> >& A)
+ {
+ mapper = new BlockMapper(A);
+ }
+
+ void exitMapper()
+ {
+ delete mapper;
+ }
+
+ protected:
+ MatrixType matrix;
+ Mapper* mapper;
};
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- template <>
- struct MTL4SolverBase<true> : ParallelSolver
+#ifdef HAVE_PARALLEL_MTL4
+ template < typename MatrixType >
+ struct MTL4SolverBase<MatrixType, typename boost::enable_if<mtl::traits::is_distributed<MatrixType> > >
+ : public ParallelSolver
{
- MTL4SolverBase(std::string name) : ParallelSolver(name, false) {}
+ typedef ParallelMapper Mapper;
+
+ MTL4SolverBase(std::string name)
+ : ParallelSolver(name, false) {}
+
+ MatrixType& getMatrix()
+ {
+ return matrix;
+ }
+
+ /// create a parallel mapper based on local-to-global mapping
+ void initMapper(const SolverMatrix<Matrix<DOFMatrix*> >& A)
+ {
+ mapper = new ParallelMapper(*ParallelSolver::getDofMapping());
+ }
+
+ void exitMapper()
+ {
+ delete mapper;
+ }
+
+ protected:
+ MatrixType matrix;
+ Mapper* mapper;
};
#endif
@@ -80,18 +121,15 @@ namespace AMDiS {
* solvers where MTL4 provides an interface, can be assigned
* by different Runner objects.
**/
- template < typename MatrixType, typename VectorType, typename Runner, bool parallel = false/*mtl::traits::is_distributed<MatrixType>::value*/ >
- class MTL4Solver : public MTL4SolverBase<parallel>
+ template < typename MatrixType, typename VectorType, typename Runner >
+ class MTL4Solver : public MTL4SolverBase<MatrixType>
{
private:
- typedef MTL4SolverBase<parallel> super;
- typedef MTL4Solver<MatrixType, VectorType, Runner> self;
-
- typedef typename boost::mpl::if_c<parallel, ParallelMapper, BlockMapper>::type Mapper;
+ typedef MTL4SolverBase<MatrixType> super;
+ typedef MTL4Solver<MatrixType, VectorType, Runner> self;
+ typedef typename super::Mapper Mapper;
Runner runner; // redirect the implementation to a runner
- MatrixType matrix;
- Mapper* mapper;
public:
/// Creator class used in the LinearSolverMap.
@@ -109,8 +147,8 @@ namespace AMDiS {
/// Constructor
MTL4Solver(std::string n)
- : super(n),
- runner(this)
+ : super(n),
+ runner(this)
{}
@@ -135,93 +173,63 @@ namespace AMDiS {
}
protected:
- /// create a sequential BlockMapper
- template<bool par>
- void initMapper(const SolverMatrix<Matrix<DOFMatrix*> >& A)
- {
- mapper = new BlockMapper(A);
- }
-
-#ifdef HAVE_PARALLEL_MTL4
- /// create a parallel mapper based on local-to-global mapping
- template<>
- void initMapper<true>(const SolverMatrix<Matrix<DOFMatrix*> >& A)
- {
- mapper = new ParallelMapper(*ParallelSolver::getDofMapping());
- }
-#endif
/// Implementation of \ref LinearSolver::solveLinearSystem()
int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
- SystemVector& x,
- SystemVector& b,
- bool createMatrixData,
- bool storeMatrixData) override
+ SystemVector& x,
+ SystemVector& b,
+ bool createMatrixData,
+ bool storeMatrixData) override
{
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
MPI::COMM_WORLD.Barrier();
#endif
- initMapper<parallel>(A);
+ super::initMapper(A);
Timer t;
if (createMatrixData) {
- initMatrix(matrix, A, *mapper);
- runner.init(A, matrix);
+ initMatrix(super::matrix, A, *super::mapper);
+ runner.init(A, super::matrix);
}
VectorType mtl_x;
- initVector(mtl_x, x, *mapper);
+ initVector(mtl_x, x, *super::mapper);
VectorType mtl_b;
- initVector(mtl_b, b, *mapper);
+ initVector(mtl_b, b, *super::mapper);
INFO(self::getInfo(), 8)("fill MTL4 matrix needed %.5f seconds\n", t.elapsed());
- int error = runner.solve(matrix ,mtl_x, mtl_b);
- VecMap< SystemVector, Mapper > xVecMap(x, *mapper);
+ int error = runner.solve(super::matrix ,mtl_x, mtl_b);
+
+ VecMap<SystemVector, Mapper> xVecMap(x, *super::mapper);
mtl_x >> xVecMap;
if (!storeMatrixData)
runner.exit();
- delete mapper;
+ super::exitMapper();
return error;
}
// functions to initialize mtl-matrix and mtl-vector
// from AMDiS matrix / vectors using mappers
- /// initialize a MTL matrix
- template< typename MatrixT, typename M >
- void initMatrix(MatrixT& m, MapperBase<M>& mapper)
- {
- dispatch::initMatrix(m, mapper.self(), typename traits::distributed_tag<MatrixT>::type() );
- }
-
/// initialize a MTL matrix and assign values from an AMDiS matrix
template< typename Matrix1, typename Matrix2, typename M >
void initMatrix(Matrix1& target, const Matrix2& source, MapperBase<M>& mapper)
{
- dispatch::initMatrix(target, mapper.self(), typename traits::distributed_tag<Matrix1>::type() );
- MatMap< const Matrix2, M > matMap(source, mapper.self());
- target << matMap;
- }
-
- /// initialize a MTL vector
- template< typename VectorT >
- void initVector(VectorT& v)
- {
- dispatch::initVector(v, matrix, typename traits::distributed_tag<VectorT>::type() );
- }
+ dispatch::initMatrix(target, mapper.self());
+ dispatch::fillMatrix(target, source, mapper.self());
+ }
/// initialize a MTL vector and assign values from an AMDiS vector
template< typename Vector1, typename Vector2, typename M >
void initVector(Vector1& target, const Vector2& source, MapperBase<M>& mapper)
{
- initVector(target);
- VecMap< const Vector2, M > srcVecMap(source, mapper.self());
- target << srcVecMap;
+ dispatch::initVector(target, super::matrix);
+ dispatch::fillVector(target, source, mapper.self());
}
};
}
diff --git a/AMDiS/src/solver/MTL4SolverBase.h b/AMDiS/src/solver/MTL4SolverBase.h
index 51616e584d8d48db042a38896b403826bc747259..198534537a83e73d7ea07151f8680b7ef9900c9b 100644
--- a/AMDiS/src/solver/MTL4SolverBase.h
+++ b/AMDiS/src/solver/MTL4SolverBase.h
@@ -54,21 +54,35 @@ namespace AMDiS {
namespace dispatch
{
+ /// init matrix
+ template< typename MatrixOut, typename M >
+ void initMatrix(MatrixOut& m, MapperBase<M>& mapper) {}
+
+ /// init vector
+ template< typename VectorOut, typename MatrixT >
+ void initVector(VectorOut& v, const MatrixT& source) {}
+
/// init systemmatrix depending on Mapper parameters.
- template< typename MatrixT, typename M >
- void initMatrix(MatrixT& m, MapperBase<M>& mapper, tag::non_distributed)
+ template< typename M >
+ void initMatrix(MTLTypes::MTLMatrix& m, MapperBase<M>& mapper)
{
m.change_dim(mapper.getNumRows(), mapper.getNumCols());
set_to_zero(m);
}
-
+
+ template< typename MatrixOut, typename MatrixIn, typename M >
+ void fillMatrix(MatrixOut& m, const MatrixIn& source, MapperBase<M>& mapper)
+ {
+ MatMap< const MatrixIn, M > matMap(source, mapper.self());
+ m << matMap;
+ }
#ifdef HAVE_PARALLEL_MTL4
/// init systemmatrix depending on Mapper parameters,
/// specialized for distributed matrices
- template< typename MatrixT >
- void initMatrix(MatrixT& m, ParallelMapper& mapper, tag::distributed)
+ template< typename M >
+ void initMatrix(MTLTypes::PMTLMatrix& m, MapperBase<M>& mapper)
{
mtl::par::block_distribution dist(mapper.getNumRows());
dist.setup_from_local_size(mapper.getMap().getLocalDofs());
@@ -76,27 +90,35 @@ namespace AMDiS {
m.init_distribution(dist, dist, mapper.getNumRows(), mapper.getNumRows());
set_to_zero(m);
}
+
+ /// init MTL-vector depending on Mapper parameters,
+ /// specialized for distributed matrices
+ template< typename MatrixT >
+ void initVector(MTLTypes::PMTLVector& v, const MatrixT& matrix)
+ {
+ v.init_distribution(row_distribution(matrix), num_rows(matrix));
+ set_to_zero(v);
+ }
#endif
-
/// init MTL-vector depending on Mapper parameters.
- template< typename VectorT, typename MatrixT >
- void initVector(VectorT& v, const MatrixT& matrix, tag::non_distributed)
+ template< typename MatrixT >
+ void initVector(MTLTypes::MTLVector& v, const MatrixT& matrix)
{
v.change_dim(num_rows(matrix));
set_to_zero(v);
}
-
- /// init MTL-vector depending on Mapper parameters,
- /// specialized for distributed matrices
- template< typename VectorT, typename MatrixT >
- void initVector(VectorT& v, const MatrixT& matrix, tag::distributed)
+ /// fill MTL-vector
+ template< typename VectorOut, typename VectorIn, typename M >
+ void fillVector(VectorOut& v, const VectorIn& source, MapperBase<M>& mapper)
{
- v.init_distribution(row_distribution(matrix), num_rows(matrix));
- set_to_zero(v);
+ VecMap< const VectorIn, M > srcVecMap(source, mapper.self());
+ v << srcVecMap;
}
+
+
} // end namespace dispatch
} // end namespace AMDiS
diff --git a/AMDiS/src/solver/Mapper.h b/AMDiS/src/solver/Mapper.h
index 22f1a7906d84589ddb78ab54e9928efa54a7a85b..960f6b7816bf4533785a2b8ea29c7c13552f5fda 100644
--- a/AMDiS/src/solver/Mapper.h
+++ b/AMDiS/src/solver/Mapper.h
@@ -61,6 +61,9 @@ namespace AMDiS {
/// return overall number of columns
inline size_type getNumCols() const { return self().getNumCols(); }
+
+ inline size_type getNumRows(unsigned int comp) const { return self().getNumRows(comp); }
+ inline size_type getNumCols(unsigned int comp) const { return self().getNumCols(comp); }
/// return number of components/blocks
inline unsigned int getNumComponents() const { return self().getNumComponents(); }
@@ -146,6 +149,9 @@ namespace AMDiS {
inline size_type getNumRows() const { return nrow; }
inline size_type getNumCols() const { return ncol; }
+ inline size_type getNumRows(unsigned int comp) const { return sizes[comp]; }
+ inline size_type getNumCols(unsigned int comp) const { return sizes[comp]; }
+
inline unsigned int getNumComponents() const { return nComp; }
private: // methods
@@ -169,6 +175,86 @@ namespace AMDiS {
std::vector< size_type > sizes;
};
+
+
+
+
+ /// Mapper implementation for non-parallel rectangular block matrices
+ struct RectangularMapper : public MapperBase< RectangularMapper >
+ {
+ typedef unsigned int size_type;
+
+ /// Constructor for block-matrices
+ RectangularMapper(const SolverMatrix<Matrix<DOFMatrix* > >& sm )
+ : nRowComp(sm.getOriginalMat()->getNumRows()),
+ nColComp(sm.getOriginalMat()->getNumCols()),
+ rowOffset(0), colOffset(0), nrow(0), ncol(0),
+ sizes_rows(nRowComp), sizes_cols(nColComp)
+ {
+ const Matrix<DOFMatrix* >& orMat(*sm.getOriginalMat());
+
+ for (int i= 0; i < nRowComp; i++) {
+ for (int j = 0; j < nColComp; j++) {
+ if (orMat[i][j]) {
+ sizes_rows[i] = orMat[i][j]->getRowFeSpace()->getAdmin()->getUsedSize();
+ sizes_cols[j] = orMat[i][j]->getColFeSpace()->getAdmin()->getUsedSize();
+ }
+ }
+ nrow += sizes_rows[i];
+ }
+ for (int j = 0; j < nColComp; j++)
+ ncol += sizes_cols[j];
+ }
+
+ /// calculate row offset for row component \param r
+ inline void setRow( unsigned int r)
+ {
+ assert( r <= sizes_rows.size() );
+ rowOffset = sum(r, sizes_rows);
+ }
+
+ /// calculate column offset for col component \param c
+ inline void setCol( unsigned int c)
+ {
+ assert( c <= sizes_cols.size() );
+ colOffset = sum(c, sizes_cols);
+ }
+
+ inline size_type row(size_type r) const { return r + rowOffset; }
+ inline size_type col(size_type c) const { return c + colOffset; }
+
+ inline size_type getNumRows() const { return nrow; }
+ inline size_type getNumCols() const { return ncol; }
+
+ inline size_type getNumRows(unsigned int comp) const { return sizes_rows[comp]; }
+ inline size_type getNumCols(unsigned int comp) const { return sizes_cols[comp]; }
+
+ inline unsigned int getNumComponents() const { return nRowComp; }
+ inline unsigned int getNumRowComponents() const { return nRowComp; }
+ inline unsigned int getNumColComponents() const { return nColComp; }
+
+ private: // methods
+
+ ///compute the sum of sizes from [0, end)
+ inline size_type sum(unsigned int end, std::vector< size_type >& sizes) const
+ {
+ unsigned int ret(0);
+ for (unsigned int i(0); i < end; ++i)
+ ret += sizes[i];
+ return ret;
+ }
+
+ private: // variables
+
+ unsigned int nRowComp, nColComp;
+ size_type rowOffset;
+ size_type colOffset;
+ unsigned int nrow;
+ unsigned int ncol;
+
+ std::vector< size_type > sizes_rows;
+ std::vector< size_type > sizes_cols;
+ };
}
#endif // AMDIS_MAPPER_H
diff --git a/AMDiS/src/solver/MatrixStreams.h b/AMDiS/src/solver/MatrixStreams.h
index 2d43dc870ca722720a2614bf5797c400973d16b9..1193da9f6ae5add2c75f67f8039814f199337356 100644
--- a/AMDiS/src/solver/MatrixStreams.h
+++ b/AMDiS/src/solver/MatrixStreams.h
@@ -30,13 +30,6 @@
#include "traits/inserter.h"
#include "traits/category.hpp"
-#ifdef HAVE_SEQ_PETSC
-#include "solver/PetscTypes.h"
-#include <petsc.h>
-#include <petscmat.h>
-#include <petscvec.h>
-#endif
-
namespace AMDiS {
using namespace mtl::operations;
@@ -181,268 +174,5 @@ namespace AMDiS {
dest >> swap;
}
}
-
-#ifdef HAVE_SEQ_PETSC
-
- /// fill PETSc matrix from DOFMatrix
- inline void operator<<(Mat& mat, const DOFMatrix& rhs)
- {
- bool initMatrix = false;
- if (mat == PETSC_NULL) {
- std::vector<PetscInt> nnz(rhs.getSize());
- for (size_t k = 0; k < static_cast<size_t>(rhs.getSize()); k++)
- nnz[k] = rhs.getBaseMatrix().nnz_local(k);
-
- MatCreateSeqAIJ(PETSC_COMM_SELF, num_rows(rhs.getBaseMatrix()), num_cols(rhs.getBaseMatrix()), 0, &(nnz[0]), &mat);
- initMatrix = true;
- }
-
- std::vector<PetscInt> indices;
- for (size_t i = 0; i < rhs.getBaseMatrix().ref_minor().size(); i++)
- indices.push_back(rhs.getBaseMatrix().ref_minor()[i]);
-
- for (PetscInt i = 0; i < static_cast<PetscInt>(num_rows(rhs.getBaseMatrix())); i++) {
- if (rhs.getBaseMatrix().nnz_local(i) > 0)
- MatSetValues(mat, 1, &i, rhs.getBaseMatrix().nnz_local(i),
- &(indices[rhs.getBaseMatrix().ref_major()[i]]),
- &(rhs.getBaseMatrix().value_from_offset(rhs.getBaseMatrix().ref_major()[i])), ADD_VALUES);
- }
-
- if (initMatrix) {
- MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
- }
- }
-
-
- /// create and fill \ref PetscMatrix
- inline void operator<<(PetscMatrix& mat, const SolverMatrix<Matrix<DOFMatrix*> >& Asolver)
- {
- if (mat.assembled)
- mat.destroy();
-
- const Matrix< DOFMatrix* >& A = *(Asolver.getOriginalMat());
- if (mat.isNested) {
-
- std::vector<PetscInt> nnz;
- mat.nestMat.resize(A.getNumRows() * A.getNumCols());
-
- for (size_t i = 0; i < static_cast<size_t>(A.getNumRows()); i++) {
- for (size_t j = 0; j < static_cast<size_t>(A.getNumCols()); j++) {
- size_t idx = i * A.getNumCols() + j;
-
- if (A[i][j] == nullptr
- || num_rows(A[i][j]->getBaseMatrix()) == 0
- || num_cols(A[i][j]->getBaseMatrix()) == 0
- || A[i][j]->getBaseMatrix().nnz() == 0) {
- mat.nestMat[idx] = PETSC_NULL;
- continue;
- }
-
- mat.nestMat[idx] << *(A[i][j]);
- }
- }
-
- // create nested matrix from a vector of block matrices
- MatCreateNest(PETSC_COMM_SELF, A.getNumRows(), PETSC_NULL, A.getNumCols(), PETSC_NULL, &(mat.nestMat[0]), &mat.matrix);
- } else {
- bool initMatrix = false;
- unsigned nRows(0);
- unsigned nEntries(0);
- std::vector< int > nRowsBlock(A.getNumRows());
- for(int i(0); i < A.getNumRows(); ++i) {
- int j(0);
- for( ; j<A.getNumCols() && A[i][j] == NULL; ++j) ;
- if( j == A.getNumCols() ) {
- std::stringstream ss; ss << "ERROR: solver matrix has empty row " << i << "\n";
- throw std::runtime_error(ss.str());
- }
- nRowsBlock[i]=num_rows(A[i][j]->getBaseMatrix());
- nRows+=nRowsBlock[i];
- }
- std::vector<PetscInt> nnz(nRows);
- //initialize the nnz (could be avoided, but gets complicated
- for (int i(0); i < nRows; ++i) nnz[i] = 0;
-
- //build a list of pairs (col, value) for each row
- std::vector< std::list< std::pair<int, double > > > rowData(nRows);
- //use the standard mapper to build global structure
- BlockMapper mapper(Asolver);
- PetscInt maxNnz(0);
- for (int i(0); i < A.getNumRows(); ++i) {
- mapper.setRow(i);
- //compute the number of nonzeros per global row
- for (int j(0); j < A.getNumCols(); ++j) {
- if (A[i][j] != NULL) {
- mapper.setCol(j);
- const DOFMatrix::base_matrix_type& mtlMat(A[i][j]->getBaseMatrix());
- for (unsigned k(0); k < nRowsBlock[i]; ++k) {
- unsigned curRow(mapper.row(k));
- nnz[curRow] += mtlMat.nnz_local(k);
- maxNnz = std::max(maxNnz, nnz[curRow]);
- unsigned curPos(mtlMat.ref_major()[k]);
- for(unsigned l(0); l < mtlMat.nnz_local(k); ++l, ++curPos) {
- rowData[curRow].push_back(std::make_pair(mapper.col(mtlMat.ref_minor()[curPos]), mtlMat.data[curPos]));
- }
- }
- }
- }
-
- }
- if (mat.matrix == PETSC_NULL) {
- MatCreateSeqAIJ(PETSC_COMM_SELF, nRows, nRows, 0, &(nnz[0]), &mat.matrix);
- initMatrix = true;
- }
- //reduce mallocs..
- std::vector< PetscInt > colIndices(maxNnz);
- std::vector< PetscScalar > values(maxNnz);
- std::list< std::pair< int, double > >::iterator it;
- unsigned j;
- for (PetscInt i(0); i < nRows; ++i) {
- j=0; it = rowData[i].begin();
- for (; it != rowData[i].end(); ++it, ++j) {
- colIndices[j] = it->first;
- values[j] = it->second;
- }
- MatSetValues(mat.matrix, 1, &i, nnz[i], &colIndices[0], &values[0], INSERT_VALUES);
- }
- if (initMatrix) {
- MatAssemblyBegin(mat.matrix, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(mat.matrix, MAT_FINAL_ASSEMBLY);
- }
- }
- mat.assembled = true;
- }
-
-
- /// fill PETSc vector from DOFVector
- inline void operator<<(Vec& petscVec, const DOFVector<double>& vec)
- {
- // Traverse all used DOFs in the dof vector.
- DOFConstIterator<double> dofIt(&vec, USED_DOFS);
- PetscInt index = 0;
- for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
- double value = *dofIt;
- VecSetValue(petscVec, index, value, ADD_VALUES);
- }
- }
-
-
- /// DOFVector from PETSc vector
- inline void operator>>(const Vec& petscVec, DOFVector<double>& vec)
- {
- // Traverse all used DOFs in the dof vector.
- DOFVector<double>::Iterator dofIt(&vec, USED_DOFS);
- PetscInt index = 0;
- for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
- double value = 0.0;
- VecGetValues(petscVec, 1, &index, &value);
- *dofIt = value;
- }
- }
-
-
- /// fill PETSc vector from DOFVector
- inline void operator<<(Vec& petscVec, const SystemVector& vec)
- {
-#ifdef VecType // PETSc uses MACROS instead of typedefs in Versions 3.3x
- const VecType vecType;
-#else
- VecType vecType;
-#endif
- VecGetType(petscVec, &vecType);
-
- if (strcmp(vecType, VECNEST) == 0) {
- for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
- Vec v;
- VecNestGetSubVec(petscVec, i, &v);
- v << *(vec.getDOFVector(i));
- }
- } else {
- PetscInt index = 0;
- for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
- DOFConstIterator<double> dofIt(vec.getDOFVector(i), USED_DOFS);
- for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
- double value = *dofIt;
- VecSetValue(petscVec, index, value, ADD_VALUES);
- }
- }
- }
- }
-
-
- /// fill SystemVector from PETSc vector
- inline void operator>>(const Vec& petscVec, SystemVector& vec)
- {
-#ifdef VecType // PETSc uses MACROS instead of typedefs in Versions 3.3x
- const VecType vecType;
-#else
- VecType vecType;
-#endif
- VecGetType(petscVec, &vecType);
-
- if (strcmp(vecType, VECNEST) == 0) {
- for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
- Vec v;
- VecNestGetSubVec(petscVec, i, &v);
- v >> *(vec.getDOFVector(i));
- }
- } else {
- PetscInt n = 0;
- for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++)
- n += vec.getDOFVector(i)->getUsedSize();
-
- PetscInt N = 0;
- VecGetSize(petscVec, &N);
- assert(n == N);
-
- PetscInt index = 0;
- for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
- DOFVector<double>::Iterator dofIt(vec.getDOFVector(i), USED_DOFS);
- for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
- double value = 0.0;
- VecGetValues(petscVec, 1, &index, &value);
- *dofIt = value;
- }
- }
- }
- }
-
-
- /// create and fill \ref PetscVector
- inline void operator<<(PetscVector& petscVec, const SystemVector& rhs)
- {
- size_t nComponents = rhs.getSize();
- if (petscVec.assembled)
- petscVec.destroy();
-
- // only create nested vector if requested, i.e. flag createNestedVector is set to true
- if (petscVec.isNested) {
- petscVec.nestVec.resize(nComponents);
-
- for (size_t i = 0; i < nComponents; i++) {
- VecCreateSeq(PETSC_COMM_SELF,rhs.getDOFVector(i)->getUsedSize(), &(petscVec.nestVec[i]));
-
- petscVec.nestVec[i] << *(rhs.getDOFVector(i));
-
- VecAssemblyBegin(petscVec.nestVec[i]);
- VecAssemblyEnd(petscVec.nestVec[i]);
- }
-
- // create nested vector from vector of block vectors
- VecCreateNest(PETSC_COMM_SELF, nComponents, PETSC_NULL, &(petscVec.nestVec[0]), &(petscVec.vector));
- } else {
- PetscInt n = 0;
- for (size_t i = 0; i < nComponents; i++)
- n += rhs.getDOFVector(i)->getUsedSize();
-
- VecCreateSeq(PETSC_COMM_SELF, n, &(petscVec.vector));
- petscVec.vector << rhs;
- }
- petscVec.assembled = true;
- }
-
-#endif
-
}
#endif // AMDIS_MATRIXSTREAMS_H
diff --git a/AMDiS/src/solver/PetscSolver.cc b/AMDiS/src/solver/PetscSolver.cc
index 256d885ebbcbe5be527696a97fee2dc9928ac5de..710cb81754175a008d4377fb5383642039c148e8 100644
--- a/AMDiS/src/solver/PetscSolver.cc
+++ b/AMDiS/src/solver/PetscSolver.cc
@@ -22,141 +22,14 @@
#ifdef HAVE_SEQ_PETSC
#include "solver/PetscSolver.h"
-
-using namespace std;
+#include "solver/PetscTypes.h"
namespace AMDiS {
- PetscRunner::PetscRunner(LinearSolver* oem_)
- : oem(*oem_),
- kspPrefix(""),
- zeroStartVector(false),
- initialized(false)
- {
- PetscParameters params;
- bool matSolverPackage = false;
-
- // set the parameter prefix
- Parameters::get(oem.getName() + "->petsc prefix", kspPrefix); // important when more than one solver to configure
-
- // set the solver
- std::string solverName = "petsc";
- Parameters::get(oem.getName(), solverName);
- if (solverName == "petsc")
- Parameters::get(oem.getName() + "->ksp_type", solverName);
-
- std::string kspSolver = params.solverMap[solverName];
-
- if (params.matSolverPackage[kspSolver]) {
- // direct solvers
- PetscOptionsInsertString(("-" + kspPrefix + "ksp_type preonly").c_str());
- PetscOptionsInsertString(("-" + kspPrefix + "pc_type lu").c_str());
- PetscOptionsInsertString(("-" + kspPrefix + "pc_factor_mat_solver_package " + (kspSolver != "direct" ? kspSolver : "umfpack")).c_str());
- oem.setMaxIterations(1);
- zeroStartVector = true;
- matSolverPackage = true;
- } else if (!params.emptyParam[kspSolver]) {
- // other solvers
- PetscOptionsInsertString(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
- }
-
- // set the preconditioner
- string precon = "";
- Parameters::get(oem.getName() + "->pc_type", precon);
- if (!precon.size())
- Parameters::get(oem.getName() + "->left precon", precon);
- if (!precon.size())
- Parameters::get(oem.getName() + "->right precon", precon);
- if (!matSolverPackage && !params.emptyParam[precon]) {
- precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
- PetscOptionsInsertString(("-" + kspPrefix + "pc_type " + precon).c_str());
- }
-
- if (oem.getInfo() >= 20)
- PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor_true_residual").c_str());
- else if (oem.getInfo() >= 10)
- PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor").c_str());
-
- // command line string
- string kspString = "";
- Parameters::get(oem.getName() + "->ksp", kspString);
- if (kspString != "")
- PetscOptionsInsertString(kspString.c_str());
- }
-
-
- void PetscRunner::init(const BlockMatrix& A, const Mat& fullMatrix)
- {
- if (initialized)
- exit();
-
- createSubSolver(ksp, fullMatrix, kspPrefix);
-
- // set tolerance options from LinearSolver-parameters
- KSPSetTolerances(ksp, oem.getRelative(), oem.getTolerance(), PETSC_DEFAULT, oem.getMaxIterations());
-
- // Do not delete the solution vector, use it for the initial guess.
- if (!zeroStartVector)
- KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
-
- KSPGetPC(ksp, &pc);
- PCSetFromOptions(pc);
-
- initialized = true;
- }
-
+ // explicit template instantiation
+ template class PetscRunner<PetscMatrix, PetscVector>;
+ template class PetscRunner<PetscMatrixNested, PetscVectorNested>;
- int PetscRunner::solve(const Mat& A, Vec& x, const Vec& b)
- {
- if (zeroStartVector)
- VecSet(x, 0.0);
-
- PetscErrorCode solverError = KSPSolve(ksp, b, x);
-
- PetscInt nIter = 0;
- KSPGetIterationNumber(ksp, &nIter);
- PetscReal residual_norm = -1.0;
- KSPGetResidualNorm(ksp, &residual_norm);
-
- if (residual_norm <= 0.0) {
- Vec r;
- VecDuplicate(b, &r);
- KSPBuildResidual(ksp, PETSC_NULL, PETSC_NULL, &r);
- VecNorm(r, NORM_2, &residual_norm);
- }
-
- oem.setErrorCode(solverError);
- oem.setIterations(nIter);
- oem.setResidual(residual_norm);
-
- return solverError;
- }
-
-
- void PetscRunner::createSubSolver(KSP &ksp_, Mat m, string kspPrefix_)
- {
- KSPCreate(PETSC_COMM_SELF, &ksp_);
- KSPSetOperators(ksp_, m, m, SAME_NONZERO_PATTERN);
- KSPSetOptionsPrefix(ksp_, kspPrefix_.c_str());
- KSPSetFromOptions(ksp_);
- }
-
-
- void PetscRunner::setSolver(KSP ksp_, std::string kspPrefix_,
- KSPType kspType, PCType pcType,
- PetscReal rtol, PetscReal atol, PetscInt maxIt)
- {
- KSPSetType(ksp_, kspType);
- KSPSetTolerances(ksp_, rtol, atol, PETSC_DEFAULT, maxIt);
- if (kspPrefix_ != "")
- KSPSetOptionsPrefix(ksp_, kspPrefix_.c_str());
- KSPSetFromOptions(ksp_);
-
- PC pc_;
- KSPGetPC(ksp_, &pc_);
- PCSetType(pc_, pcType);
- PCSetFromOptions(pc_);
- }
}
#endif // HAVE_SEQ_PETSC
diff --git a/AMDiS/src/solver/PetscSolver.h b/AMDiS/src/solver/PetscSolver.h
index 1ee3e5e0cfda218ba42f92834dc83ca82de5638e..c750c484ef6ac3864d58d799abc6987e4886a6b6 100644
--- a/AMDiS/src/solver/PetscSolver.h
+++ b/AMDiS/src/solver/PetscSolver.h
@@ -26,9 +26,9 @@
#ifdef HAVE_SEQ_PETSC
-#include "solver/LinearSolver.h"
-#include "solver/MatrixStreams.h"
+#include "solver/MTL4Solver.h"
#include "solver/PetscTypes.h"
+#include "solver/MatrixStreams.h"
#include "Timer.h"
#include <vector>
#include <iostream>
@@ -42,42 +42,90 @@
#include <petscao.h>
namespace AMDiS {
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief Common base class for wrappers to use Petsc preconditioners in AMDiS.
+ */
+ template< class MatrixType, class VectorType >
+ struct Petsc_BasePreconditioner : public OEMPreconditioner
+ {
+ Petsc_BasePreconditioner(std::string prefix, std::string name) : prefix(prefix), name(name) {}
+ virtual ~Petsc_BasePreconditioner() {}
+
+ virtual void init(PC pc, const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix)
+ {
+ PetscOptionsInsertString(("-" + prefix + "pc_type " + name).c_str());
+ PCSetFromOptions(pc);
+ }
+
+ virtual void exit() {}
+
+ protected:
+ std::string prefix;
+ std::string name;
+ };
+
+
+ typedef Petsc_BasePreconditioner<PetscMatrix, PetscVector> PetscPreconditioner;
+ typedef Petsc_BasePreconditioner<PetscMatrixNested, PetscVectorNested> PetscPreconditionerNested;
+
+ template<typename MatrixType, typename VectorType>
class PetscRunner : public OEMRunner
{
public:
/// Constructor of standard PETSc runner. Reads ksp and pc parameters from initfile.
- PetscRunner(LinearSolver* oem_);
-
- typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
+ PetscRunner(LinearSolver* oem);
+
+ /// Destructor, deletes preconditioner object \ref preconditioner
+ ~PetscRunner()
+ {
+ if (preconditioner) {
+ delete preconditioner;
+ preconditioner = NULL;
+ }
+ }
- /// initialize the solver \ref ksp and preconditioner \ref pc
- void init(const BlockMatrix& A, const Mat& fullMatrix);
+ /// Initialize the solver \ref ksp and preconditioner \ref pc
+ void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix);
- /// solve the linear equation \f$ A\cdot x = b \f$ by applying the PETSc solver \ref ksp
- int solve(const Mat& A, Vec& x, const Vec& b);
+ /// Solve the linear equation \f$ A\cdot x = b \f$ by applying the PETSc solver \ref ksp
+ int solve(const MatrixType& A, VectorType& x, const VectorType& b);
- /// destroy solver \ref ksp and preconditioner \ref pc
+ /// Destroy solver \ref ksp and preconditioner \ref pc
virtual void exit()
{
+ preconditioner->exit();
KSPDestroy(&ksp);
}
- /// get the PETSc solver \ref ksp
+ /// Get the PETSc solver \ref ksp
KSP getSolver()
{
return ksp;
}
- /// get the PETSc preconditioner \ref pc
+ /// Get the PETSc preconditioner \ref pc
PC getPc()
{
return pc;
}
- ~PetscRunner() { }
+ /// Returns preconditioner object \ref preconditioner
+ OEMPreconditioner* getLeftPrecon()
+ {
+ return preconditioner;
+ }
+
+ /// Returns preconditioner object \ref preconditioner
+ OEMPreconditioner* getRightPrecon()
+ {
+ return preconditioner;
+ }
+
- protected:
static void createSubSolver(KSP &ksp_, Mat m, std::string kspPrefix_);
static void setSolver(KSP ksp_, std::string kspPrefix_,
@@ -85,10 +133,15 @@ namespace AMDiS {
PetscReal rtol = PETSC_DEFAULT,
PetscReal atol = PETSC_DEFAULT,
PetscInt maxIt = PETSC_DEFAULT);
-
+
+
+ protected:
+ /// Creates a preconditioner object \ref preconditioner
+ void setPrecon();
+
+ protected:
LinearSolver& oem;
- //private:
/// PETSc solver object
KSP ksp;
@@ -98,8 +151,12 @@ namespace AMDiS {
/// KSP database prefix
std::string kspPrefix;
- bool zeroStartVector;
- bool initialized;
+ bool zeroStartVector; ///< initialize startsolution for iterative solver with 0
+ bool initialized; ///< true when \ref init() was called
+ bool matSolverPackage; ///< true when using a direct solver
+
+ /// Preconditioner object, that set parameters to \ref pc
+ Petsc_BasePreconditioner<MatrixType, VectorType>* preconditioner;
};
@@ -113,100 +170,13 @@ namespace AMDiS {
* This is a suite of data structures and routines for the
* scalable (parallel) solution of scientific applications.
*/
- template< typename Runner = PetscRunner >
- class PetscSolver : public LinearSolver
- {
- public:
- /// Creator class used in the LinearSolverMap.
- class Creator : public LinearSolverCreator
- {
- public:
- virtual ~Creator() {}
-
- /// Returns a new PetscSolver object.
- LinearSolver* create()
- {
- return new PetscSolver<Runner>(this->name);
- }
- };
-
- PetscSolver(std::string n)
- : LinearSolver(n),
- runner(this)
- {}
-
- ~PetscSolver() {}
-
- virtual OEMRunner* getRunner()
- {
- return &runner;
- }
-
- void setNestedVectors(bool nested = true)
- {
- vecSol.isNested = nested;
- vecRhs.isNested = nested;
- }
-
- void setNestedMatrix(bool nested = true)
- {
- petscMat.isNested=nested;
- }
-
- void setNested(bool n)
- {
- setNestedVectors(n);
- setNestedMatrix(n);
- }
-
- protected:
- Runner runner;
-
- int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
- SystemVector& x,
- SystemVector& b,
- bool createMatrixData,
- bool storeMatrixData)
- { FUNCNAME("PetscSolver::solveLinearSystem()");
-
- Timer t;
- // transfer matrix and rhs-vector to PETSc data-structures
- if (createMatrixData) {
- petscMat << A;
- runner.init(A, petscMat.matrix);
- }
-
- vecSol << x;
- vecRhs << b;
- INFO(info, 8)("fill PETSc matrix needed %.5f seconds\n", t.elapsed());
-
- // solve the linear system using PETSc solvers
- t.reset();
- error = runner.solve(petscMat.matrix, vecSol.vector, vecRhs.vector);
-
- // transfer solution from PETSc vector to SystemVector
- vecSol.vector >> x;
-
- vecSol.destroy();
- vecRhs.destroy();
-
- if (!storeMatrixData) {
- petscMat.destroy();
- runner.exit();
- }
- return error;
- }
-
- private:
- /// PETSc System-Matrix
- PetscMatrix petscMat;
-
- /// Solution and RHS vectors.
- PetscVector vecSol, vecRhs;
- };
+ typedef MTL4Solver< PetscMatrix, PetscVector, PetscRunner<PetscMatrix, PetscVector> > PetscSolver;
+ typedef MTL4Solver< PetscMatrixNested, PetscVectorNested, PetscRunner<PetscMatrixNested, PetscVectorNested> > PetscSolverNested;
} // end namespace AMDiS
#endif // HAVE_SEQ_PETSC
+#include "solver/PetscSolver.hh"
+
#endif // AMDIS_PETSC_SOLVER_SEQ_H
diff --git a/AMDiS/src/solver/PetscSolver.hh b/AMDiS/src/solver/PetscSolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..096178c7f3b4fb3d3732b19b8936fab9d665095c
--- /dev/null
+++ b/AMDiS/src/solver/PetscSolver.hh
@@ -0,0 +1,205 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+#ifdef HAVE_SEQ_PETSC
+
+namespace AMDiS {
+
+ template<typename M, typename V>
+ PetscRunner<M,V>::PetscRunner(LinearSolver* oem_)
+ : oem(*oem_),
+ kspPrefix(""),
+ zeroStartVector(false),
+ initialized(false),
+ preconditioner(NULL)
+ {
+ PetscParameters params;
+ matSolverPackage = false;
+
+ // set the parameter prefix
+ Parameters::get(oem.getName() + "->petsc prefix", kspPrefix); // important when more than one solver to configure
+
+ // set the solver
+ std::string solverName = "petsc";
+ Parameters::get(oem.getName(), solverName);
+ if (solverName == "petsc")
+ Parameters::get(oem.getName() + "->ksp_type", solverName);
+
+ std::string kspSolver = params.solverMap[solverName];
+
+ if (params.matSolverPackage[kspSolver]) {
+ // direct solvers
+ PetscOptionsInsertString(("-" + kspPrefix + "ksp_type preonly").c_str());
+ PetscOptionsInsertString(("-" + kspPrefix + "pc_type lu").c_str());
+ PetscOptionsInsertString(("-" + kspPrefix + "pc_factor_mat_solver_package " + (kspSolver != "direct" ? kspSolver : "umfpack")).c_str());
+ oem.setMaxIterations(1);
+ zeroStartVector = true;
+ matSolverPackage = true;
+ } else if (!params.emptyParam[kspSolver]) {
+ // other solvers
+ PetscOptionsInsertString(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
+ }
+
+ // set the preconditioner
+ setPrecon();
+
+ if (oem.getInfo() >= 20)
+ PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor_true_residual").c_str());
+ else if (oem.getInfo() >= 10)
+ PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor").c_str());
+
+ // command line string
+ std::string kspString = "";
+ Parameters::get(oem.getName() + "->ksp", kspString);
+ if (kspString != "")
+ PetscOptionsInsertString(kspString.c_str());
+ }
+
+
+ template<typename MatrixType, typename VectorType>
+ void PetscRunner<MatrixType, VectorType>::init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix)
+ {
+ if (initialized)
+ exit();
+
+ createSubSolver(ksp, fullMatrix.matrix, kspPrefix);
+
+ // set tolerance options from LinearSolver-parameters
+ KSPSetTolerances(ksp, oem.getRelative(), oem.getTolerance(), PETSC_DEFAULT, oem.getMaxIterations());
+
+ // Do not delete the solution vector, use it for the initial guess.
+ if (!zeroStartVector)
+ KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
+
+ KSPGetPC(ksp, &pc);
+
+ preconditioner->init(pc, A, fullMatrix);
+ initialized = true;
+ }
+
+
+ template<typename MatrixType, typename VectorType>
+ int PetscRunner<MatrixType, VectorType>::solve(const MatrixType& A, VectorType& x, const VectorType& b)
+ {
+ if (zeroStartVector)
+ VecSet(x.vector, 0.0);
+
+ PetscErrorCode solverError = KSPSolve(ksp, b.vector, x.vector);
+
+ PetscInt nIter = 0;
+ KSPGetIterationNumber(ksp, &nIter);
+ PetscReal residual_norm = -1.0;
+ KSPGetResidualNorm(ksp, &residual_norm);
+
+ if (residual_norm <= 0.0) {
+ Vec r;
+ VecDuplicate(b.vector, &r);
+ KSPBuildResidual(ksp, PETSC_NULL, PETSC_NULL, &r);
+ VecNorm(r, NORM_2, &residual_norm);
+ }
+
+ oem.setErrorCode(solverError);
+ oem.setIterations(nIter);
+ oem.setResidual(residual_norm);
+
+ return solverError;
+ }
+
+
+ template<typename M, typename V>
+ void PetscRunner<M,V>::createSubSolver(KSP &ksp_, Mat m, std::string kspPrefix_)
+ {
+ KSPCreate(PETSC_COMM_SELF, &ksp_);
+#if (PETSC_VERSION_MINOR >= 5)
+ KSPSetOperators(ksp_, m, m);
+#else
+ KSPSetOperators(ksp_, m, m, SAME_NONZERO_PATTERN);
+#endif
+ KSPSetOptionsPrefix(ksp_, kspPrefix_.c_str());
+ KSPSetFromOptions(ksp_);
+ }
+
+
+ template<typename M, typename V>
+ void PetscRunner<M,V>::setSolver(KSP ksp_, std::string kspPrefix_,
+ KSPType kspType, PCType pcType,
+ PetscReal rtol, PetscReal atol, PetscInt maxIt)
+ {
+ KSPSetType(ksp_, kspType);
+ KSPSetTolerances(ksp_, rtol, atol, PETSC_DEFAULT, maxIt);
+ if (kspPrefix_ != "")
+ KSPSetOptionsPrefix(ksp_, kspPrefix_.c_str());
+ KSPSetFromOptions(ksp_);
+
+ PC pc_;
+ KSPGetPC(ksp_, &pc_);
+ PCSetType(pc_, pcType);
+ PCSetFromOptions(pc_);
+ }
+
+
+ template<typename MatrixType, typename VectorType>
+ void PetscRunner<MatrixType, VectorType>::setPrecon()
+ {
+ std::string precon = "";
+ std::string initFileStr = oem.getName() + "->pc_type";
+ Parameters::get(initFileStr, precon);
+ if (!precon.size() || precon == "no" || precon == "0") {
+ initFileStr = oem.getName() + "->left precon";
+ Parameters::get(initFileStr, precon);
+ }
+ if (!precon.size() || precon == "no" || precon == "0") {
+ initFileStr = oem.getName() + "->right precon";
+ Parameters::get(initFileStr, precon);
+ }
+ // no preconditioner should be created
+ if (!precon.size() || precon == "no" || precon == "0")
+ return;
+
+ PetscParameters params;
+ if (preconditioner) {
+ delete preconditioner;
+ preconditioner = NULL;
+ }
+
+ // Creator for the left preconditioner
+ CreatorInterfaceName< Petsc_BasePreconditioner<MatrixType, VectorType> >* creator(NULL);
+
+ try {
+ creator = dynamic_cast<CreatorInterfaceName< Petsc_BasePreconditioner<MatrixType, VectorType> >*>(
+ CreatorMap<Petsc_BasePreconditioner<MatrixType, VectorType> >::getCreator(precon, initFileStr) );
+ } catch (...) {
+ creator = NULL;
+ }
+
+ if (creator) {
+ creator->setName(initFileStr);
+ preconditioner = creator->create();
+ } else if (!matSolverPackage && !params.emptyParam[precon]) {
+ precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
+ preconditioner = new Petsc_BasePreconditioner<MatrixType, VectorType>(kspPrefix, precon);
+ } else {
+ ERROR_EXIT((std::string("There is no creator for the given preconditioner '") + precon + "'").c_str());
+ }
+ }
+}
+
+#endif // HAVE_SEQ_PETSC
diff --git a/AMDiS/src/solver/PetscTypes.h b/AMDiS/src/solver/PetscTypes.h
index 8802c4195d4a55fa592513cb4fe1954698f43f1d..504017af6c4b77ba3934a80ee238659be23452e8 100644
--- a/AMDiS/src/solver/PetscTypes.h
+++ b/AMDiS/src/solver/PetscTypes.h
@@ -23,11 +23,14 @@
#define AMDIS_PETSCTYPES_H
#include <vector>
+#include <map>
+#include <string>
+
+#include "MatrixStreams.h"
+
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
-#include <map>
-#include <string>
namespace AMDiS {
@@ -38,6 +41,8 @@ namespace AMDiS {
: matrix(PETSC_NULL),
isNested(false),
assembled(false) {}
+
+ ~PetscMatrix() { destroy(); }
std::vector<Mat> nestMat;
Mat matrix;
@@ -47,6 +52,14 @@ namespace AMDiS {
bool assembled;
};
+ struct PetscMatrixNested : PetscMatrix
+ {
+ PetscMatrixNested() : PetscMatrix()
+ {
+ isNested = true;
+ }
+ };
+
/// pair of nested vectors and blockvector
struct PetscVector
@@ -55,6 +68,8 @@ namespace AMDiS {
: vector(PETSC_NULL),
isNested(false),
assembled(false) {}
+
+ ~PetscVector() { destroy(); }
std::vector<Vec> nestVec;
Vec vector;
@@ -64,6 +79,13 @@ namespace AMDiS {
bool assembled;
};
+ struct PetscVectorNested : PetscVector
+ {
+ PetscVectorNested() : PetscVector()
+ {
+ isNested = true;
+ }
+ };
/// structure for parameter mapping
struct PetscParameters
@@ -76,6 +98,53 @@ namespace AMDiS {
PetscParameters();
};
-}
+
+ // matrix-streams for Petsc Types
+ // ___________________________________________________________________________
+
+ /// fill Petsc Mat datastructure directly
+ inline void operator<<(Mat& mat, const DOFMatrix& rhs);
+
+ /// fill PetscMatrix
+ template< typename Mapper >
+ inline void operator<<(PetscMatrix& mat, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& rhs);
+
+ /// fill nested PetscMatrix
+ template< typename Mapper >
+ inline void operator<<(PetscMatrixNested& mat, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& mapper);
+
+
+ /// fill Petsc Vec datastructure directly from DOFVector
+ inline void operator<<(Vec& petscVec, const DOFVector<double>& vec);
+
+ /// fill Petsc Vec datastructure directly from SystemVector
+ inline void operator<<(Vec& petscVec, const SystemVector& vec);
+
+ inline void operator<<(PetscVector& petscVec, const SystemVector& rhs);
+
+ template< typename Mapper >
+ inline void operator<<(PetscVector& petscVec, VecMap<const SystemVector, Mapper>& rhs);
+
+ template< typename Mapper >
+ inline void operator<<(PetscVectorNested& petscVec, VecMap<const SystemVector, Mapper>& rhs);
+
+
+ /// fill AMDiS DOFVector from Petsc Vec datastructur
+ inline void operator>>(const Vec& petscVec, DOFVector<double>& vec);
+
+ /// fill AMDiS SystemVector from Petsc Vec datastructur
+ inline void operator>>(const Vec& petscVec, SystemVector& vec);
+
+ /// fill AMDiS SystemVector from PetscVector
+ template< typename Mapper >
+ void operator>>(const PetscVector& dest, VecMap<SystemVector, Mapper>& rhs);
+
+ /// fill AMDiS SystemVector from nested PetscVector
+ template< typename Mapper >
+ void operator>>(const PetscVectorNested& dest, VecMap<SystemVector, Mapper>& rhs);
+
+} // end namespace AMDiS
+
+#include "solver/PetscTypes.hh"
#endif // AMDIS_PETSCTYPES_H
diff --git a/AMDiS/src/solver/PetscTypes.hh b/AMDiS/src/solver/PetscTypes.hh
new file mode 100644
index 0000000000000000000000000000000000000000..774650909683305d6c09f1056f274489825fca18
--- /dev/null
+++ b/AMDiS/src/solver/PetscTypes.hh
@@ -0,0 +1,295 @@
+#ifdef HAVE_SEQ_PETSC
+
+namespace AMDiS {
+
+ inline void operator<<(Mat& mat, const DOFMatrix& rhs)
+ {
+ bool initMatrix = false;
+ if (mat == PETSC_NULL) {
+ std::vector<PetscInt> nnz(rhs.getSize());
+ for (size_t k = 0; k < static_cast<size_t>(rhs.getSize()); k++)
+ nnz[k] = rhs.getBaseMatrix().nnz_local(k);
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF, num_rows(rhs.getBaseMatrix()), num_cols(rhs.getBaseMatrix()), 0, &(nnz[0]), &mat);
+ initMatrix = true;
+ }
+
+ std::vector<PetscInt> indices;
+ for (size_t i = 0; i < rhs.getBaseMatrix().ref_minor().size(); i++)
+ indices.push_back(rhs.getBaseMatrix().ref_minor()[i]);
+
+ for (PetscInt i = 0; i < static_cast<PetscInt>(num_rows(rhs.getBaseMatrix())); i++) {
+ if (rhs.getBaseMatrix().nnz_local(i) > 0)
+ MatSetValues(mat, 1, &i, rhs.getBaseMatrix().nnz_local(i),
+ &(indices[rhs.getBaseMatrix().ref_major()[i]]),
+ &(rhs.getBaseMatrix().value_from_offset(rhs.getBaseMatrix().ref_major()[i])), ADD_VALUES);
+ }
+
+ if (initMatrix) {
+ MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
+ }
+ }
+
+
+ template< typename Mapper >
+ inline void operator<<(PetscMatrix& mat, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& rhs)
+ {
+ if (mat.assembled)
+ mat.destroy();
+
+ const Matrix< DOFMatrix* >& A = *(rhs.mat.getOriginalMat());
+
+ bool initMatrix = false;
+ unsigned nRows(0);
+ unsigned nEntries(0);
+ std::vector< int > nRowsBlock(A.getNumRows());
+ for(int i(0); i < A.getNumRows(); ++i) {
+ int j(0);
+ for( ; j<A.getNumCols() && A[i][j] == NULL; ++j) ;
+ if( j == A.getNumCols() ) {
+ std::stringstream ss; ss << "ERROR: solver matrix has empty row " << i << "\n";
+ throw std::runtime_error(ss.str());
+ }
+ nRowsBlock[i]=num_rows(A[i][j]->getBaseMatrix());
+ nRows+=nRowsBlock[i];
+ }
+ std::vector<PetscInt> nnz(nRows);
+
+ //initialize the nnz (could be avoided, but gets complicated)
+ for (int i(0); i < nRows; ++i)
+ nnz[i] = 0;
+
+ //build a list of pairs (col, value) for each row
+ std::vector< std::list< std::pair<int, double > > > rowData(nRows);
+
+ PetscInt maxNnz(0);
+ Mapper &mapper = rhs.mapper;
+ for (int i(0); i < A.getNumRows(); ++i) {
+ mapper.setRow(i);
+ //compute the number of nonzeros per global row
+ for (int j(0); j < A.getNumCols(); ++j) {
+ if (A[i][j] != NULL) {
+ mapper.setCol(j);
+ const DOFMatrix::base_matrix_type& mtlMat(A[i][j]->getBaseMatrix());
+ for (unsigned k(0); k < nRowsBlock[i]; ++k) {
+ unsigned curRow(mapper.row(k));
+ nnz[curRow] += mtlMat.nnz_local(k);
+ maxNnz = std::max(maxNnz, nnz[curRow]);
+ unsigned curPos(mtlMat.ref_major()[k]);
+ for(unsigned l(0); l < mtlMat.nnz_local(k); ++l, ++curPos) {
+ rowData[curRow].push_back(std::make_pair(mapper.col(mtlMat.ref_minor()[curPos]), mtlMat.data[curPos]));
+ }
+ }
+ }
+ }
+
+ }
+ if (mat.matrix == PETSC_NULL) {
+ MatCreateSeqAIJ(PETSC_COMM_SELF, nRows, nRows, 0, &(nnz[0]), &mat.matrix);
+ initMatrix = true;
+ }
+ //reduce mallocs..
+ std::vector< PetscInt > colIndices(maxNnz);
+ std::vector< PetscScalar > values(maxNnz);
+ std::list< std::pair< int, double > >::iterator it;
+ unsigned j;
+ for (PetscInt i(0); i < nRows; ++i) {
+ j=0; it = rowData[i].begin();
+ for (; it != rowData[i].end(); ++it, ++j) {
+ colIndices[j] = it->first;
+ values[j] = it->second;
+ }
+ MatSetValues(mat.matrix, 1, &i, nnz[i], &colIndices[0], &values[0], INSERT_VALUES);
+ }
+ if (initMatrix) {
+ MatAssemblyBegin(mat.matrix, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat.matrix, MAT_FINAL_ASSEMBLY);
+ }
+
+ mat.assembled = true;
+ }
+
+
+ template< typename Mapper >
+ inline void operator<<(PetscMatrixNested& mat, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& mapper)
+ {
+ if (mat.assembled)
+ mat.destroy();
+
+ const Matrix< DOFMatrix* >& A = *(mapper.mat.getOriginalMat());
+
+ std::vector<PetscInt> nnz;
+ mat.nestMat.resize(A.getNumRows() * A.getNumCols());
+
+ for (unsigned int i = 0; i < A.getNumRows(); i++) {
+ for (unsigned int j = 0; j < A.getNumCols(); j++) {
+ size_t idx = i * A.getNumCols() + j;
+
+ if (!(A[i][j])
+ || num_rows(A[i][j]->getBaseMatrix()) == 0
+ || num_cols(A[i][j]->getBaseMatrix()) == 0
+ || A[i][j]->getBaseMatrix().nnz() == 0) {
+ mat.nestMat[idx] = PETSC_NULL;
+ continue;
+ }
+
+ mat.nestMat[idx] << *(A[i][j]);
+ }
+ }
+
+ // create nested matrix from a vector of block matrices
+ MatCreateNest(PETSC_COMM_SELF, A.getNumRows(), PETSC_NULL, A.getNumCols(), PETSC_NULL, &(mat.nestMat[0]), &mat.matrix);
+ mat.assembled = true;
+ }
+
+
+ inline void operator<<(Vec& petscVec, const DOFVector<double>& vec)
+ {
+ // Traverse all used DOFs in the dof vector.
+ DOFConstIterator<double> dofIt(&vec, USED_DOFS);
+ PetscInt index = 0;
+ for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
+ double value = *dofIt;
+ VecSetValue(petscVec, index, value, ADD_VALUES);
+ }
+ }
+
+
+ inline void operator>>(const Vec& petscVec, DOFVector<double>& vec)
+ {
+ // Traverse all used DOFs in the dof vector.
+ DOFVector<double>::Iterator dofIt(&vec, USED_DOFS);
+ PetscInt index = 0;
+ for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
+ double value = 0.0;
+ VecGetValues(petscVec, 1, &index, &value);
+ *dofIt = value;
+ }
+ }
+
+
+ inline void operator<<(Vec& petscVec, const SystemVector& vec)
+ {
+#ifdef VecType // PETSc uses MACROS instead of typedefs in Versions 3.3x
+ const VecType vecType;
+#else
+ VecType vecType;
+#endif
+ VecGetType(petscVec, &vecType);
+
+ if (strcmp(vecType, VECNEST) == 0) {
+ for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
+ Vec v;
+ VecNestGetSubVec(petscVec, i, &v);
+ v << *(vec.getDOFVector(i));
+ }
+ } else {
+ PetscInt index = 0;
+ for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
+ DOFConstIterator<double> dofIt(vec.getDOFVector(i), USED_DOFS);
+ for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
+ double value = *dofIt;
+ VecSetValue(petscVec, index, value, ADD_VALUES);
+ }
+ }
+ }
+ }
+
+
+ inline void operator>>(const Vec& petscVec, SystemVector& vec)
+ {
+#ifdef VecType // PETSc uses MACROS instead of typedefs in Versions 3.3x
+ const VecType vecType;
+#else
+ VecType vecType;
+#endif
+ VecGetType(petscVec, &vecType);
+
+ if (strcmp(vecType, VECNEST) == 0) {
+ for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
+ Vec v;
+ VecNestGetSubVec(petscVec, i, &v);
+ v >> *(vec.getDOFVector(i));
+ }
+ } else {
+ PetscInt n = 0;
+ for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++)
+ n += vec.getDOFVector(i)->getUsedSize();
+
+ PetscInt N = 0;
+ VecGetSize(petscVec, &N);
+ assert(n == N);
+
+ PetscInt index = 0;
+ for (size_t i = 0; i < static_cast<size_t>(vec.getSize()); i++) {
+ DOFVector<double>::Iterator dofIt(vec.getDOFVector(i), USED_DOFS);
+ for (dofIt.reset(); !dofIt.end(); ++dofIt, ++index) {
+ double value = 0.0;
+ VecGetValues(petscVec, 1, &index, &value);
+ *dofIt = value;
+ }
+ }
+ }
+ }
+
+
+ template< typename Mapper >
+ void operator>>(const PetscVector& dest, VecMap<SystemVector, Mapper>& rhs)
+ {
+ dest.vector >> rhs.vec;
+ }
+
+
+ template< typename Mapper >
+ void operator>>(const PetscVectorNested& dest, VecMap<SystemVector, Mapper>& rhs)
+ {
+ dest.vector >> rhs.vec;
+ }
+
+
+ template< typename Mapper >
+ inline void operator<<(PetscVector& petscVec, VecMap<const SystemVector, Mapper>& rhs)
+ {
+ size_t nComponents = rhs.vec.getSize();
+ if (petscVec.assembled)
+ petscVec.destroy();
+
+ PetscInt n = 0;
+ for (size_t i = 0; i < nComponents; i++)
+ n += rhs.vec.getDOFVector(i)->getUsedSize();
+
+ VecCreateSeq(PETSC_COMM_SELF, n, &(petscVec.vector));
+ petscVec.vector << rhs.vec;
+
+ petscVec.assembled = true;
+ }
+
+
+ template< typename Mapper >
+ inline void operator<<(PetscVectorNested& petscVec, VecMap<const SystemVector, Mapper>& rhs)
+ {
+ size_t nComponents = rhs.mapper.getNumComponents();
+ if (petscVec.assembled)
+ petscVec.destroy();
+
+ petscVec.nestVec.resize(nComponents);
+
+ for (size_t i = 0; i < nComponents; i++) {
+ VecCreateSeq(PETSC_COMM_SELF,rhs.mapper.getNumRows(i), &(petscVec.nestVec[i]));
+
+ petscVec.nestVec[i] << *(rhs.vec.getDOFVector(i));
+
+ VecAssemblyBegin(petscVec.nestVec[i]);
+ VecAssemblyEnd(petscVec.nestVec[i]);
+ }
+
+ // create nested vector from vector of block vectors
+ VecCreateNest(PETSC_COMM_SELF, nComponents, PETSC_NULL, &(petscVec.nestVec[0]), &(petscVec.vector));
+
+ petscVec.assembled = true;
+ }
+
+} // end namespace AMDiS
+
+#endif
diff --git a/AMDiS/src/solver/TriangularPreconditioner.h b/AMDiS/src/solver/TriangularPreconditioner.h
new file mode 100644
index 0000000000000000000000000000000000000000..b79dae83413211bc96c39f4e9d23008045140e1d
--- /dev/null
+++ b/AMDiS/src/solver/TriangularPreconditioner.h
@@ -0,0 +1,104 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+/** \file TriangularPreconditioner.h */
+
+#ifndef AMDIS_TRIANGULAR_PRECONDITIONER_H
+#define AMDIS_TRIANGULAR_PRECONDITIONER_H
+
+#include "solver/CombinedPreconditioner.h"
+
+namespace AMDiS {
+
+ /// preconditioner structure that combines various block-preconditioners
+ template<typename MatrixType>
+ struct TriangularPreconditioner : CombinedPreconditioner<MatrixType>
+ {
+ typedef CombinedPreconditioner<MatrixType> super;
+ typedef ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector> precon_base;
+ typedef TriangularPreconditioner<MatrixType> self;
+
+ class Creator : public CreatorInterfaceName<precon_base>
+ {
+ public:
+ Creator(int l0)
+ { l.push_back(l0); };
+
+ Creator(int l0, int l1)
+ { l.push_back(l0); l.push_back(l1); };
+
+ Creator(int l0, int l1, int l2)
+ { l.push_back(l0); l.push_back(l1); l.push_back(l2); };
+
+ Creator(const std::vector<int>& l) : l(l) {}
+
+ virtual ~Creator() {}
+
+ precon_base* create() {
+ return new self(l);
+ }
+ private:
+ std::vector<int> l;
+ };
+
+ /// Constructor
+ TriangularPreconditioner(const std::vector<int>& parts_)
+ : super(parts_)
+ { }
+
+ virtual ~TriangularPreconditioner() {}
+
+ /// Apply the preconditioners block-wise
+ /**
+ * solve Px = b, with
+ * P = diag(P1, P2, ...)
+ **/
+ void solve(const MTLTypes::MTLVector& b, MTLTypes::MTLVector& x) const override
+ { FUNCNAME("CombinedPreconditioner::solve()");
+
+ x.change_dim(num_rows(b));
+ y.change_dim(num_rows(b));
+ c.change_dim(num_rows(b));
+ set_to_zero(c);
+
+ for (size_t i = super::precon.size()-1; i >= 0; i--) {
+ MTLTypes::MTLVector y_i(b[super::rows[i]]);
+ if (i < super::precon.size()-1)
+ y_i -= y[super::rows[i]];
+
+ MTLTypes::MTLVector x_i(x[super::rows[i]]);
+ super::precon[i]->solve(y_i, x_i);
+
+ if (i > 0) {
+ c[super::rows[i]] = x_i;
+ y = (*super::fullMatrix) * c;
+ }
+ }
+ }
+
+ private:
+ mutable MTLTypes::MTLVector y;
+ mutable MTLTypes::MTLVector c;
+ };
+
+} // end namespace AMDiS
+
+#endif // AMDIS_COMBINED_PRECONDITIONER_H
diff --git a/AMDiS/src/solver/UmfPackSolver.h b/AMDiS/src/solver/UmfPackSolver.h
index 9c6818455d558a22cced06b7986c5738522a5611..81c8784d076ae5247c53d38ce384414764759daf 100644
--- a/AMDiS/src/solver/UmfPackSolver.h
+++ b/AMDiS/src/solver/UmfPackSolver.h
@@ -51,7 +51,7 @@ namespace AMDiS {
}
/// Implementation of \ref MTL4Runner::init()
- void init(const typename super::BlockMatrix& A, const MatrixType& fullMatrix) override
+ void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override
{
if (solver != nullptr) {
delete solver;
diff --git a/AMDiS/src/solver/itl/block_diagonal.hpp b/AMDiS/src/solver/itl/block_diagonal.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..71d8b92873328e08f96494c53082ebafd351ca83
--- /dev/null
+++ b/AMDiS/src/solver/itl/block_diagonal.hpp
@@ -0,0 +1,77 @@
+#ifndef ITL_PC_BLOCK_DIAGONAL_INCLUDE
+#define ITL_PC_BLOCK_DIAGONAL_INCLUDE
+
+#include <boost/numeric/itl/pc/diagonal.hpp>
+
+#include "solver/BlockMTLMatrix.h"
+
+namespace itl { namespace pc {
+
+/// Diagonal Preconditioner
+template <>
+class diagonal<AMDiS::BlockMTLMatrix, mtl::Collection<AMDiS::BlockMTLMatrix>::value_type>
+{
+ public:
+ typedef mtl::Collection<AMDiS::BlockMTLMatrix>::value_type Value;
+ typedef Value value_type;
+ typedef mtl::Collection<AMDiS::BlockMTLMatrix>::size_type size_type;
+ typedef diagonal self;
+
+ /// Constructor takes matrix reference
+ explicit diagonal(const AMDiS::BlockMTLMatrix& A) : inv_diag(num_rows(A))
+ {
+ MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
+ using math::reciprocal;
+
+ for (size_t i = 0; i < A.n_rows; i++)
+ inv_diag[A.getRowRange(i)] = mtl::matrix::diagonal(A.getSubMatrix(i,i));
+
+ for (size_type i= 0; i < num_rows(A); ++i)
+ inv_diag[i]= reciprocal(inv_diag[i]);
+ }
+
+ /// Member function solve, better use free function solve
+ template <typename Vector>
+ Vector solve(const Vector& x) const
+ {
+ Vector y(resource(x));
+ solve(x, y);
+ return y;
+ }
+
+ template <typename VectorIn, typename VectorOut>
+ void solve(const VectorIn& x, VectorOut& y) const
+ {
+ y.checked_change_resource(x);
+ MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
+ for (size_type i= 0; i < size(inv_diag); ++i)
+ y[i]= inv_diag[i] * x[i];
+ }
+
+ /// Member function for solving adjoint problem, better use free function adjoint_solve
+ template <typename Vector>
+ Vector adjoint_solve(const Vector& x) const
+ {
+ Vector y(resource(x));
+ adjoint_solve(x, y);
+ return y;
+ }
+
+ template <typename VectorIn, typename VectorOut>
+ void adjoint_solve(const VectorIn& x, VectorOut& y) const
+ {
+ using mtl::conj;
+ y.checked_change_resource(x);
+ MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
+ for (size_type i= 0; i < size(inv_diag); ++i)
+ y[i]= conj(inv_diag[i]) * x[i];
+ }
+
+ protected:
+ mtl::vector::dense_vector<value_type> inv_diag;
+};
+
+
+}} // namespace itl::pc
+
+#endif // ITL_PC_DIAGONAL_INCLUDE
diff --git a/AMDiS/src/traits/category.hpp b/AMDiS/src/traits/category.hpp
index e52059e986abb76d13dbc640d1ef55c1af08d6dc..be27b0a00de9e667a1d834c2ed16f4be944487df 100644
--- a/AMDiS/src/traits/category.hpp
+++ b/AMDiS/src/traits/category.hpp
@@ -30,6 +30,7 @@
#include "types.hpp"
#include "DOFMatrix.h"
+#include "solver/BlockMTLMatrix.h"
#include <boost/mpl/if.hpp>
#include <boost/mpl/or.hpp>
@@ -206,6 +207,14 @@ namespace AMDiS
typedef AMDiS::DOFMatrix::size_type size_type;
};
+ template <>
+ struct category< AMDiS::BlockMTLMatrix >
+ {
+ typedef tag::matrix tag;
+ typedef AMDiS::DOFMatrix::value_type value_type;
+ typedef AMDiS::DOFMatrix::size_type size_type;
+ };
+
template < class Matrix >
struct category< AMDiS::SolverMatrix<Matrix> >
{
diff --git a/AMDiS/src/traits/num_cols.hpp b/AMDiS/src/traits/num_cols.hpp
index 08f44137280cdff8ab04eb0e32691f5003caaa15..178de24c89acd2de797b1a19d0648db87191221f 100644
--- a/AMDiS/src/traits/num_cols.hpp
+++ b/AMDiS/src/traits/num_cols.hpp
@@ -139,6 +139,14 @@ namespace AMDiS
typedef category< AMDiS::DOFMatrix >::size_type type;
type operator()(const AMDiS::DOFMatrix& v) const { return mtl::matrix::num_cols(v.getBaseMatrix()); }
};
+
+ /// size implementation for AMDiS::BlockMTLMatrix
+ template <>
+ struct num_cols< AMDiS::BlockMTLMatrix >
+ {
+ typedef category< AMDiS::BlockMTLMatrix >::size_type type;
+ type operator()(const AMDiS::BlockMTLMatrix& M) const { return M.m; }
+ };
// === multi-vector or multi-multi-vector ===
diff --git a/AMDiS/src/traits/num_rows.hpp b/AMDiS/src/traits/num_rows.hpp
index e1372215957bddec352da24ad09d91d88113f1c8..e244582b87ba756fff04169f9d5be2a70344082f 100644
--- a/AMDiS/src/traits/num_rows.hpp
+++ b/AMDiS/src/traits/num_rows.hpp
@@ -139,6 +139,14 @@ namespace AMDiS
typedef category< AMDiS::DOFMatrix >::size_type type;
type operator()(const AMDiS::DOFMatrix& v) const { return mtl::matrix::num_rows(v.getBaseMatrix()); }
};
+
+ /// size implementation for AMDiS::BlockMTLMatrix
+ template <>
+ struct num_rows< AMDiS::BlockMTLMatrix >
+ {
+ typedef category< AMDiS::BlockMTLMatrix >::size_type type;
+ type operator()(const AMDiS::BlockMTLMatrix& M) const { return M.n; }
+ };
// === multi-vector or multi-multi-vector ===
diff --git a/AMDiS/src/traits/size.hpp b/AMDiS/src/traits/size.hpp
index 503bc6f7b03be7f7e8136d94b180f074900e6447..7d6e1b8148986b1548f79c9b4c46760edcd18abe 100644
--- a/AMDiS/src/traits/size.hpp
+++ b/AMDiS/src/traits/size.hpp
@@ -139,6 +139,14 @@ namespace AMDiS
typedef category< AMDiS::DOFMatrix >::size_type type;
type operator()(const AMDiS::DOFMatrix& v) const { return mtl::matrix::size(v.getBaseMatrix()); }
};
+
+ /// size implementation for AMDiS::BlockMTLMatrix
+ template <>
+ struct size< AMDiS::BlockMTLMatrix >
+ {
+ typedef category< AMDiS::BlockMTLMatrix >::size_type type;
+ type operator()(const AMDiS::BlockMTLMatrix& M) const { return (M.n) * (M.m); }
+ };
// === multi-vector or multi-matrix ===
diff --git a/AMDiS/src/utility/calculate_nnz.hpp b/AMDiS/src/utility/calculate_nnz.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d456739cc4931517bc8d0db447bcd93c014c479
--- /dev/null
+++ b/AMDiS/src/utility/calculate_nnz.hpp
@@ -0,0 +1,64 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+
+/** \file calculate_nnz.hpp */
+
+#ifndef AMDIS_UTILITY_CALCULATE_NNZ_HPP
+#define AMDIS_UTILITY_CALCULATE_NNZ_HPP
+
+namespace AMDiS {
+
+ inline void calculate_nnz(const DOFMatrix& matrix, std::vector<size_t>& nnz_per_row)
+ {
+ std::vector<std::set<DegreeOfFreefom> > nnz(nnz_per_row.size());
+ std::vector<DegreeOfFreedom> rowIndices, colIndices;
+ const FiniteElemSpace* rowFeSpace = matrix.getRowFeSpace(),
+ colFeSpace = matrix.getColFeSpace();
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(rowFeSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ rowFeSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
+ rowFeSpace->getAdmin(),
+ rowIndices);
+ if (rowFeSpace == colFeSpace) {
+ colIndices = rowIndices;
+ } else {
+ colFeSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
+ colFeSpace->getAdmin(),
+ colIndices);
+ }
+
+ for (size_t r = 0; r < rowIndices.size(); r++)
+ for (size_t c = 0; c < rowIndices.size(); c++)
+ nnz[rowIndices[r]].insert(colIndices[c]);
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ for (size_t r = 0; r < nnz.size(); r++)
+ nnz_per_row[r] = nnz[r].size();
+ }
+
+} // end namspace AMDiS
+
+#endif
\ No newline at end of file
diff --git a/AMDiS/src/utility/is_symmetric.hpp b/AMDiS/src/utility/is_symmetric.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f030ccf85be33b9206c2f5e30ee6cd52a7eccc8
--- /dev/null
+++ b/AMDiS/src/utility/is_symmetric.hpp
@@ -0,0 +1,53 @@
+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+
+
+/** \file is_symmetric.hpp */
+
+#ifndef AMDIS_UTILITY_IS_SYMMETRIC_HPP
+#define AMDIS_UTILITY_IS_SYMMETRIC_HPP
+
+namespace AMDiS {
+
+ template< class Matrix >
+ inline bool is_symmetric(const Matrix& matrix, double tol = 1.e-5)
+ {
+ using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
+ namespace traits= mtl::traits;
+
+ traits::row<Matrix>::type row(matrix);
+ traits::col<Matrix>::type col(matrix);
+ traits::const_value<Matrix>::type value(matrix);
+
+ typedef traits::range_generator<major, Matrix>::type cursor_type;
+ typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+ for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
+ for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
+ // Compare each non-zero entry with its transposed
+ if (abs(value(*icursor) - matrix[col(*icursor)][row(*icursor)]) > tol)
+ return false;
+ return true;
+ }
+
+} // end namspace AMDiS
+
+#endif
\ No newline at end of file