diff --git a/AMDiS/AMDiSConfig.cmake.in b/AMDiS/AMDiSConfig.cmake.in
index 878394e9c8cb74c52ec8b2c17a8025a64eed1475..113c0b97ff99401ae2ee9b2945568bfd3371fd8a 100644
--- a/AMDiS/AMDiSConfig.cmake.in
+++ b/AMDiS/AMDiSConfig.cmake.in
@@ -61,6 +61,9 @@ set(AMDiS_NEED_UMFPACK @ENABLE_UMFPACK@)
set(AMDiS_NEED_MKL @ENABLE_MKL@)
set(AMDIS_USE_FILE ${AMDIS_DIR}/AMDiSUse.cmake)
set(AMDiS_COMPILEFLAGS "@COMPILEFLAGS@")
+set(AMDIS_VERSION @CurrentRevision@)
+set(AMDIS_MAJOR_VERSION @AMDIS_MAJOR@)
+set(AMDIS_MINOR_VERSION @AMDIS_MINOR@)
if(AMDiS_NEED_UMFPACK)
#look for umfpack, we need the headers only
diff --git a/AMDiS/AMDiSUse.cmake b/AMDiS/AMDiSUse.cmake
index f1b3c4733be4261064d19732a30f538be6b9731e..d5d664b7c91174f4ef51b10751a194cf2d21fd1a 100644
--- a/AMDiS/AMDiSUse.cmake
+++ b/AMDiS/AMDiSUse.cmake
@@ -12,13 +12,19 @@ if(AMDiS_HAS_PARALLEL_DOMAIN)
list(APPEND AMDiS_COMPILEFLAGS ${MPI_COMPILE_FLAGS})
list(APPEND AMDIS_INCLUDE_DIRS ${MPI_INCLUDE_PATH})
endif(MPI_FOUND)
- set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${AMDIS_DIR})
- set(PETSC_EXECUTABLE_RUNS ON)
- find_package(PETSc REQUIRED)
- if(PETSC_FOUND)
+ if(AMDiS_HAS_PARALLEL_DOMAIN STREQUAL "PETSC")
+ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${AMDIS_DIR})
+ set(PETSC_EXECUTABLE_RUNS ON)
+ find_package(PETSc REQUIRED)
+ if(PETSC_FOUND)
list(APPEND AMDIS_LIBRARIES ${PETSC_LIBRARY_SYS} ${PETSC_LIBRARIES})
list(APPEND AMDIS_INCLUDE_DIRS ${PETSC_INCLUDES})
- endif(PETSC_FOUND)
+ endif(PETSc_FOUND)
+ elseif(AMDiS_HAS_PARALLEL_DOMAIN STREQUAL "PMTL")
+ find_package(MTL REQUIRED)
+ list(APPEND AMDIS_LIBRARIES ${MTL_LIBRARIES})
+# set(AMDiS_COMPILEFLAGS "${AMDIS_COMPILEFLAGS} -DMTL_HAS_MPI -DHAVE_PARALLEL_MTL4")
+ endif()
endif(NOT AMDiS_NEED_PARMETIS)
endif(AMDiS_HAS_PARALLEL_DOMAIN)
#thats bad because it affects each target
diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index df04e9336a146a73c75f43edd7178a671c6f34a4..bbf3e5f17fac14e448f99cc7d34098f5c99c061a 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -1,17 +1,31 @@
project(AMDiS)
cmake_minimum_required(VERSION 2.6)
-
-#needs: --enable-debug, enable-intel --enable-openmp --enable-parmetis --enable-parallel-domain --enable-umfpack --enable-mkl --enable-boost --enable-marmot
+#if revision is not set, let cmake look for the svn-version
+if( NOT CurrentRevision )
+ find_package(Subversion)
+ if(Subversion_FOUND)
+ Subversion_WC_INFO(${AMDiS_SOURCE_DIR} mtlSubinfo)
+ set(CurrentRevision "0.9${mtlSubinfo_WC_REVISION}")
+ # message("current revision: ${mtlSubinfo_WC_REVISION}")
+ endif(Subversion_FOUND)
+endif( NOT CurrentRevision )
+#split the current revision in MAJOR.MINOR
+string(REGEX MATCH "([0-9]+).([0-9]+)" AMDIS_MAJOR "${CurrentRevision}" )
+set(AMDIS_MINOR ${CMAKE_MATCH_2})
+set(AMDIS_MAJOR ${CMAKE_MATCH_1})
+message("major: ${AMDIS_MAJOR}")
+message("minor: ${AMDIS_MINOR}")
SET(LIB_DIR ${AMDiS_SOURCE_DIR}/lib)
SET(SOURCE_DIR ${AMDiS_SOURCE_DIR}/src)
#TODO: use the cmake build type
-SET(MTL_DIR ${LIB_DIR}/mtl4/ CACHE PATH "mtl4 directory")
+SET(MTL_INCLUDE_DIR ${LIB_DIR}/mtl4/ CACHE PATH "mtl4 directory")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -g -Wall -DDEBUG=0")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -DDEBUG=1 -O0")
+set(COMPILEFLAGS "-DAMDIS_VERSION=${CurrentRevision}")
#define the build type, empty can be everything and nothing
if(CMAKE_BUILD_TYPE STREQUAL "")
@@ -19,13 +33,12 @@ if(CMAKE_BUILD_TYPE STREQUAL "")
endif()
#option(ENABLE_INTEL "use intel compiler" false)
-option(ENABLE_PARALLEL_DOMAIN "use parallel domain decomposition" false)
+SET(ENABLE_PARALLEL_DOMAIN "OFF" CACHE STRING "use parallel domain decomposition. please set to one of: PMTL, PETSC, OFF" )
option(USE_PETSC_DEV false)
#parmetis is not optional
set(ENABLE_PARMETIS false)
option(ENABLE_ZOLTAN false)
option(ENABLE_UMFPACK "use umfpack solver" false)
-option(ENABLE_TESTS "compile the tests" false)
if(ENABLE_INTEL)
Message("please set the icc manually")
@@ -37,14 +50,16 @@ endif(ENABLE_INTEL)
find_package(Boost 1.42 REQUIRED)
if(Boost_FOUND)
- include_directories(${Boost_INCLUDE_DIR})
+# include_directories(${Boost_INCLUDE_DIR})
+ list(APPEND AMDiS_INCLUDE_DIRS ${Boost_INCLUDE_DIR})
message("boost lib-dirs: ${Boost_LIBRARY_DIRS}")
message("use the boost dir: ${Boost_INCLUDE_DIR}")
if(WIN32)
message("the windows find_boost does not set the boost library paths, please set it")
SET(Boost_LIBRARY_DIRS CACHE PATH "The directory containing the boost libraries")
- endif(WIN32)
- link_directories(${Boost_LIBRARY_DIRS})
+ endif(WIN32)
+ list(APPEND AMDiS_LINK_DIRS ${Boost_LIBRARY_DIRS})
+# link_directories(${Boost_LIBRARY_DIRS})
endif(Boost_FOUND)
SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
@@ -120,6 +135,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
${SOURCE_DIR}/SecondOrderTerm.cc
${SOURCE_DIR}/Serializer.cc
${SOURCE_DIR}/StandardProblemIteration.cc
+ ${SOURCE_DIR}/SolverMatrix.cc
${SOURCE_DIR}/SubAssembler.cc
${SOURCE_DIR}/SubElInfo.cc
${SOURCE_DIR}/SubQuadrature.cc
@@ -154,10 +170,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
if(ENABLE_PARALLEL_DOMAIN)
find_package(MPI REQUIRED)
if(MPI_FOUND)
- #INCLUDE (CMakeForceCompiler)
- #CMAKE_FORCE_C_COMPILER(mpicc "MPI C compiler")
- #CMAKE_FORCE_CXX_COMPILER(mpiCC "MPI C++ compiler")
- list(APPEND AMDiS_COMPILEFLAGS ${MPI_COMPILE_FLAGS})
+ list(APPEND COMPILEFLAGS "${MPI_COMPILE_FLAGS}")
include_directories(${MPI_INCLUDE_PATH})
endif(MPI_FOUND)
@@ -165,7 +178,10 @@ if(ENABLE_PARALLEL_DOMAIN)
make
WORKING_DIRECTORY ${LIB_DIR}/ParMetis-3.1
)
- include_directories( ${LIB_DIR}/ParMetis-3.1)
+#include_directories( ${LIB_DIR}/ParMetis-3.1)
+ list(APPEND AMDiS_INCLUDE_DIRS ${LIB_DIR}/ParMetis-3.1)
+ list(APPEND AMDiS_LIBS ${LIB_DIR}/ParMetis-3.1/libparmetis.a ${LIB_DIR}/ParMetis-3.1/libmetis.a)
+# list(APPEND AMDiS_LINK_DIRS ${LIB_DIR}/ParMetis-3.1/)
INSTALL(FILES ${LIB_DIR}/ParMetis-3.1/parmetis.h
${LIB_DIR}/ParMetis-3.1/libparmetis.a
${LIB_DIR}/ParMetis-3.1/libmetis.a
@@ -180,14 +196,12 @@ if(ENABLE_PARALLEL_DOMAIN)
else()
message( FATAL_ERROR "could not find zoltan.")
endif( ZOLTAN_HEADER_FILE )
- set(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_ZOLTAN")
+ list(APPEND COMPILEFLAGS "-DHAVE_ZOLTAN")
endif(ENABLE_ZOLTAN)
- set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};/usr/share/cmake-2.8/Modules/;${CMAKE_SOURCE_DIR}/")
- set(PETSC_EXECUTABLE_RUNS ON)
- find_package(PETSc REQUIRED)
- include_directories(${PETSC_DIR}/include ${PETSC_DIR}/${PETSC_ARCH}/include)
+ list(APPEND COMPILEFLAGS "-DHAVE_PARALLEL_DOMAIN_AMDIS=1")
SET(PARALLEL_DOMAIN_AMDIS_SRC
+ ${SOURCE_DIR}/parallel/ParMetisPartitioner.cc
${SOURCE_DIR}/parallel/CheckerPartitioner.cc
${SOURCE_DIR}/parallel/ElementObjectData.cc
${SOURCE_DIR}/parallel/MeshDistributor.cc
@@ -197,22 +211,41 @@ if(ENABLE_PARALLEL_DOMAIN)
${SOURCE_DIR}/parallel/ParallelDebug.cc
${SOURCE_DIR}/parallel/ParallelProblemStatBase.cc
${SOURCE_DIR}/parallel/ParMetisPartitioner.cc
- ${SOURCE_DIR}/parallel/PetscProblemStat.cc
- ${SOURCE_DIR}/parallel/PetscSolver.cc
- ${SOURCE_DIR}/parallel/PetscSolverFeti.cc
- ${SOURCE_DIR}/parallel/PetscSolverGlobalMatrix.cc
- ${SOURCE_DIR}/parallel/PetscSolverSchur.cc
${SOURCE_DIR}/parallel/StdMpi.cc
${SOURCE_DIR}/parallel/ZoltanPartitioner.cc)
- SET(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_PARALLEL_DOMAIN_AMDIS=1")
+
+ if(ENABLE_PARALLEL_DOMAIN STREQUAL "PETSC")
+ set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};/usr/share/cmake-2.8/Modules/;${CMAKE_SOURCE_DIR}/")
+ set(PETSC_EXECUTABLE_RUNS ON)
+ find_package(PETSc REQUIRED)
+ include_directories(${PETSC_DIR}/include ${PETSC_DIR}/${PETSC_ARCH}/include)
+ list(APPEND AMDiS_INCLUDE_DIRS ${PETSC_DIR}/include ${PETSC_DIR}/${PETSC_ARCH}/include)
+ list(APPEND PARALLEL_DOMAIN_AMDIS_SRC ${SOURCE_DIR}/parallel/PetscSolver.cc
+ ${SOURCE_DIR}/parallel/PetscProblemStat.cc
+ ${SOURCE_DIR}/parallel/PetscSolverFeti.cc
+ ${SOURCE_DIR}/parallel/PetscSolverGlobalMatrix.cc
+ ${SOURCE_DIR}/parallel/PetscSolverSchur.cc)
+ elseif(ENABLE_PARALLEL_DOMAIN STREQUAL "PMTL")
+ set(MTL_INCLUDE_DIR "")
+ find_package(MTL REQUIRED)
+ list(APPEND COMPILEFLAGS "-DMTL_HAS_MPI" "-DHAVE_PARALLEL_MTL4")
+ list(APPEND PARALLEL_DOMAIN_AMDIS_SRC ${SOURCE_DIR}/parallel/Mtl4Solver.cc)
+ list(APPEND AMDiS_LIBS ${MTL_LIBRARIES})
+#message("using the parallel mtl4 version")
+
+ else()
+ message(SEND_ERROR "please set the domain-decomposition version correctly")
+ endif()
+
INSTALL(FILES ResolveCompilerPaths.cmake FindPackageMultipass.cmake FindPETSc.cmake
DESTINATION share/amdis/)
if( USE_PETSC_DEV )
- SET(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_PETSC_DEV")
+ list(APPEND COMPILEFLAGS "-DHAVE_PETSC_DEV")
endif( USE_PETSC_DEV)
endif(ENABLE_PARALLEL_DOMAIN)
if(ENABLE_UMFPACK)
+
# include_directories(${LIB_DIR}/UFconfig ${LIB_DIR}/AMD/Include ${LIB_DIR}/UMFPACK/Include)
# FILE(GLOB UMFPACK_HEADERS "${LIB_DIR}/UMFPACK/Include/*.h")
# INSTALL(FILES ${UMFPACK_HEADERS}
@@ -281,7 +314,7 @@ if(ENABLE_UMFPACK)
#add the compileflags and directories
include_directories(${UMFPACK_PATH} ${UFCONFIG_PATH} ${AMD_PATH})
- SET(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_UMFPACK=1 -DMTL_HAS_UMFPACK")
+ list(APPEND COMPILEFLAGS "-DHAVE_UMFPACK=1" "-DMTL_HAS_UMFPACK")
else()
message(FATAL_ERROR "Could not find the umfpack-headers. Please install umfpack")
endif()
@@ -304,20 +337,25 @@ file(GLOB REINIT_SRC ${REINIT_SOURCE_DIR}/*.cc)
include_directories(${REINIT_SOURCE_DIR})
#mtl4 includes
-include_directories(${MTL_DIR})
-include_directories(${SOURCE_DIR})
+list(APPEND AMDiS_INCLUDE_DIRS ${MTL_INCLUDE_DIR})
+#include_directories(${MTL_INCLUDE_DIR})
+list(APPEND AMDiS_INCLUDE_DIRS ${SOURCE_DIR})
+#include_directories(${SOURCE_DIR})
+include_directories(${AMDiS_INCLUDE_DIRS})
add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC})
add_library(compositeFEM SHARED ${COMPOSITE_FEM_SRC})
add_library(reinit STATIC ${REINIT_SRC})
target_link_libraries(compositeFEM amdis)
+list(APPEND AMDiS_LIBS amdis ${Boost_LIBRARIES})
if(WIN32)
- SET(COMPILEFLAGS "${COMPILEFLAGS} -D_SCL_SECURE_NO_WARNINGS -D_CRT_SECURE_NO_WARNINGS")
+ list(APPEND COMPILEFLAGS "-D_SCL_SECURE_NO_WARNINGS" "-D_CRT_SECURE_NO_WARNINGS")
endif(WIN32)
message("compileflags: ${COMPILEFLAGS}")
#SET_TARGET_PROPERTIES(amdis PROPERTIES COMPILE_FLAGS "${COMPILEFLAGS}")
add_definitions(${COMPILEFLAGS})
+
if(ENABLE_MARMOT)
Message("please set marmotcc manually")
endif(ENABLE_MARMOT)
@@ -391,15 +429,6 @@ endforeach(mtl4_file)
list(REMOVE_DUPLICATES deb_add_dirs)
#PACKAGES
-#let cmake look for the svn-version
-if( NOT CurrentRevision )
- find_package(Subversion)
- if(Subversion_FOUND)
- Subversion_WC_INFO(${AMDiS_SOURCE_DIR} mtlSubinfo)
- set(CurrentRevision "0.9.${mtlSubinfo_WC_REVISION}")
- # message("current revision: ${mtlSubinfo_WC_REVISION}")
- endif(Subversion_FOUND)
-endif( NOT CurrentRevision )
set(CPACK_PACKAGE_NAME "AMDiS")
set(CPACK_PACKAGE_CONTACT "Andreas.Naumann@tu-dresden.de")
set(CPACK_PACKAGE_VERSION "${CurrentRevision}")
@@ -418,7 +447,3 @@ set(CPACK_DEBIAN_PACKAGE_DEPENDS "libboost-dev (>= 1.42), libboost-iostreams-dev
set(CPACK_DEBIAN_PACKAGE_CONTROL_EXTRA "${AMDiS_BINARY_DIR}/preinst" "${AMDiS_BINARY_DIR}/postrm")
include(CPack)
-if(ENABLE_TESTS)
-ENABLE_TESTING()
-add_subdirectory(test)
-endif(ENABLE_TESTS)
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 6b052c5ec9a4d290c075cc25fe7fa83b2965220c..a33921472fd35e0fc0b3c14d99c79b40ce2ed07c 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -23,6 +23,8 @@
#ifndef AMDIS_H
#define AMDIS_H
+#include "stdint.h"
+#include "MTL4Types.h"
#include "AbstractFunction.h"
#include "AdaptInfo.h"
#include "AdaptInstationary.h"
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index 266bf61337695035d2c90cf54fe14a35328cfe77..a107029556c2f3aaf6bcb865d9b4c8c16883b6c2 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -53,7 +53,6 @@ namespace AMDiS {
class Flag;
class IdentityPreconditioner;
class InteriorBoundary;
- class ITL_BasePreconditioner;
class LeafDataPeriodic;
class LevelAdmin;
class Line;
@@ -120,7 +119,7 @@ namespace AMDiS {
template<typename T> class DimVec;
template<typename T> class DimMat;
template<typename ITLSolver> class ITL_OEMSolver;
- template<typename T> class ITL_Preconditioner;
+ template<typename T, typename Vec, typename Mat > class ITL_Preconditioner;
template<typename T> class Matrix;
template<typename T> class MatrixOfFixVecs;
template<typename T> class MatVecMultiplier;
diff --git a/AMDiS/src/BasePreconditioner.h b/AMDiS/src/BasePreconditioner.h
new file mode 100644
index 0000000000000000000000000000000000000000..0161895c6cbff29f3f3e02c2f775c22e6f0fea78
--- /dev/null
+++ b/AMDiS/src/BasePreconditioner.h
@@ -0,0 +1,46 @@
+#ifndef BASEPRECONDITIONER_H
+#define BASEPRECONDITIONER_H
+
+#include <stdexcept>
+#include "CreatorInterface.h"
+#include "ITL_Preconditioner.h"
+#include "Collection.h"
+
+namespace AMDiS {
+
+ template< typename Vector >
+ struct CreatorInterface< ITL_BasePreconditioner< Vector > > {
+ typedef typename AMDiS::Collection< Vector >::PreconditionMatrix PreconditionMatrix;
+ virtual ITL_BasePreconditioner< Vector >* create(const PreconditionMatrix& )=0 /*{ return NULL; }*/ ;
+ };
+
+ template< typename Vector >
+ class NullCreator< ITL_BasePreconditioner< Vector > > :
+ public CreatorInterface< ITL_BasePreconditioner< Vector > >
+ {
+ typedef typename AMDiS::Collection< Vector >::PreconditionMatrix PreconditionMatrix;
+ typedef ITL_BasePreconditioner< Vector > BaseClass;
+ /// Creates no object.
+ BaseClass* create()
+ {
+ TEST_EXIT(false)("should never class this function");
+ return NULL;
+ }
+
+ BaseClass* create(const PreconditionMatrix& ) {
+ TEST_EXIT(false)("this is the wrong null creator\n");
+ return NULL;
+ }
+ ///
+ virtual bool isNullCreator()
+ {
+ return true;
+ }
+ };
+
+ typedef CreatorInterface< ITL_BasePreconditioner< MTLTypes::MTLVector > > PreconditionCreator;
+ typedef ITL_BasePreconditioner< MTLTypes::MTLVector > BasePreconditioner;
+
+}
+
+#endif
diff --git a/AMDiS/src/Collection.h b/AMDiS/src/Collection.h
new file mode 100644
index 0000000000000000000000000000000000000000..4e2c5c1be3e006aa2e5087fc783eba0ce6250390
--- /dev/null
+++ b/AMDiS/src/Collection.h
@@ -0,0 +1,41 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == TU Dresden ==
+// == ==
+// == Institut für Wissenschaftliches Rechnen ==
+// == Zellescher Weg 12-14 ==
+// == 01069 Dresden ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
+// == ==
+// ============================================================================
+
+#ifndef AMDIS_COLLECTION_H
+#define AMDIS_COLLECTION_H
+
+#include "MTL4Types.h"
+
+namespace AMDiS {
+ template< typename T >
+ struct Collection {};
+
+ using namespace MTLTypes;
+ template< >
+ struct Collection< MTLMatrix > {
+ typedef mtl::matrix::inserter< MTLMatrix > Inserter;
+ };
+
+ template< >
+ struct Collection< MTLVector > {
+ typedef mtl::vector::inserter< MTLVector > Inserter;
+ typedef MTLMatrix PreconditionMatrix;
+ };
+}
+#endif
diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc
index f9e434590fd330f99d2cd3dbb31984c68746a507..70091af41569a29112569ae36ca91eb1b85268ce 100644
--- a/AMDiS/src/CreatorMap.cc
+++ b/AMDiS/src/CreatorMap.cc
@@ -11,6 +11,7 @@
#include "CreatorMap.h"
+#include "MTL4Types.h"
#include "OEMSolver.h"
#include "ITL_Solver.h"
#include "ITL_Preconditioner.h"
@@ -74,9 +75,9 @@ namespace AMDiS {
template<>
- void CreatorMap<ITL_BasePreconditioner>::addDefaultCreators()
+ void CreatorMap<BasePreconditioner>::addDefaultCreators()
{
- CreatorInterface<ITL_BasePreconditioner> *creator;
+ PreconditionCreator *creator;
creator = new DiagonalPreconditioner::Creator;
addCreator("diag", creator);
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index 63e4d733727f807be44b9c0e42fdfbfb18f5d64a..8e1bd8a2f39ea8af80f797ece412ee2bd59a514d 100644
--- a/AMDiS/src/Global.h
+++ b/AMDiS/src/Global.h
@@ -37,7 +37,9 @@
#endif
/** \brief current AMDiS version */
-#define AMDIS_VERSION "AMDiS: Version 0.100"
+#ifndef AMDIS_VERSION
+#define AMDIS_VERSION "AMDiS: Version 0.9.1"
+#endif
#include <string>
#include <vector>
diff --git a/AMDiS/src/ITL_OEMSolver.h b/AMDiS/src/ITL_OEMSolver.h
index 7af26e23d5962fcf857736299c0792f8734971d8..6889a75efedb1f2824e312fb923d1bc5e55e2b50 100644
--- a/AMDiS/src/ITL_OEMSolver.h
+++ b/AMDiS/src/ITL_OEMSolver.h
@@ -24,55 +24,36 @@
#ifndef AMDIS_ITL_OEM_SOLVER_H
#define AMDIS_ITL_OEM_SOLVER_H
-#include <iostream>
-#include "OEMSolver.h"
-#include "ITL_Preconditioner.h"
-#include "SolverMatrix.h"
-#include "DOFVector.h"
-#include "SystemVector.h"
-#include "DOFMatrix.h"
-#include "Parameters.h"
+#include "MTL4Solver.h"
+#include "ITL_OEMSolver.hh"
+#include "MTL4Types.h"
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/mtl/mtl.hpp>
namespace AMDiS {
-
/**
* \ingroup Solver
*
*\brief
- * Template for using solvers from ITL in AMDiS.
+ * Template for using solvers from ITL in AMDiS.
+ * Uses a previous definition of MTLMatrix and MTLVector.
+ * Additional template generalization of ITL_OEMSolver destroys class hierarchy
*/
+
template <typename ITLSolver>
- class ITL_OEMSolver : public OEMSolver
- {
- typedef ITL_Preconditioner<itl::pc::identity<DOFMatrix> > id_t;
+ class ITL_OEMSolver : public MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > > {
- public:
+ protected:
+ public:
/// The constructor reads needed parameters and sets solvers \ref name.
- ITL_OEMSolver(std::string name) : OEMSolver(name) {}
+ ITL_OEMSolver(std::string name) :
+ MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > >(name) {}
+
~ITL_OEMSolver() {}
-
- /// Solves the system iteratively
- // This function is virtual because derived class calls solver with parameter
- virtual int solveSystem(const DOFMatrix::base_matrix_type& A,
- mtl::dense_vector<value_type>& x,
- const mtl::dense_vector<value_type>& b)
- {
- itl::cyclic_iteration<value_type> iter(b, this->max_iter, this->relative,
- this->tolerance, this->print_cycle);
-
- error = ITLSolver()(A, x, b, *this->leftPrecon, *this->rightPrecon, iter);
- iterations = iter.iterations();
- residual = iter.resid();
- return error;
- }
-
-
/// Creator class used in the OEMSolverMap.
class Creator : public OEMSolverCreator
{
@@ -88,7 +69,6 @@ namespace AMDiS {
};
-
/**
* \ingroup Solver
*
@@ -97,34 +77,19 @@ namespace AMDiS {
*/
template <typename ITLSolver>
class ITL_OEMSolver_para
- : public OEMSolver // ITL_OEMSolver<ITLSolver>
+ : public MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > > // ITL_OEMSolver<ITLSolver>
{
- typedef DOFMatrix::value_type value_type;
public:
/// The constructor reads needed parameters and sets solvers \ref name.
ITL_OEMSolver_para(std::string name)
- : OEMSolver(name), ell(OEMSolver::max_iter)
+ : MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > >(name)
{
- Parameters::get(name + "->ell", ell);
+
}
~ITL_OEMSolver_para() {}
/// Set parameter of iterative solver
- void setEll(int p) { ell= p; }
-
- /// Solves the system iteratively
- int solveSystem(const DOFMatrix::base_matrix_type& A,
- mtl::dense_vector<value_type>& x,
- const mtl::dense_vector<value_type>& b)
- {
- itl::cyclic_iteration<value_type> iter(b, this->max_iter, this->relative,
- this->tolerance, this->print_cycle);
- error = ITLSolver()(A, x, b, *this->leftPrecon, *this->rightPrecon, iter, ell);
- iterations = iter.iterations();
- residual = iter.resid();
- return error;
- }
class Creator : public OEMSolverCreator
{
@@ -139,11 +104,7 @@ namespace AMDiS {
};
private:
- /// parameter for bicgstab_ell iterative solver
- int ell;
};
-
-
} // namespace AMDiS
#endif // AMDIS_ITL_OEM_SOLVER_H
diff --git a/AMDiS/src/ITL_OEMSolver.hh b/AMDiS/src/ITL_OEMSolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..248584d3545f1884faff5d90da55ecbd989db4e0
--- /dev/null
+++ b/AMDiS/src/ITL_OEMSolver.hh
@@ -0,0 +1,107 @@
+#include "ITL_Preconditioner.h"
+
+namespace AMDiS {
+ template< typename Vector >
+ struct PreconPair {
+ /// Pointer to the left preconditioner
+ ITL_BasePreconditioner< Vector >* l;
+
+ /// Pointer to the right preconditioner
+ ITL_BasePreconditioner< Vector >* r;
+ PreconPair():
+ l(NULL), r(NULL)
+ {}
+
+ };
+
+ struct BaseITL_runner {
+ OEMSolver& oem;
+ BaseITL_runner(OEMSolver* oemPtr):
+ oem(*oemPtr)
+ {
+ TEST_EXIT(oemPtr != NULL)("need a real oemsolver\n");
+ }
+
+ template< typename Vector , typename Matrix >
+ void setPrecon(PreconPair< Vector >& pair, const Matrix& matrix)
+ {
+ /// Creator for the left preconditioner
+ CreatorInterface< ITL_BasePreconditioner< Vector > >* leftCreator(NULL);
+
+ /// Creator for the right preconditioner
+ CreatorInterface< ITL_BasePreconditioner< Vector > >* rightCreator(NULL);
+
+ std::string preconType("no");
+ Parameters::get(oem.getName() + "->left precon", preconType);
+
+ leftCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
+ TEST_EXIT(leftCreator != NULL)
+ ("there is no creator for the given left preconditioner");
+
+ preconType = "no";
+ Parameters::get(oem.getName() + "->right precon", preconType);
+ rightCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
+ TEST_EXIT(rightCreator != NULL)
+ ("there is no creator for the given right preconditioner");
+
+ pair.l = leftCreator->create(matrix);
+ pair.r = rightCreator->create(matrix);
+ }
+ };
+
+ template< typename ITLSolver >
+ struct ITL_OEMSolver_para_runner : BaseITL_runner {
+ ITL_OEMSolver_para_runner(OEMSolver* oem_):
+ BaseITL_runner(oem_),
+ ell(1)
+ {
+ ell=oem.getMaxIterations();
+ Parameters::get(oem.getName() + "->ell", ell);
+ }
+
+ template< typename Matrix, typename Vector>
+ int solve(const Matrix& A, Vector& x, Vector& b)
+ {
+ itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(),
+ oem.getTolerance(), oem.getPrint_cycle());
+
+ static PreconPair< Vector > preconPair;
+ if (preconPair.l == NULL || preconPair.r == NULL || !oem.getMultipleRhs())
+ BaseITL_runner::setPrecon(preconPair, A);
+
+ int error = ITLSolver()(A, x, b, *(preconPair.l), *(preconPair.r), iter, ell);
+ oem.setIterations(iter.iterations());
+ oem.setResidual(iter.resid());
+ oem.setErrorCode(error);
+ return error;
+ }
+ private:
+ int ell;
+ };
+
+ template< typename ITLSolver >
+ struct ITL_OEMSolver_runner : BaseITL_runner {
+
+
+ ITL_OEMSolver_runner(OEMSolver* oem):
+ BaseITL_runner(oem)
+ {}
+
+ template< typename Matrix, typename Vector>
+ int solve(const Matrix& A, Vector& x, Vector& b)
+ {
+ itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(),
+ oem.getTolerance(), oem.getPrint_cycle());
+
+ static PreconPair< Vector > preconPair;
+ if (preconPair.l == NULL || preconPair.r == NULL || !oem.getMultipleRhs())
+ BaseITL_runner::setPrecon(preconPair, A);
+
+ int error = ITLSolver()(A, x, b, *(preconPair.l), *(preconPair.r), iter);
+ oem.setErrorCode(error);
+ oem.setIterations(iter.iterations());
+ oem.setResidual(iter.resid());
+ return error;
+ }
+ };
+}
diff --git a/AMDiS/src/ITL_Preconditioner.h b/AMDiS/src/ITL_Preconditioner.h
index d62f87f5d4daae0dc76b29de1bfef7c84306ca54..39a30f8d7a35f94e461d32eb2cccec664ffb619b 100644
--- a/AMDiS/src/ITL_Preconditioner.h
+++ b/AMDiS/src/ITL_Preconditioner.h
@@ -23,6 +23,7 @@
#ifndef AMDIS_ITL_PRECONDITIONER_H
#define AMDIS_ITL_PRECONDITIONER_H
+#include "MTL4Types.h"
#include "DOFMatrix.h"
#include "CreatorInterface.h"
@@ -42,58 +43,56 @@ namespace AMDiS {
*
* \brief Common base class for wrappers to use ITL preconditioners in AMDiS.
*/
+ template< class Vec >
class ITL_BasePreconditioner
{
public:
- typedef DOFMatrix::value_type value_type;
-
virtual ~ITL_BasePreconditioner() {}
- virtual mtl::dense_vector<value_type>
- member_solve(const mtl::dense_vector<value_type>& vin) const = 0;
+ virtual Vec member_solve(const Vec& vin) const =0 ;
- virtual mtl::dense_vector<value_type>
- member_adjoint_solve(const mtl::dense_vector<value_type>& vin) const = 0;
+ virtual Vec member_adjoint_solve(const Vec& vin) const =0 ;
};
- inline mtl::dense_vector<DOFMatrix::value_type>
- solve(const ITL_BasePreconditioner& P, const mtl::dense_vector<DOFMatrix::value_type>& vin)
- {
- return P.member_solve(vin);
- }
+ template< typename Vec >
+ inline Vec
+ solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
+ {
+ return P.member_solve(vin);
+ }
- inline mtl::dense_vector<DOFMatrix::value_type>
- adjoint_solve(const ITL_BasePreconditioner& P, const mtl::dense_vector<DOFMatrix::value_type>& vin)
- {
- return P.member_adjoint_solve(vin);
- }
+ template< typename Vec >
+ inline Vec
+ adjoint_solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
+ {
+ return P.member_adjoint_solve(vin);
+ }
/**
* \ingroup Solver
*
* \brief Wrapper for using ITL preconditioners in AMDiS.
*/
- template <typename Preconditioner>
- class ITL_Preconditioner : public ITL_BasePreconditioner
+ template <typename Preconditioner, typename MTLVector, typename MTLMatrix >
+ class ITL_Preconditioner : public ITL_BasePreconditioner< MTLVector >
{
- typedef DOFMatrix::value_type value_type;
public:
- ITL_Preconditioner(const DOFMatrix::base_matrix_type& A) : precon(A) {}
+ ITL_Preconditioner(const MTLMatrix& A) : precon(A) {}
- mtl::dense_vector<value_type>
- member_solve(const mtl::dense_vector<value_type>& vin) const
+ MTLVector
+ member_solve(const MTLVector& vin) const
{
return solve(precon, vin);
}
- mtl::dense_vector<value_type>
- member_adjoint_solve(const mtl::dense_vector<value_type>& vin) const
+ MTLVector
+ member_adjoint_solve(const MTLVector& vin) const
{
return adjoint_solve(precon, vin);
}
/// Creator class used in the OEMSolverMap.
- class Creator : public CreatorInterface<ITL_BasePreconditioner>
+ class Creator : public CreatorInterface< ITL_BasePreconditioner< MTLVector > >
{
public:
virtual ~Creator() {}
@@ -101,12 +100,8 @@ namespace AMDiS {
/** \brief
* Creates an ITL preconditioner
*/
- ITL_BasePreconditioner *create(const DOFMatrix::base_matrix_type& A) {
- return new ITL_Preconditioner<Preconditioner>(A);
- }
-
- ITL_BasePreconditioner *create() {
- ERROR_EXIT("Must not be called! Only defined to avoid abstract function."); return 0;
+ ITL_BasePreconditioner< MTLVector >* create(const MTLMatrix& A) {
+ return new ITL_Preconditioner<Preconditioner, MTLVector, MTLMatrix>(A);
}
};
@@ -126,12 +121,12 @@ namespace AMDiS {
* Diagonal preconditioner.
*/
class DiagonalPreconditioner
- : public ITL_Preconditioner<itl::pc::diagonal<DOFMatrix::base_matrix_type> >
+ : public ITL_Preconditioner<itl::pc::diagonal<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
{
- typedef ITL_Preconditioner<itl::pc::diagonal<DOFMatrix::base_matrix_type> > base;
+ typedef ITL_Preconditioner<itl::pc::diagonal<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix > base;
public:
- DiagonalPreconditioner(const DOFMatrix::base_matrix_type& A) : base(A) {}
+ DiagonalPreconditioner(const MTLTypes::MTLMatrix& A) : base(A) {}
};
/**
@@ -141,12 +136,12 @@ namespace AMDiS {
* Identity preconditioner. Behaves like no preconditioning.
*/
class IdentityPreconditioner
- : public ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type> >
+ : public ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
{
- typedef ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type> > base;
+ typedef ITL_Preconditioner<itl::pc::identity<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix > base;
public:
- IdentityPreconditioner(const DOFMatrix::base_matrix_type& A) : base(A) {}
+ IdentityPreconditioner(const MTLTypes::MTLMatrix& A) : base(A) {}
};
@@ -165,11 +160,11 @@ namespace AMDiS {
* described in chapter 10.3 (algorithm 10.4).
*/
class ILUPreconditioner
- : public ITL_Preconditioner< itl::pc::ilu_0<DOFMatrix::base_matrix_type> >
+ : public ITL_Preconditioner< itl::pc::ilu_0<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
{};
class ICPreconditioner
- : public ITL_Preconditioner< itl::pc::ic_0<DOFMatrix::base_matrix_type> >
+ : public ITL_Preconditioner< itl::pc::ic_0<DOFMatrix::base_matrix_type>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
{};
diff --git a/AMDiS/src/ITL_Solver.h b/AMDiS/src/ITL_Solver.h
index 787903da60917538cea2e8904ec5a7aa29866699..9b283fbb26a21cf14ef270750ded4af9b38fc8b2 100644
--- a/AMDiS/src/ITL_Solver.h
+++ b/AMDiS/src/ITL_Solver.h
@@ -16,9 +16,7 @@
//
// See also license.opensource.txt in the distribution.
-
-
-/** \file ITL_Preconditioner.h */
+/** \file ITL_Solver.h */
#ifndef AMDIS_ITL_SOLVER_H
#define AMDIS_ITL_SOLVER_H
diff --git a/AMDiS/src/MTL4Solver.h b/AMDiS/src/MTL4Solver.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6ba1907409fae38e74baf150f56128a822be840
--- /dev/null
+++ b/AMDiS/src/MTL4Solver.h
@@ -0,0 +1,79 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == TU Dresden ==
+// == ==
+// == Institut für Wissenschaftliches Rechnen ==
+// == Zellescher Weg 12-14 ==
+// == 01069 Dresden ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
+// == ==
+// ============================================================================
+
+#ifndef MTL4SOLVER_H
+#define MTL4SOLVER_H
+
+#include "OEMSolver.h"
+#include "MatrixStreams.h"
+#include <iostream>
+
+namespace AMDiS {
+ template< typename MTLMatrix, typename MTLVector, typename Worker >
+ class MTL4Solver : public OEMSolver {
+
+ protected:
+
+ Worker worker;
+
+ template< typename Matrix, typename Vector, typename Mapper >
+ int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper)
+ {
+ MTLMatrix matrix(mapper.nRow(), mapper.nCol());
+ set_to_zero(matrix);
+ MatMap< const Matrix, Mapper > matMap(A,mapper);
+ matrix << matMap;
+
+ MTLVector mtl_x(mapper.nRow());
+ set_to_zero(mtl_x);
+ VecMap< Vector, Mapper > xVecMap(x, mapper);
+ mtl_x << xVecMap;
+
+ MTLVector mtl_b(mapper.nRow());
+ set_to_zero(mtl_b);
+ VecMap< Vector, Mapper> bVecMap(b, mapper);
+ mtl_b << bVecMap;
+
+ error = worker.solve(matrix, mtl_x, mtl_b);
+
+ mtl_x >> xVecMap;
+ return error;
+ }
+
+ public:
+ MTL4Solver(std::string n):
+ OEMSolver(n),
+ worker(this)
+ {
+ FUNCNAME("MTL4Solver()");
+
+ }
+
+ virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+ SystemVector& x,
+ SystemVector& b,
+ VectorialMapper& m)
+ {
+ return solve(A,x,b,m);
+ }
+
+
+ };
+}
+#endif
diff --git a/AMDiS/src/MTL4Types.h b/AMDiS/src/MTL4Types.h
new file mode 100644
index 0000000000000000000000000000000000000000..be018fd1a0a8de4a91e3aaa83f527ef1146a7273
--- /dev/null
+++ b/AMDiS/src/MTL4Types.h
@@ -0,0 +1,12 @@
+#ifndef MTL4TYPES_H
+#define MTL4TYPES_H
+#include <boost/numeric/mtl/mtl.hpp>
+namespace AMDiS {
+ namespace MTLTypes {
+ typedef double value_type;
+ typedef mtl::matrix::compressed2D< value_type > MTLMatrix;
+ typedef mtl::vector::dense_vector< value_type > MTLVector;
+ }
+}
+#endif
+
diff --git a/AMDiS/src/Mapper.h b/AMDiS/src/Mapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..e16230ce0259b99bc6b9ba33adcd318fde206e54
--- /dev/null
+++ b/AMDiS/src/Mapper.h
@@ -0,0 +1,78 @@
+#ifndef MAPPER_H
+#define MAPPER_H
+
+#include "SolverMatrix.h"
+#include <vector>
+namespace AMDiS {
+/* struct ScalarMapper {
+ ScalarMapper(int nrow, int ncol):
+ nrow(nrow), ncol(ncol) {}
+ inline int row(int r) { return r; }
+ inline int col(int c) { return c; }
+ inline int nRow() { return nrow; }
+ inline int nCol() { return ncol; }
+ protected:
+ int nrow;
+ int ncol;
+ };
+*/
+ struct VectorialMapper
+ {
+
+ VectorialMapper(const SolverMatrix<Matrix<DOFMatrix* > >& sm ) :
+ nComp(sm.getOriginalMat()->getSize()),
+ rowOffset(0), colOffset(0), nrow(0), ncol(0), sizes(nComp)
+ {
+ const Matrix<DOFMatrix* >& orMat(*sm.getOriginalMat());
+ const int ns = orMat.getSize();
+ for (int i= 0; i < ns; i++) {
+ sizes[i] = orMat[i][i]->getFeSpace()->getAdmin()->getUsedSize();
+ nrow += sizes[i];
+ }
+ ncol = nrow;
+ }
+
+ VectorialMapper(unsigned int nComp, int nDOFperComp) :
+ nComp(nComp), rowOffset(0), colOffset(0),
+ nrow(nComp*nDOFperComp), ncol(nrow),sizes(nComp)
+ {
+ for (unsigned int i(0); i < nComp; ++i)
+ sizes[i] = nDOFperComp;
+ }
+
+ inline void setRow( unsigned int r) {
+ assert( r <= sizes.size() );
+ rowOffset = sum( r);
+ }
+
+ inline void setCol( unsigned int c) { colOffset = sum( c); }
+
+ inline int row(int r) const { return r + rowOffset; }
+
+ inline int col(int c) const { return c + colOffset; }
+
+ inline int nCol() const { return ncol; }
+
+ inline int nRow() const { return nrow; }
+
+ protected:
+
+ ///compute the sum of sizes from [0, end)
+ inline unsigned int sum(unsigned int end) const
+ {
+ unsigned int ret(0);
+ for (unsigned int i(0); i < end; ++i)
+ ret += sizes[i];
+ return ret;
+ }
+
+ int nComp;
+ int rowOffset;
+ int colOffset;
+ int nrow;
+ int ncol;
+ std::vector< unsigned int > sizes;
+ };
+}
+
+#endif
diff --git a/AMDiS/src/MatrixStreams.h b/AMDiS/src/MatrixStreams.h
new file mode 100644
index 0000000000000000000000000000000000000000..98ccbcab3a1779fc52757032d282c48727dfd0ad
--- /dev/null
+++ b/AMDiS/src/MatrixStreams.h
@@ -0,0 +1,141 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == TU Dresden ==
+// == ==
+// == Institut für Wissenschaftliches Rechnen ==
+// == Zellescher Weg 12-14 ==
+// == 01069 Dresden ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
+// == ==
+// ============================================================================
+
+/** \file MatrixStreams.h */
+
+#ifndef MATRIXSTREAMS_H
+#define MATRIXSTREAMS_H
+
+#include "DOFMatrix.h"
+#include "SolverMatrix.h"
+#include "SystemVector.h"
+#include "Collection.h"
+
+using namespace mtl::operations;
+namespace AMDiS {
+ template< typename VecT, typename CurMap>
+ struct VecMap {
+ VecT& vec;
+ CurMap& mapper;
+ VecMap(VecT& vec, CurMap& mapper):
+ vec(vec),mapper(mapper) {}
+ };
+
+ template< typename MatT, typename CurMap>
+ struct MatMap {
+ MatT& mat;
+ CurMap& mapper;
+ MatMap(MatT& mat, CurMap& m):
+ mat(mat), mapper(m) {}
+ };
+
+ template< typename MTLMatrix, typename Mapper >
+ void operator<<(MTLMatrix& matrix, MatMap<const SolverMatrix<Matrix<DOFMatrix* > >, Mapper >& Asolver)
+ {
+ const Matrix< DOFMatrix* >& A = *(Asolver.mat.getOriginalMat());
+ int ns = A.getSize();
+
+ Mapper& mapper(Asolver.mapper);
+ set_to_zero(matrix);
+
+ int nnz = 0;
+ for (int rb = 0; rb < ns; ++rb)
+ for (int cb = 0; cb < ns; ++cb)
+ if (A[rb][cb])
+ nnz += A[rb][cb]->getBaseMatrix().nnz();
+
+ {
+ typedef mtl::matrix::mapped_inserter< typename Collection< MTLMatrix >::Inserter, Mapper > Inserter;
+ Inserter ins(matrix, mapper, int(1.2 * nnz / matrix.dim1()));
+ for (int rb = 0; rb < ns; ++rb) {
+ mapper.setRow(rb);
+ for (int cb = 0; cb < ns; ++cb) {
+ mapper.setCol(cb);
+ if (A[rb][cb])
+ ins << A[rb][cb]->getBaseMatrix();
+ }
+ }
+ }
+ }
+
+ template< typename Vector, typename CurMap >
+ void operator<<(Vector& dest, VecMap<DOFVector<double >, CurMap >& rhs)
+ {
+ DOFVector<double>::Iterator it_x(&rhs.vec, USED_DOFS);
+ int counter(0);
+ typedef typename Collection< Vector >::Inserter Inserter ;
+
+ Inserter swapIns(dest);
+ mtl::vector::mapped_inserter< Inserter, CurMap > ins(swapIns, rhs.mapper);
+
+ for (it_x.reset(); !it_x.end(); ++it_x) {
+ ins[counter] << *it_x;
+ counter++;
+ }
+ }
+
+ template< typename Vector, typename CurMap >
+ inline void operator>>(const Vector& dest, VecMap<DOFVector< double >, CurMap >& rhs)
+ {
+ DOFVector<double>::Iterator it_x(&rhs.vec, USED_DOFS);
+ int counter(0);
+ {
+ mtl::vector::extracter< Vector > extracter(dest);
+ for (it_x.reset(); !it_x.end(); ++it_x) {
+ extracter[rhs.mapper.row(counter)] >> *it_x ;
+ counter++;
+ }
+ }
+ }
+
+ template< typename Vector, typename CurMap >
+ void operator<<(Vector& dest, VecMap<SystemVector, CurMap>& rhs)
+ {
+ int ns = rhs.vec.getSize(); // Number of systems.
+
+ // Copy vectors
+ typedef typename Collection< Vector >::Inserter Inserter;
+ Inserter swapInserter(dest);
+ mtl::vector::mapped_inserter< Inserter, CurMap > ins(swapInserter, rhs.mapper);
+ for (int i = 0; i < ns; i++) {
+ DOFVector<double>::Iterator it_source(rhs.vec.getDOFVector(i), USED_DOFS);
+ int counter(0);
+ rhs.mapper.setRow(i);
+ for (it_source.reset(); !it_source.end(); ++it_source) {
+ ins[counter] << *it_source;
+ counter++;
+ }
+ }
+ }
+
+ template< typename Vector, typename CurMap >
+ void operator>>(const Vector& dest, VecMap<SystemVector, CurMap>& rhs)
+ {
+ int ns = rhs.vec.getSize(); // Number of systems.
+
+ // Copy vectors
+ for (int i = 0; i < ns; i++) {
+ DOFVector< double >& dofvec(*(rhs.vec.getDOFVector(i)));
+ rhs.mapper.setRow(i);
+ VecMap< DOFVector< double >, CurMap > swap(dofvec, rhs.mapper);
+ dest >> swap;
+ }
+ }
+}
+#endif
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index 6c126a0f469ff883305be1bf954f7adf033eba95..ef7fb8950a7ce79e82e8fdf854146ceccf300092 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -39,7 +39,13 @@
#include "DOFVector.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
-#include "ITL_Preconditioner.h"
+#include "BasePreconditioner.h"
+#include "Mapper.h"
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#include "parallel/ParallelMapper.h"
+#endif
+//#include "ITL_Preconditioner.h"
namespace AMDiS {
@@ -67,9 +73,9 @@ namespace AMDiS {
print_cycle(100),
iterations(-1),
error(-1),
- multipleRhs(false),
+ multipleRhs(false)/*,
leftPrecon(NULL),
- rightPrecon(NULL)
+ rightPrecon(NULL)*/
{}
///
@@ -87,75 +93,22 @@ namespace AMDiS {
Parameters::get(name + "->info", info);
}
- /** Set left Preconditioner
- *
- * It is assumed that the solver owns the preconditioner from this call.
- * That means that the preconditioner is deleted by the solver when not needed.
- */
- void setLeftPrecon(ITL_BasePreconditioner* p)
- {
- if (leftPrecon)
- delete leftPrecon;
- leftPrecon = p;
- }
-
- /** Set right Preconditioner
- *
- * It is assumed that the solver owns the preconditioner from this call.
- * That means that the preconditioner is deleted by the solver when not needed.
- */
- void setRightPrecon(ITL_BasePreconditioner* p)
- {
- if (rightPrecon)
- delete rightPrecon;
- rightPrecon = p;
- }
-
- /// Linear System to be solved in the derived class
- virtual int solveSystem(const DOFMatrix::base_matrix_type& A,
- mtl::dense_vector<value_type>& x,
- const mtl::dense_vector<value_type>& b) = 0;
-
-
- /// Solve a linear system for a scalar problem.
- int solveSystem(const SolverMatrix<DOFMatrix>& A,
- DOFVector<double>& x,
- DOFVector<double>& b)
- {
- mtl::dense_vector<value_type> xx(x.getUsedSize()), bb(b.getUsedSize());
-
- // Copy rhs vector
- int counter = 0;
- DOFVector<double>::Iterator it_b(&b, USED_DOFS);
- DOFVector<double>::Iterator it_x(&x, USED_DOFS);
-
- for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {
- bb[counter] = *it_b;
- xx[counter] = *it_x;
- counter++;
- }
-
- int r = solveSystem(A.getMatrix(), xx, bb);
-
- // Copy solution vector to DOFVector
- counter = 0;
- for (it_x.reset(); !it_x.end(); ++it_x)
- *it_x = xx[counter++];
-
- return r;
- }
/// Solve a linear system for a vectorial problem.
virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
- SystemVector& b)
+ SystemVector& b,
+ VectorialMapper&) = 0;
+#if 0
{
int ns = x.getSize(), // Number of systems.
size = x.getUsedSize(); // Size of all DOFVectors
// Copy vectors
mtl::dense_vector<value_type> xx(size), bb(size);
- for (int i = 0, counter = 0; i < ns; i++) {
+ xx << x;
+ bb << b
+/* for (int i = 0, counter = 0; i < ns; i++) {
DOFVector<double>::Iterator it_b(b.getDOFVector(i), USED_DOFS);
DOFVector<double>::Iterator it_x(x.getDOFVector(i), USED_DOFS);
@@ -165,9 +118,11 @@ namespace AMDiS {
counter++;
}
}
-
+ */
// Solver on DOFVector for single system
- int r = solveSystem(A.getMatrix(), xx, bb);
+ DOFMatrix::base_matrix_type baseA;
+ baseA << A;
+ int r = solveSystem(baseA, xx, bb);
for (int i = 0, counter = 0; i < ns; i++) {
DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);
@@ -177,6 +132,27 @@ namespace AMDiS {
return r;
}
+#endif
+
+ inline int solveSystem(const SolverMatrix< Matrix< DOFMatrix* > >& A,
+ SystemVector& x,
+ SystemVector& b)
+ {
+ VectorialMapper mapper(A);
+ return solveSystem(A, x, b, mapper);
+ }
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+ SystemVector& x,
+ SystemVector& b,
+ ParallelMapper&)
+ {
+ FUNCNAME("OEMSolver::solveSystem()");
+ TEST_EXIT(false)("This linear solver is not suitable for domaindecomposition problems\n");
+ return -1;
+ }
+#endif
/** \name getting methods
* \{
@@ -200,12 +176,49 @@ namespace AMDiS {
return max_iter;
}
+ /// Returns number of iterations in last run of an iterative solver
+ inline int getIterations()
+ {
+ return iterations;
+ }
+
+ /// Returns error code in last run of an iterative solver
+ inline int getErrorCode()
+ {
+ return error;
+ }
+
+ /// Returns info
+ inline int getInfo()
+ {
+ return info;
+ }
+
+ /// Returns \ref print_cycle
+ inline int getPrint_cycle() const
+ {
+ return print_cycle;
+ }
+
/// Returns \ref residual
inline double getResidual() const
{
return residual;
}
+ /// Returns \ref relative
+ inline double getRelative() const
+ {
+ return relative;
+ }
+
+ /// Returns \ref multipleRhs
+ inline bool getMultipleRhs() const
+ {
+ return multipleRhs;
+ }
+
+
/** \} */
/** \name setting methods
@@ -218,6 +231,11 @@ namespace AMDiS {
tolerance = tol;
}
+ inline void setResidual(double r)
+ {
+ residual = r;
+ }
+
/// Sets \ref relative
inline void setRelative(double rel)
{
@@ -230,25 +248,24 @@ namespace AMDiS {
max_iter = i;
}
- /// Returns number of iterations in last run of an iterative solver
- inline int getIterations()
+ inline void setIterations(int i)
{
- return iterations;
+ iterations=i;
}
- /// Returns error code in last run of an iterative solver
- inline int getErrorCode()
+ /// set the \ref error
+ inline void setErrorCode(int code)
{
- return error;
+ error=code;
}
-
+
/// Sets \ref info
inline void setInfo(int i)
{
info = i;
}
- void setMultipleRhs(bool b)
+ inline void setMultipleRhs(bool b)
{
multipleRhs = b;
}
@@ -289,12 +306,13 @@ namespace AMDiS {
* as they can do the factorization of the matrix only once.
*/
bool multipleRhs;
-
+#if 0
/// Left preconditioner (ignored by some iterative solvers and by UmfPack)
- ITL_BasePreconditioner *leftPrecon;
+ PreconCreator *leftPrecon;
/// Right preconditioner (ignored by some iterative solvers and by UmfPack)
- ITL_BasePreconditioner *rightPrecon;
+ PreconCreator *rightPrecon;
+#endif
};
/**
@@ -323,6 +341,6 @@ namespace AMDiS {
};
}
-#include "ITL_Solver.h"
+//#include "ITL_Solver.h"
#endif // AMDIS_OEM_SOLVER_H
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index c2e65056dda18d39de31ccd887300039476b7808..711f6a397db4434cdbacf8c716aa3fb277fecda7 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -1290,7 +1290,7 @@ namespace AMDiS {
void ProblemStatSeq::createPrecon()
{
- string preconType("no");
+ /*string preconType("no");
Parameters::get(name + "->solver->left precon", preconType);
CreatorInterface<ITL_BasePreconditioner> *preconCreator =
@@ -1302,7 +1302,7 @@ namespace AMDiS {
Parameters::get(name + "->solver->right precon", preconType);
preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
- solver->setRightPrecon(preconCreator->create(solverMatrix.getMatrix()));
+ solver->setRightPrecon(preconCreator->create(solverMatrix.getMatrix()));*/
}
diff --git a/AMDiS/src/SolverMatrix.cc b/AMDiS/src/SolverMatrix.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c7a1c0b24a51c89c0600b279f0953917316326b7
--- /dev/null
+++ b/AMDiS/src/SolverMatrix.cc
@@ -0,0 +1,15 @@
+#include "SolverMatrix.h"
+#include "Mapper.h"
+#include "MatrixStreams.h"
+
+namespace AMDiS {
+ void SolverMatrix< Matrix<DOFMatrix*> >::buildMatrix() const
+ {
+ VectorialMapper mapper(*this);
+ matrix.change_dim(mapper.nRow(), mapper.nCol());
+ set_to_zero(matrix);
+ MatMap< const SolverMatrix< Matrix<DOFMatrix*> >, VectorialMapper > matMap(*this,mapper);
+ matrix << matMap;
+
+ }
+}
diff --git a/AMDiS/src/SolverMatrix.h b/AMDiS/src/SolverMatrix.h
index 1fe1ae321387cdde86fad71d7e8da88ffa9c8873..9a9e8bde9d9dbad4e2a45a84db160635164bfb6c 100644
--- a/AMDiS/src/SolverMatrix.h
+++ b/AMDiS/src/SolverMatrix.h
@@ -40,8 +40,11 @@ namespace AMDiS {
template <>
class SolverMatrix<DOFMatrix>
{
- public :
- SolverMatrix() : matrix(0) {}
+ public :
+
+ SolverMatrix() :
+ matrix(0)
+ {}
void setMatrix(const DOFMatrix& A)
{
@@ -62,36 +65,22 @@ namespace AMDiS {
template <>
class SolverMatrix<Matrix<DOFMatrix*> >
{
- public :
+ void buildMatrix() const;
+ public :
+ SolverMatrix<Matrix<DOFMatrix* > >():
+ originalMat(NULL)
+ {}
+
void setMatrix(const Matrix<DOFMatrix*>& A)
{
- int ns = A.getSize();
- std::vector<int> block_starts(ns + 1);
-
- block_starts[0] = 0;
- for (int i= 0; i < ns; i++)
- block_starts[i + 1] = block_starts[i] + A[i][i]->getFeSpace()->getAdmin()->getUsedSize();
-
- matrix.change_dim(block_starts[ns], block_starts[ns]);
- set_to_zero(matrix);
-
- int nnz = 0;
- for (int rb = 0; rb < ns; ++rb)
- for (int cb = 0; cb < ns; ++cb)
- if (A[rb][cb])
- nnz += A[rb][cb]->getBaseMatrix().nnz();
-
- DOFMatrix::inserter_type ins(matrix, int(1.2 * nnz / matrix.dim1()));
- for (int rb = 0; rb < ns; ++rb)
- for (int cb = 0; cb < ns; ++cb)
- if (A[rb][cb])
- ins[block_starts[rb]][block_starts[cb]] << A[rb][cb]->getBaseMatrix();
-
originalMat = &A;
+ matrix.change_dim(0,0);
}
-
+
const DOFMatrix::base_matrix_type& getMatrix() const
{
+ if( num_rows(matrix) == 0)
+ buildMatrix();
return matrix;
}
@@ -101,12 +90,11 @@ namespace AMDiS {
}
private:
- DOFMatrix::base_matrix_type matrix;
+ mutable DOFMatrix::base_matrix_type matrix;
const Matrix<DOFMatrix*>* originalMat;
-
- };
+ };
} // namespace AMDiS
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index 520388b478b9ef61d3caee67b325193d5ffd0d20..a07f1b7953bc72e896bb85052660655bc0a3883d 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -28,7 +28,7 @@
#include <iostream>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
-#include "OEMSolver.h"
+#include "MTL4Solver.h"
namespace AMDiS {
@@ -41,81 +41,97 @@ namespace AMDiS {
*
* This is a direct solver for large sparse matrices.
*/
-
- class UmfPackSolver : public OEMSolver
+ template< typename Matrix >
+ struct Umfpack_runner
{
- public:
- /// Creator class used in the OEMSolverMap.
- class Creator : public OEMSolverCreator
- {
- public:
- virtual ~Creator() {}
-
- /// Returns a new UmfPackSolver object.
- OEMSolver* create()
- {
- return new UmfPackSolver(this->name);
- }
- };
-
- /// Constructor
- UmfPackSolver(std::string name)
- : OEMSolver(name),
- solver(0),
- store_symbolic(0),
- symmetric_strategy(0)
+ typedef Matrix matrix_type ;
+ OEMSolver& oem;
+ Umfpack_runner(OEMSolver* oem_):
+ oem(*oem_),
+ solver(NULL),
+ store_symbolic(0),
+ symmetric_strategy(0)
{
- Parameters::get(name + "->store symbolic", store_symbolic);
- Parameters::get(name + "->symmetric strategy", symmetric_strategy);
- }
-
- /// Destructor
- ~UmfPackSolver()
- {
- if (solver)
- delete solver;
+ TEST_EXIT_DBG(oem_ != NULL)("Need real OEMSolver\n");
+ Parameters::get(oem.getName() + "->store symbolic", store_symbolic);
+ Parameters::get(oem.getName() + "->symmetric strategy", symmetric_strategy);
}
- /// Solves the system directly
- int solveSystem(const DOFMatrix::base_matrix_type& A,
- mtl::dense_vector<value_type>& x,
- const mtl::dense_vector<value_type>& b)
- {
+ template< typename Vector>
+ int solve(const Matrix& A, Vector& x, Vector& b)
+ {
if (!solver) {
if (!symmetric_strategy)
solver = new mtl::matrix::umfpack::solver<matrix_type>(A);
else
solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
} else {
- if (!multipleRhs) {
+ if (!oem.getMultipleRhs()) {
+
if (store_symbolic)
solver->update_numeric();
else
solver->update();
}
- }
-
+ }
+
int code = (*solver)(x, b);
- if (info > 0) {
- mtl::dense_vector<value_type> r(b);
+ if( oem.getInfo() > 0) {
+ mtl::dense_vector<typename matrix_type::value_type> r(b);
r -= A * x;
- residual = two_norm(r);
+ double residual = two_norm(r);
+ oem.setResidual(residual);
+ oem.setErrorCode(code);
std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
- if (residual > tolerance) {
- WARNING("Tolerance tol=%e could not be reached!\n",tolerance);
+ if (residual > oem.getTolerance()) {
+ WARNING("Tolerance tol=%e could not be reached!\n", oem.getTolerance());
}
}
return code;
+
}
-
- private:
+ ~Umfpack_runner()
+ {
+ if(!solver)
+ delete solver;
+ }
+
+ private:
mtl::matrix::umfpack::solver<matrix_type> *solver;
int store_symbolic;
-
int symmetric_strategy;
};
+
+ using namespace MTLTypes;
+ class UmfPackSolver : public MTL4Solver<MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >
+ {
+ protected:
+
+ public:
+ /// Creator class used in the OEMSolverMap.
+ class Creator : public OEMSolverCreator
+ {
+ public:
+ virtual ~Creator() {}
+
+ /// Returns a new UmfPackSolver object.
+ OEMSolver* create()
+ {
+ return new UmfPackSolver(this->name);
+ }
+ };
+
+ /// Constructor
+ UmfPackSolver(std::string name)
+ : MTL4Solver< MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >(name)
+ {
+ }
+
+
+ private:
+ };
}
diff --git a/AMDiS/src/parallel/Mtl4Solver.cc b/AMDiS/src/parallel/Mtl4Solver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b7b20d0541b6711d4d3867bb39ae27d7e19795d0
--- /dev/null
+++ b/AMDiS/src/parallel/Mtl4Solver.cc
@@ -0,0 +1,133 @@
+#include "parallel/Mtl4Solver.h"
+#include "parallel/PITL_Solver.h"
+#include "parallel/ParallelMapper.h"
+
+namespace AMDiS {
+ namespace Parallel {
+
+ void Mtl4Solver::addPMTLSolvers()
+ {
+ OEMSolverCreator *creator;
+
+ creator = new CGSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pcg", creator);
+
+ creator = new CGSSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pcgs", creator);
+
+ creator = new BiCGSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pbicg", creator);
+
+ creator = new BiCGStabSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pbicgstab", creator);
+
+ creator = new BiCGStab2Solver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pbicgstab2", creator);
+
+ creator = new BiCGStabEllSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pbicgstab_ell", creator);
+
+ creator = new QMRSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pqmr", creator);
+
+ creator = new TFQMRSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("ptfqmr", creator);
+
+ creator = new GMResSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pgmres", creator);
+
+ creator = new IDRsSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pidr_s", creator);
+
+ creator = new MinResSolver::Creator;
+ CreatorMap< OEMSolver >::addCreator("pminres", creator);
+ }
+
+/* void Mtl4Solver::addPMTLPrecons()
+ {
+ ParallelPreconditionCreator *creator;
+
+ creator = new DiagonalPreconditioner::Creator;
+ CreatorMap<ParallelPreconditioner >::addCreator("pdiag", creator);
+
+ creator = new ILUPreconditioner::Creator;
+ CreatorMap<ParallelPreconditioner >::addCreator("pilu", creator);
+
+ creator = new ICPreconditioner::Creator;
+ CreatorMap<ParallelPreconditioner >::addCreator("pic", creator);
+
+ creator = new IdentityPreconditioner::Creator;
+ CreatorMap<ParallelPreconditioner >::addCreator("pno", creator);
+
+ }*/
+
+ void Mtl4Solver::createSolver()
+ {
+ FUNCNAME("Parallel::Mtl4Solver::createSolver()");
+
+ // === create solver ===
+ std::string solverType("0");
+ Parameters::get(name + "->solver", solverType);
+ solverType = "p" + solverType;
+ OEMSolverCreator *solverCreator =
+ dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
+ TEST_EXIT(solverCreator)("no solver type\n");
+ solverCreator->setName(name + "->solver");
+ solver = solverCreator->create();
+ solver->initParameters();
+ }
+
+ //TODO: replace the name in the map
+ void Mtl4Solver::createPrecon()
+ {
+ /*std::string preconType("no");
+ GET_PARAMETER(0, name + "->solver->left precon", &preconType);
+ preconType = "p" + preconType ;
+
+ CreatorInterface<BasePreconditioner> *preconCreator =
+ CreatorMap<BasePreconditioner>::getCreator(preconType);
+
+ // solver->setLeftPrecon(preconCreator->create(solverMatrix.getMatrix()));
+ solver->setLeftPrecon(preconCreator);
+
+ preconType= "no";
+ GET_PARAMETER(0, name + "->solver->right precon", &preconType);
+ preconType = "p" + preconType;
+
+ preconCreator = CreatorMap<BasePreconditioner>::getCreator(preconType);
+ // solver->setRightPrecon(preconCreator->create(solverMatrix.getMatrix()));
+ solver->setRightPrecon(preconCreator);*/
+ }
+
+ void Mtl4Solver::solve(AdaptInfo* adaptInfo, bool fixedMatrix)
+ {
+/* FUNCNAME("Mtl4Solver::solve()");
+ std::string solverName("");
+ GET_PARAMETER(0, name+"->solver", &solverName);
+ TEST_EXIT(solverName != "")("need solver name");
+ OEMSolver* solver = getPMTLSolver(solverName);
+
+ delete solver;
+ */
+ ParallelMapper pmapper(*meshDistributor, nComponents);
+ solver->solveSystem(solverMatrix, *solution, *rhs, pmapper);
+ }
+ }
+
+ template< >
+ void CreatorMap< Parallel::ParallelPreconditioner >::addDefaultCreators() {
+ Parallel::ParallelPreconditionCreator *creator;
+
+ creator = new Parallel::DiagonalPreconditioner::Creator;
+ addCreator("diag", creator);
+
+ creator = new Parallel::ILUPreconditioner::Creator;
+ addCreator("ilu", creator);
+
+ creator = new Parallel::ICPreconditioner::Creator;
+ addCreator("ic", creator);
+
+ creator = new Parallel::IdentityPreconditioner::Creator;
+ addCreator("no", creator);
+ }
+}
diff --git a/AMDiS/src/parallel/Mtl4Solver.h b/AMDiS/src/parallel/Mtl4Solver.h
index e3043edce8834faab3a026c9b59237404c6186de..30c9d3830b7a7dd19ec23b4d47ad7d77cc0dd551 100644
--- a/AMDiS/src/parallel/Mtl4Solver.h
+++ b/AMDiS/src/parallel/Mtl4Solver.h
@@ -27,26 +27,33 @@
namespace AMDiS {
- class Mtl4Solver : public ParallelProblemStatBase
- {
- public:
- Mtl4Solver(std::string nameStr,
- ProblemIterationInterface *problemIteration = NULL)
- : ParallelProblemStatBase(nameStr, problemIteration)
- {}
+ namespace Parallel {
+ class Mtl4Solver : public ParallelProblemStatBase
+ {
+ void addPMTLSolvers();
+ //void addPMTLPrecons();
+ public:
+ Mtl4Solver(std::string nameStr,
+ ProblemIterationInterface *problemIteration = NULL)
+ : ParallelProblemStatBase(nameStr, problemIteration)
+ {
+ addPMTLSolvers();
+ //addPMTLPrecons();
+ }
- ~Mtl4Solver()
- {}
+ ~Mtl4Solver()
+ {}
- void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false)
- {
- ERROR_EXIT("TO BE IMPLEMENTED!\n");
- }
+ void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false);
+ void createPrecon();
+ void createSolver();
+ };
}
- typedef Mtl4Solver ParallelProblemStat;
+ typedef Parallel::Mtl4Solver ParallelProblemStat;
+ typedef ParallelProblemStat ProblemStat;
}
#endif
diff --git a/AMDiS/src/parallel/PITL_Solver.h b/AMDiS/src/parallel/PITL_Solver.h
new file mode 100644
index 0000000000000000000000000000000000000000..2754dd0e290b059772bafee68b7df02455f1d3c7
--- /dev/null
+++ b/AMDiS/src/parallel/PITL_Solver.h
@@ -0,0 +1,318 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == TU Dresden ==
+// == ==
+// == Institut für Wissenschaftliches Rechnen ==
+// == Zellescher Weg 12-14 ==
+// == 01069 Dresden ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
+// == ==
+// ============================================================================
+
+/** \file PITL_Solver.h */
+
+#ifndef AMDIS_PITL_SOLVER_H
+#define AMDIS_PITL_SOLVER_H
+
+#include "ITL_Solver.h"
+
+namespace AMDiS {
+ namespace Parallel {
+ typedef mtl::matrix::distributed< MTLTypes::MTLMatrix > MTLMatrix;
+ typedef mtl::vector::distributed< MTLTypes::MTLVector > MTLVector;
+ typedef CreatorInterface< ITL_BasePreconditioner< MTLVector > > ParallelPreconditionCreator;
+ typedef ITL_BasePreconditioner< MTLVector > ParallelPreconditioner;
+
+ template< typename ITLSolver >
+ class PITL_Solver : public AMDiS::MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_runner< ITLSolver > > {
+ public:
+ PITL_Solver(std::string name):
+ AMDiS::MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_runner< ITLSolver > >(name) {}
+
+ int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+ SystemVector& x,
+ SystemVector& b,
+ ParallelMapper& m)
+ {
+ return AMDiS::MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_runner< ITLSolver > >::solve(A,x,b,m);
+ }
+
+
+ /// Creator class used in the OEMSolverMap.
+ class Creator : public OEMSolverCreator
+ {
+ public:
+ virtual ~Creator() {}
+
+ /// Returns a new CGSolver object.
+ OEMSolver* create()
+ {
+ return new PITL_Solver<ITLSolver>(this->name);
+ }
+ };
+ };
+
+ template< typename ITLSolver >
+ class PITL_Solver_para : public AMDiS::MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > > {
+ public:
+ PITL_Solver_para(std::string name):
+ AMDiS::MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > >(name) {}
+
+ /// Creator class used in the OEMSolverMap.
+ class Creator : public OEMSolverCreator
+ {
+ public:
+ virtual ~Creator() {}
+
+ /// Returns a new CGSolver object.
+ OEMSolver* create()
+ {
+ return new PITL_Solver_para<ITLSolver>(this->name);
+ }
+ };
+ };
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by the conjugate gradient method (CG) and can be used for
+ * symmetric positive definite system matrices.
+ * Right preconditioner is ignored.
+ */
+ class CGSolver : public PITL_Solver<cg_solver_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ CGSolver(std::string name) : PITL_Solver<cg_solver_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by the squared conjugate gradient method (CGS).
+ * Right preconditioner is ignored.
+ */
+ class CGSSolver : public PITL_Solver<cgs_solver_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ CGSSolver(std::string name) : PITL_Solver<cgs_solver_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by a stabilized BiCG method and can be used for
+ * system matrices.
+ */
+ class BiCGSolver : public PITL_Solver<bicg_solver_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ BiCGSolver(std::string name) : PITL_Solver<bicg_solver_type>(name) {}
+ };
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by a stabilized BiCG method and can be used for
+ * system matrices.
+ */
+ class BiCGStabSolver : public PITL_Solver<bicgstab_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ BiCGStabSolver(std::string name) : PITL_Solver<bicgstab_type>(name) {}
+ };
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by a stabilized BiCG method and can be used for
+ * system matrices.
+ */
+ class BiCGStab2Solver : public PITL_Solver<bicgstab2_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ BiCGStab2Solver(std::string name) : PITL_Solver<bicgstab2_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by the Quasi-Minimal Residual method (QMR).
+ */
+ class QMRSolver : public PITL_Solver<qmr_solver_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ QMRSolver(std::string name) : PITL_Solver<qmr_solver_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by the Transposed-Free Quasi-Minimal Residual method (TFQMR).
+ * Does not use preconditioning currently.
+ */
+ class TFQMRSolver : public PITL_Solver<tfqmr_solver_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ TFQMRSolver(std::string name) : PITL_Solver<tfqmr_solver_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by a stabilized BiCG(l) method and can be used for
+ * system matrices.
+ */
+ class BiCGStabEllSolver : public PITL_Solver_para<bicgstab_ell_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ BiCGStabEllSolver(std::string name) : PITL_Solver_para<bicgstab_ell_type>(name) {}
+ };
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by the GMRES method.
+ * The parameter ell is the maximal number of orthogonalized vectors.
+ * The method is not preconditioned
+ */
+ class GMResSolver : public PITL_Solver_para<gmres_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ GMResSolver(std::string name) : PITL_Solver_para<gmres_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by an Induced Dimension Reduction method and can be used for
+ * system matrices.
+ * The parameter s can be specified as ell.
+ * Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems.
+ * SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
+ *
+ */
+ class IDRsSolver : public PITL_Solver_para<idr_s_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ IDRsSolver(std::string name) : PITL_Solver_para<idr_s_type>(name) {}
+ };
+
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Solves a linear system by the Minres method. Can be used for symmetric
+ * indefinite systems.
+ */
+ class MinResSolver : public PITL_Solver<minres_solver_type>
+ {
+ public:
+ /// The constructor reads required parameters and sets solvers \ref name.
+ MinResSolver(std::string name) : PITL_Solver<minres_solver_type>(name) {}
+ };
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Diagonal preconditioner.
+ */
+ class DiagonalPreconditioner
+ : public ITL_Preconditioner<itl::pc::diagonal<MTLMatrix >, MTLVector, MTLMatrix >
+ {
+ typedef ITL_Preconditioner<itl::pc::diagonal<MTLMatrix >, MTLVector, MTLMatrix > base;
+
+ public:
+ DiagonalPreconditioner(const MTLMatrix& A) : base(A) {}
+ };
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * Identity preconditioner. Behaves like no preconditioning.
+ */
+ class IdentityPreconditioner
+ : public ITL_Preconditioner<itl::pc::identity<MTLMatrix >, MTLVector, MTLMatrix >
+ {
+ typedef ITL_Preconditioner<itl::pc::identity<MTLMatrix >, MTLVector, MTLMatrix > base;
+
+ public:
+ IdentityPreconditioner(const MTLMatrix& A) : base(A) {}
+ };
+
+
+ // ============================================================================
+ // ===== class ILUPreconditioner ==============================================
+ // ============================================================================
+
+ /**
+ * \ingroup Solver
+ *
+ * \brief
+ * ILU (Incomplete LU factorization) preconditioner.
+ *
+ * The preconditioner is used from ITL. It corresponds for instance to "Iterative Methods for
+ * Sparce Linear Systems", second edition, Yousef Saad. The preconditioner is
+ * described in chapter 10.3 (algorithm 10.4).
+ */
+ class ILUPreconditioner
+ : public ITL_Preconditioner< itl::pc::ilu_0<MTLMatrix >,MTLVector, MTLMatrix >
+ {};
+
+ class ICPreconditioner
+ : public ITL_Preconditioner< itl::pc::ic_0<MTLMatrix >, MTLVector, MTLMatrix >
+ {};
+
+
+ }
+ template< >
+ struct Collection< Parallel::MTLMatrix > {
+ typedef mtl::matrix::inserter< Parallel::MTLMatrix , mtl::update_plus< Parallel::MTLMatrix::value_type > > Inserter;
+ };
+
+ template< >
+ struct Collection< Parallel::MTLVector > {
+ typedef mtl::vector::inserter< Parallel::MTLVector , mtl::update_plus< Parallel::MTLVector::value_type > > Inserter;
+ typedef Parallel::MTLMatrix PreconditionMatrix;
+ };
+
+} // namespace AMDiS
+
+#endif // AMDIS_PITL_SOLVER
+
diff --git a/AMDiS/src/parallel/ParallelMapper.h b/AMDiS/src/parallel/ParallelMapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..5b9c36e144d884e943a2ac88238bcae06428f42b
--- /dev/null
+++ b/AMDiS/src/parallel/ParallelMapper.h
@@ -0,0 +1,42 @@
+#ifndef PARALLEL_MAPPER_H
+#define PARALLEL_MAPPER_H
+
+#include "parallel/MeshDistributor.h"
+namespace AMDiS {
+ //todo: request from a problem -> distributor mapper?
+ template< typename size_type >
+ class ParallelMapper_base {
+ MeshDistributor& md;
+ /// the current row in the problem system
+ size_type r;
+ /// the current column in the problem system
+ size_type c;
+ int nComponents;
+
+ public:
+
+ ParallelMapper_base(MeshDistributor& md, int nComponents):
+ md(md),r(0),c(0),nComponents(nComponents)
+ {
+ TEST_EXIT_DBG(nComponents > 0)("the system must have at least one component");
+ }
+ size_type row(size_type r_)
+ {
+ size_type ret = md.mapLocalToGlobal(r_)*nComponents+r;
+ return ret;
+ }
+ size_type col(size_type c_)
+ {
+ size_type ret = md.mapLocalToGlobal(c_)*nComponents+c;
+ return ret ;
+ }
+
+ inline void setRow(size_type r_) { r = r_ ;}
+ inline void setCol(size_type c_) { c = c_; }
+ inline size_type nRow() { return md.getNumberOverallDofs() * nComponents; }
+ inline size_type nCol() { return md.getNumberOverallDofs() * nComponents; }
+ };
+
+ typedef ParallelMapper_base< unsigned int > ParallelMapper;
+}
+#endif