diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index 197e8ecc3d82b8182aaf83b87404acd268aaab47..5e5ae77862efc1febdda152e8073cd15c78d52e6 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -28,8 +28,9 @@ endif
if USE_PARALLEL_DOMAIN_AMDIS
PARALLEL_AMDIS_SOURCES += \
+ $(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
$(SOURCE_DIR)/parallel/ParallelDomainBase.h $(SOURCE_DIR)/parallel/ParallelDomainBase.cc \
- $(SOURCE_DIR)/parallel/ParallelDomainVec.h $(SOURCE_DIR)/parallel/ParallelDomainVec.cc \
+ $(SOURCE_DIR)/parallel/GlobalMatrixSolver.h $(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc \
$(SOURCE_DIR)/parallel/MpiHelper.h $(SOURCE_DIR)/parallel/MpiHelper.cc
libamdis_la_CXXFLAGS += -DHAVE_PARALLEL_DOMAIN_AMDIS=1
AMDIS_INCLUDES += -I$(PETSC_DIR)/include -I$(PETSC_DIR)/$(PETSC_ARCH)/include
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 1861189a7d751bb335ffb019c25065879554f4d8..94777480f03113c35bd995e01349a9af9ef0874f 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -38,8 +38,9 @@ build_triplet = @build@
host_triplet = @host@
@USE_PARALLEL_AMDIS_TRUE@am__append_1 = -DHAVE_PARALLEL_AMDIS=1
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_2 = \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/ParallelDomainBase.h $(SOURCE_DIR)/parallel/ParallelDomainBase.cc \
-@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/ParallelDomainVec.h $(SOURCE_DIR)/parallel/ParallelDomainVec.cc \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/GlobalMatrixSolver.h $(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/MpiHelper.h $(SOURCE_DIR)/parallel/MpiHelper.cc
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_3 = -DHAVE_PARALLEL_DOMAIN_AMDIS=1
@@ -73,11 +74,12 @@ am__installdirs = "$(DESTDIR)$(libdir)"
libLTLIBRARIES_INSTALL = $(INSTALL)
LTLIBRARIES = $(lib_LTLIBRARIES)
libamdis_la_LIBADD =
-am__libamdis_la_SOURCES_DIST = \
+am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
+ $(SOURCE_DIR)/parallel/StdMpi.cc \
$(SOURCE_DIR)/parallel/ParallelDomainBase.h \
$(SOURCE_DIR)/parallel/ParallelDomainBase.cc \
- $(SOURCE_DIR)/parallel/ParallelDomainVec.h \
- $(SOURCE_DIR)/parallel/ParallelDomainVec.cc \
+ $(SOURCE_DIR)/parallel/GlobalMatrixSolver.h \
+ $(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc \
$(SOURCE_DIR)/parallel/MpiHelper.h \
$(SOURCE_DIR)/parallel/MpiHelper.cc \
$(PARALLEL_DIR)/ConditionalEstimator.h \
@@ -231,8 +233,9 @@ am__libamdis_la_SOURCES_DIST = \
$(SOURCE_DIR)/parallel/InteriorBoundary.h \
$(SOURCE_DIR)/parallel/InteriorBoundary.cc \
$(SOURCE_DIR)/Debug.h $(SOURCE_DIR)/Debug.cc
-@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__objects_1 = libamdis_la-ParallelDomainBase.lo \
-@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelDomainVec.lo \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__objects_1 = libamdis_la-StdMpi.lo \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelDomainBase.lo \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-GlobalMatrixSolver.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-MpiHelper.lo
@USE_PARALLEL_AMDIS_FALSE@am__objects_2 = $(am__objects_1)
@USE_PARALLEL_AMDIS_TRUE@am__objects_2 = \
@@ -767,6 +770,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-FixVec.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-GNUPlotWriter.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Global.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-GlobalMatrixSolver.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-InteriorBoundary.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Lagrange.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-LeafData.Plo@am__quote@
@@ -782,7 +786,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Operator.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParMetisPartitioner.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParallelDomainBase.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParallelDomainVec.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParallelProblem.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Parameters.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Parametric.Plo@am__quote@
@@ -812,6 +815,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ScalableQuadrature.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-SecondOrderAssembler.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-StandardProblemIteration.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-StdMpi.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-SubAssembler.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-SubElInfo.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-SurfaceQuadrature.Plo@am__quote@
@@ -856,6 +860,13 @@ distclean-compile:
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LTCXXCOMPILE) -c -o $@ $<
+libamdis_la-StdMpi.lo: $(SOURCE_DIR)/parallel/StdMpi.cc
+@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-StdMpi.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-StdMpi.Tpo" -c -o libamdis_la-StdMpi.lo `test -f '$(SOURCE_DIR)/parallel/StdMpi.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/StdMpi.cc; \
+@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-StdMpi.Tpo" "$(DEPDIR)/libamdis_la-StdMpi.Plo"; else rm -f "$(DEPDIR)/libamdis_la-StdMpi.Tpo"; exit 1; fi
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/parallel/StdMpi.cc' object='libamdis_la-StdMpi.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-StdMpi.lo `test -f '$(SOURCE_DIR)/parallel/StdMpi.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/StdMpi.cc
+
libamdis_la-ParallelDomainBase.lo: $(SOURCE_DIR)/parallel/ParallelDomainBase.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ParallelDomainBase.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-ParallelDomainBase.Tpo" -c -o libamdis_la-ParallelDomainBase.lo `test -f '$(SOURCE_DIR)/parallel/ParallelDomainBase.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParallelDomainBase.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-ParallelDomainBase.Tpo" "$(DEPDIR)/libamdis_la-ParallelDomainBase.Plo"; else rm -f "$(DEPDIR)/libamdis_la-ParallelDomainBase.Tpo"; exit 1; fi
@@ -863,12 +874,12 @@ libamdis_la-ParallelDomainBase.lo: $(SOURCE_DIR)/parallel/ParallelDomainBase.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ParallelDomainBase.lo `test -f '$(SOURCE_DIR)/parallel/ParallelDomainBase.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParallelDomainBase.cc
-libamdis_la-ParallelDomainVec.lo: $(SOURCE_DIR)/parallel/ParallelDomainVec.cc
-@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ParallelDomainVec.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-ParallelDomainVec.Tpo" -c -o libamdis_la-ParallelDomainVec.lo `test -f '$(SOURCE_DIR)/parallel/ParallelDomainVec.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParallelDomainVec.cc; \
-@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-ParallelDomainVec.Tpo" "$(DEPDIR)/libamdis_la-ParallelDomainVec.Plo"; else rm -f "$(DEPDIR)/libamdis_la-ParallelDomainVec.Tpo"; exit 1; fi
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/parallel/ParallelDomainVec.cc' object='libamdis_la-ParallelDomainVec.lo' libtool=yes @AMDEPBACKSLASH@
+libamdis_la-GlobalMatrixSolver.lo: $(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc
+@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-GlobalMatrixSolver.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-GlobalMatrixSolver.Tpo" -c -o libamdis_la-GlobalMatrixSolver.lo `test -f '$(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc; \
+@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-GlobalMatrixSolver.Tpo" "$(DEPDIR)/libamdis_la-GlobalMatrixSolver.Plo"; else rm -f "$(DEPDIR)/libamdis_la-GlobalMatrixSolver.Tpo"; exit 1; fi
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc' object='libamdis_la-GlobalMatrixSolver.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ParallelDomainVec.lo `test -f '$(SOURCE_DIR)/parallel/ParallelDomainVec.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParallelDomainVec.cc
+@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-GlobalMatrixSolver.lo `test -f '$(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/GlobalMatrixSolver.cc
libamdis_la-MpiHelper.lo: $(SOURCE_DIR)/parallel/MpiHelper.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-MpiHelper.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-MpiHelper.Tpo" -c -o libamdis_la-MpiHelper.lo `test -f '$(SOURCE_DIR)/parallel/MpiHelper.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/MpiHelper.cc; \
diff --git a/AMDiS/libtool b/AMDiS/libtool
index b24ee23d4db5c89b8becf2ee49cfa31486adc673..d6ae6aec3bdfe931f7738cb747117b7a8a279992 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -30,10 +30,10 @@
# the same distribution terms that you use for the rest of that program.
# A sed program that does not truncate output.
-SED="/bin/sed"
+SED="/usr/bin/sed"
# Sed that helps us avoid accidentally triggering echo(1) options like -n.
-Xsed="/bin/sed -e 1s/^X//"
+Xsed="/usr/bin/sed -e 1s/^X//"
# The HP-UX ksh and POSIX shell print the target directory to stdout
# if CDPATH is set.
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host deimos102:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -66,12 +66,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -82,13 +82,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="gcc"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -97,7 +97,7 @@ with_gcc=yes
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -171,7 +171,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -325,10 +325,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib64/gcc/x86_64-suse-linux/4.1.2/ /usr/lib/gcc/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/ /lib/x86_64-suse-linux/4.1.2/ /lib/../lib64/ /usr/lib/x86_64-suse-linux/4.1.2/ /usr/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../ /lib/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host deimos102:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6782,12 +6782,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -6798,13 +6798,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="g++"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpiCC"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -6813,7 +6813,7 @@ with_gcc=yes
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -6887,7 +6887,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6942,11 +6942,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects="/usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crti.o /usr/lib/gcc/i386-redhat-linux/4.1.2/crtbeginS.o"
+predep_objects="/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/crti.o /usr/lib64/gcc/x86_64-suse-linux/4.1.2/crtbeginS.o"
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects="/usr/lib/gcc/i386-redhat-linux/4.1.2/crtendS.o /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crtn.o"
+postdep_objects="/usr/lib64/gcc/x86_64-suse-linux/4.1.2/crtendS.o /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/crtn.o"
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -6954,11 +6954,11 @@ predeps=""
# Dependencies to place after the objects being linked to create a
# shared library.
-postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
+postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -libverbs -lrt -lnuma -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path="-L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.."
+compiler_lib_search_path="-L/usr/lib64 -L/licsoft/libraries/openmpi/1.2.6/64bit/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7038,10 +7038,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib64/gcc/x86_64-suse-linux/4.1.2/ /usr/lib/gcc/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/ /lib/x86_64-suse-linux/4.1.2/ /lib/../lib64/ /usr/lib/x86_64-suse-linux/4.1.2/ /usr/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../ /lib/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host deimos102:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7087,12 +7087,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -7103,7 +7103,7 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
@@ -7112,13 +7112,13 @@ LTCFLAGS="-g -O2"
CC="g77"
# Is the compiler the GNU C compiler?
-with_gcc=yes
+with_gcc=
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -7346,10 +7346,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../ /lib/i386-redhat-linux/3.4.6/ /lib/ /usr/lib/i386-redhat-linux/3.4.6/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/ /usr/lib/gcc/x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../ /lib/x86_64-suse-linux/3.3.5/ /lib/ /usr/lib/x86_64-suse-linux/3.3.5/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 2b8a9cfe25496edf0af721e672a4863e23d32b06..3b4849c5a180493bf975ed58b35501060b1479bd 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -95,7 +95,7 @@
#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/InteriorBoundary.h"
#include "parallel/MpiHelper.h"
-#include "parallel/ParallelDomainVec.h"
+#include "parallel/GlobalMatrixSolver.h"
#endif
#endif
diff --git a/AMDiS/src/parallel/GlobalMatrixSolver.cc b/AMDiS/src/parallel/GlobalMatrixSolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..57838d521fab4bc210ce43aad14118465d8f9eb1
--- /dev/null
+++ b/AMDiS/src/parallel/GlobalMatrixSolver.cc
@@ -0,0 +1,488 @@
+#include "GlobalMatrixSolver.h"
+#include "DOFVector.h"
+#include "Debug.h"
+#include "SystemVector.h"
+#include "parallel/StdMpi.h"
+
+#include "petscksp.h"
+
+namespace AMDiS {
+
+ PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
+ {
+ if (iter % 10 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
+ std::cout << "[0] Petsc-Iteration " << iter << ": " << rnorm << std::endl;
+
+ return 0;
+ }
+
+ void GlobalMatrixSolver::solve()
+ {
+ FUNCNAME("GlobalMatrixSolver::solve()");
+
+#ifdef _OPENMP
+ double wtime = omp_get_wtime();
+#endif
+ clock_t first = clock();
+
+ fillPetscMatrix(probStat->getSystemMatrix(), probStat->getRHS());
+ solvePetscMatrix(*(probStat->getSolution()));
+
+#ifdef _OPENMP
+ INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
+ TIME_USED(first, clock()),
+ omp_get_wtime() - wtime);
+#else
+ INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
+ TIME_USED(first, clock()));
+#endif
+ }
+
+
+ void GlobalMatrixSolver::setDofMatrix(DOFMatrix* mat, int dispMult,
+ int dispAddRow, int dispAddCol)
+ {
+ FUNCNAME("GlobalMatrixSolver::setDofMatrix()");
+
+ TEST_EXIT(mat)("No DOFMatrix!\n");
+
+ using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
+ namespace traits= mtl::traits;
+ typedef DOFMatrix::base_matrix_type Matrix;
+
+ traits::col<Matrix>::type col(mat->getBaseMatrix());
+ traits::const_value<Matrix>::type value(mat->getBaseMatrix());
+
+ typedef traits::range_generator<row, Matrix>::type cursor_type;
+ typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+ std::vector<int> cols;
+ std::vector<double> values;
+ cols.reserve(300);
+ values.reserve(300);
+
+ // === Traverse all rows of the dof matrix and insert row wise the values ===
+ // === to the petsc matrix. ===
+
+ for (cursor_type cursor = begin<row>(mat->getBaseMatrix()),
+ cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {
+
+ cols.clear();
+ values.clear();
+
+ // Global index of the current row dof.
+ int globalRowDof = mapLocalGlobalDofs[*cursor];
+ // Test if the current row dof is a periodic dof.
+ bool periodicRow = (periodicDof.count(globalRowDof) > 0);
+
+
+ // === Traverse all non zero entries of the row and produce vector cols ===
+ // === with the column indices of all row entries and vector values ===
+ // === with the corresponding values. ===
+
+ for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
+ icursor != icend; ++icursor) {
+
+ // Set only non null values.
+ if (value(*icursor) != 0.0) {
+ // Global index of the current column index.
+ int globalColDof = mapLocalGlobalDofs[col(*icursor)];
+ // Calculate the exact position of the column index in the petsc matrix.
+ int colIndex = globalColDof * dispMult + dispAddCol;
+
+ // If the current row is not periodic, but the current dof index is periodic,
+ // we have to duplicate the value to the other corresponding periodic columns.
+ if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
+ // The value is assign to n matrix entries, therefore, every entry
+ // has only 1/n value of the original entry.
+ double scalFactor = 1.0 / (periodicDof[globalColDof].size() + 1.0);
+
+ // Insert original entry.
+ cols.push_back(colIndex);
+ values.push_back(value(*icursor) * scalFactor);
+
+ // Insert the periodic entries.
+ for (std::set<DegreeOfFreedom>::iterator it =
+ periodicDof[globalColDof].begin();
+ it != periodicDof[globalColDof].end(); ++it) {
+ cols.push_back(*it * dispMult + dispAddCol);
+ values.push_back(value(*icursor) * scalFactor);
+ }
+ } else {
+ // Neigher row nor column dof index is periodic, simple add entry.
+ cols.push_back(colIndex);
+ values.push_back(value(*icursor));
+ }
+ }
+ }
+
+
+ // === Up to now we have assembled on row. Now, the row must be send to the ===
+ // === corresponding rows to the petsc matrix. ===
+
+ // Calculate petsc row index.
+ int rowIndex = globalRowDof * dispMult + dispAddRow;
+
+ if (periodicRow) {
+ // The row dof is periodic, so send dof to all the corresponding rows.
+
+ double scalFactor = 1.0 / (periodicDof[globalRowDof].size() + 1.0);
+
+ int diagIndex = -1;
+ for (int i = 0; i < static_cast<int>(values.size()); i++) {
+ // Change only the non diagonal values in the col. For the diagonal test
+ // we compare the global dof indices of the dof matrix (not of the petsc
+ // matrix!).
+ if ((cols[i] - dispAddCol) / dispMult != globalRowDof)
+ values[i] *= scalFactor;
+ else
+ diagIndex = i;
+ }
+
+ // Send the main row to the petsc matrix.
+ MatSetValues(petscMatrix, 1, &rowIndex, cols.size(),
+ &(cols[0]), &(values[0]), ADD_VALUES);
+
+ // Set diagonal element to zero, i.e., the diagonal element of the current
+ // row is not send to the periodic row indices.
+ if (diagIndex != -1)
+ values[diagIndex] = 0.0;
+
+ // Send the row to all periodic row indices.
+ for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRowDof].begin();
+ it != periodicDof[globalRowDof].end(); ++it) {
+ int perRowIndex = *it * dispMult + dispAddRow;
+ MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(),
+ &(cols[0]), &(values[0]), ADD_VALUES);
+ }
+
+ } else {
+ // The row dof is not periodic, simply send the row to the petsc matrix.
+ MatSetValues(petscMatrix, 1, &rowIndex, cols.size(),
+ &(cols[0]), &(values[0]), ADD_VALUES);
+ }
+ }
+ }
+
+
+ void GlobalMatrixSolver::setDofVector(Vec& petscVec, DOFVector<double>* vec,
+ int dispMult, int dispAdd)
+ {
+ // Traverse all used dofs in the dof vector.
+ DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
+ for (dofIt.reset(); !dofIt.end(); ++dofIt) {
+ // Calculate global row index of the dof.
+ int globalRow = mapLocalGlobalDofs[dofIt.getDOFIndex()];
+ // Calculate petsc index of the row dof.
+ int index = globalRow * dispMult + dispAdd;
+
+ if (periodicDof.count(globalRow) > 0) {
+ // The dof index is periodic, so devide the value to all dof entries.
+
+ double value = *dofIt / (periodicDof[globalRow].size() + 1.0);
+ VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+
+ for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRow].begin();
+ it != periodicDof[globalRow].end(); ++it) {
+ index = *it * dispMult + dispAdd;
+ VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+ }
+
+ } else {
+ // The dof index is not periodic.
+ double value = *dofIt;
+ VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+ }
+ }
+ }
+
+
+ void GlobalMatrixSolver::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
+ {
+ FUNCNAME("GlobalMatrixSolver::createPetscNnzStructure()");
+
+ TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
+ TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
+
+ d_nnz = new int[nRankRows];
+ o_nnz = new int[nRankRows];
+ for (int i = 0; i < nRankRows; i++) {
+ d_nnz[i] = 0;
+ o_nnz[i] = 0;
+ }
+
+ using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
+ namespace traits = mtl::traits;
+ typedef DOFMatrix::base_matrix_type Matrix;
+ typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
+
+ // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
+ // that this rank will send to. This nnz entries will be assembled on this rank,
+ // but because the row DOFs are not DOFs of this rank they will be send to the
+ // owner of the row DOFs.
+ std::map<int, MatrixNnzEntry> sendMatrixEntry;
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ if ((*mat)[i][j]) {
+ Matrix bmat = (*mat)[i][j]->getBaseMatrix();
+
+ traits::col<Matrix>::type col(bmat);
+ traits::const_value<Matrix>::type value(bmat);
+
+ typedef traits::range_generator<row, Matrix>::type cursor_type;
+ typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+ for (cursor_type cursor = begin<row>(bmat),
+ cend = end<row>(bmat); cursor != cend; ++cursor) {
+
+ // Map the local row number to the global DOF index and create from it
+ // the global PETSc row index of this DOF.
+ int petscRowIdx = mapLocalGlobalDofs[*cursor] * nComponents + i;
+
+ if (isRankDof[*cursor]) {
+
+ // === The current row DOF is a rank dof, so create the corresponding ===
+ // === nnz values directly on rank's nnz data. ===
+
+ // This is the local row index of the local PETSc matrix.
+ int localPetscRowIdx = petscRowIdx - rstart * nComponents;
+
+#if (DEBUG != 0)
+ if (localPetscRowIdx < 0 || localPetscRowIdx >= nRankRows) {
+ std::cout << "ERROR in rank: " << mpiRank << std::endl;
+ std::cout << " Wrong r = " << localPetscRowIdx << " " << *cursor
+ << " " << mapLocalGlobalDofs[*cursor] << " "
+ << nRankRows << std::endl;
+ ERROR_EXIT("Should not happen!\n");
+ }
+#endif
+
+ // Traverse all non zero entries in this row.
+ for (icursor_type icursor = begin<nz>(cursor),
+ icend = end<nz>(cursor); icursor != icend; ++icursor) {
+ if (value(*icursor) != 0.0) {
+ int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
+
+ // The row DOF is a rank DOF, if also the column is a rank DOF,
+ // increment the d_nnz values for this row, otherwise the o_nnz value.
+ if (petscColIdx >= rstart * nComponents &&
+ petscColIdx < rstart * nComponents + nRankRows)
+ d_nnz[localPetscRowIdx]++;
+ else
+ o_nnz[localPetscRowIdx]++;
+ }
+ }
+ } else {
+
+ // === The current row DOF is not a rank dof, i.e., it will be created ===
+ // === on this rank, but after this it will be send to another rank ===
+ // === matrix. So we need to send also the corresponding nnz structure ===
+ // === of this row to the corresponding rank. ===
+
+ // Find out who is the member of this DOF.
+ int sendToRank = -1;
+ for (RankToDofContainer::iterator it = recvDofs.begin();
+ it != recvDofs.end(); ++it) {
+ for (DofContainer::iterator dofIt = it->second.begin();
+ dofIt != it->second.end(); ++dofIt) {
+ if (**dofIt == *cursor) {
+
+ if (petscRowIdx == 6717) {
+ debug::writeDofMesh(mpiRank, *cursor, feSpace);
+ MSG("SEND DOF TO: %d/%d\n", it->first, *cursor);
+ }
+
+ sendToRank = it->first;
+ break;
+ }
+ }
+
+ if (sendToRank != -1)
+ break;
+ }
+
+ TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n");
+
+ // Send all non zero entries to the member of the row DOF.
+ for (icursor_type icursor = begin<nz>(cursor),
+ icend = end<nz>(cursor); icursor != icend; ++icursor) {
+ if (value(*icursor) != 0.0) {
+ int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
+
+ sendMatrixEntry[sendToRank].
+ push_back(std::make_pair(petscRowIdx, petscColIdx));
+ }
+ }
+
+ } // if (isRankDof[*cursor]) ... else ...
+ } // for each row in mat[i][j]
+ } // if mat[i][j]
+ }
+ }
+
+ // === Send and recv the nnz row structure to/from other ranks. ===
+
+ StdMpi<MatrixNnzEntry> stdMpi(mpiComm, true);
+ stdMpi.send(sendMatrixEntry);
+ stdMpi.recv(sendDofs);
+ stdMpi.startCommunication<int>(MPI_INT);
+
+ // === Evaluate the nnz structure this rank got from other ranks and add it to ===
+ // === the PETSc nnz data structure. ===
+
+ for (std::map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
+ it != stdMpi.getRecvData().end(); ++it) {
+ if (it->second.size() > 0) {
+ for (unsigned int i = 0; i < it->second.size(); i++) {
+ int r = it->second[i].first;
+ int c = it->second[i].second;
+
+ int localRowIdx = r - rstart * nComponents;
+
+ TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
+ ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
+ r, localRowIdx, nRankRows, it->first);
+
+ if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
+ o_nnz[localRowIdx]++;
+ else
+ d_nnz[localRowIdx]++;
+ }
+ }
+ }
+ }
+
+
+ void GlobalMatrixSolver::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
+ {
+ FUNCNAME("GlobalMatrixSolver::fillPetscMatrix()");
+
+ clock_t first = clock();
+
+ // === Create PETSc vector (rhs, solution and a temporary vector). ===
+
+ VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
+ VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
+ VecSetType(petscRhsVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscSolVec);
+ VecSetSizes(petscSolVec, nRankRows, nOverallRows);
+ VecSetType(petscSolVec, VECMPI);
+
+ VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
+ VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
+ VecSetType(petscTmpVec, VECMPI);
+
+ if (!d_nnz)
+ createPetscNnzStructure(mat);
+
+ // === Create PETSc matrix with the computed nnz data structure. ===
+
+ MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
+ 0, d_nnz, 0, o_nnz, &petscMatrix);
+
+ INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", TIME_USED(first, clock()));
+
+#if (DEBUG != 0)
+ int a, b;
+ MatGetOwnershipRange(petscMatrix, &a, &b);
+ TEST_EXIT(a == rstart * nComponents)("Wrong matrix ownership range!\n");
+ TEST_EXIT(b == rstart * nComponents + nRankRows)("Wrong matrix ownership range!\n");
+#endif
+
+ // === Transfer values from DOF matrices to the PETSc matrix. ===
+
+ for (int i = 0; i < nComponents; i++)
+ for (int j = 0; j < nComponents; j++)
+ if ((*mat)[i][j])
+ setDofMatrix((*mat)[i][j], nComponents, i, j);
+
+ INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", TIME_USED(first, clock()));
+
+ MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
+
+ // === Transfer values from DOF vector to the PETSc vector. ===
+
+ for (int i = 0; i < nComponents; i++)
+ setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
+
+ VecAssemblyBegin(petscRhsVec);
+ VecAssemblyEnd(petscRhsVec);
+
+ INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
+ }
+
+
+ void GlobalMatrixSolver::solvePetscMatrix(SystemVector &vec)
+ {
+ FUNCNAME("GlobalMatrixSolver::solvePetscMatrix()");
+
+#if 0
+ // Set old solution to be initiual guess for petsc solver.
+ for (int i = 0; i < nComponents; i++)
+ setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i);
+
+ VecAssemblyBegin(petscSolVec);
+ VecAssemblyEnd(petscSolVec);
+#endif
+
+ // === Init Petsc solver. ===
+
+ KSP solver;
+ KSPCreate(PETSC_COMM_WORLD, &solver);
+ KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN);
+ KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
+ KSPSetType(solver, KSPBCGS);
+ KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
+ KSPSetFromOptions(solver);
+ // Do not delete the solution vector, use it for the initial guess.
+ // KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
+
+
+ // === Run Petsc. ===
+
+ KSPSolve(solver, petscRhsVec, petscSolVec);
+
+ // === Transfere values from Petsc's solution vectors to the dof vectors.
+ PetscScalar *vecPointer;
+ VecGetArray(petscSolVec, &vecPointer);
+
+ for (int i = 0; i < nComponents; i++) {
+ DOFVector<double> *dofvec = vec.getDOFVector(i);
+ for (int j = 0; j < nRankDofs; j++)
+ (*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];
+ }
+
+ VecRestoreArray(petscSolVec, &vecPointer);
+
+
+ // === Synchronize dofs at common dofs, i.e., dofs that correspond to more ===
+ // === than one partition. ===
+ synchVector(vec);
+
+
+ // === Print information about solution process. ===
+
+ int iterations = 0;
+ KSPGetIterationNumber(solver, &iterations);
+ MSG(" Number of iterations: %d\n", iterations);
+
+ double norm = 0.0;
+ MatMult(petscMatrix, petscSolVec, petscTmpVec);
+ VecAXPY(petscTmpVec, -1.0, petscRhsVec);
+ VecNorm(petscTmpVec, NORM_2, &norm);
+ MSG(" Residual norm: %e\n", norm);
+
+
+ // === Destroy Petsc's variables. ===
+
+ MatDestroy(petscMatrix);
+ VecDestroy(petscRhsVec);
+ VecDestroy(petscSolVec);
+ VecDestroy(petscTmpVec);
+ KSPDestroy(solver);
+ }
+
+}
diff --git a/AMDiS/src/parallel/GlobalMatrixSolver.h b/AMDiS/src/parallel/GlobalMatrixSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..c5087f46f36111829b05660366b83dc5ae1aac1e
--- /dev/null
+++ b/AMDiS/src/parallel/GlobalMatrixSolver.h
@@ -0,0 +1,93 @@
+// ============================================================================
+// == ==
+// == 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 GlobalMatrixSolver.h */
+
+#ifndef AMDIS_GLOBALMATRIXSOLVER_H
+#define AMDIS_GLOBALMATRIXSOLVER_H
+
+#include "AMDiS_fwd.h"
+#include "Global.h"
+#include "ParallelDomainBase.h"
+#include "ProblemVec.h"
+#include "ProblemInstat.h"
+
+#include "petsc.h"
+#include "petscsys.h"
+#include "petscao.h"
+
+namespace AMDiS {
+
+ class GlobalMatrixSolver : public ParallelDomainBase
+ {
+ public:
+ GlobalMatrixSolver(ProblemVec *problemStat, ProblemInstatVec *problemInstat)
+ : ParallelDomainBase(problemStat, problemInstat),
+ d_nnz(NULL),
+ o_nnz(NULL)
+ {}
+
+ ~GlobalMatrixSolver()
+ {}
+
+ void solve();
+
+ protected:
+ /// Creates a new non zero pattern structure for the Petsc matrix.
+ void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
+
+ /** \brief
+ * Create a PETSc matrix and PETSc vectors. The given DOF matrices are used to
+ * create the nnz structure of the PETSc matrix and the values are transfered to it.
+ * The given DOF vectors are used to the the values of the PETSc rhs vector.
+ *
+ * \param[in] mat
+ * \param[in] vec
+ */
+ void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
+
+ /// Takes a dof matrix and sends the values to the global petsc matrix.
+ void setDofMatrix(DOFMatrix* mat, int dispMult = 1,
+ int dispAddRow = 0, int dispAddCol = 0);
+
+ /// Takes a dof vector and sends its values to a given petsc vector.
+ void setDofVector(Vec& petscVec, DOFVector<double>* vec,
+ int disMult = 1, int dispAdd = 0);
+
+ void solvePetscMatrix(SystemVector &vec);
+
+ protected:
+ /// Petsc's matrix structure.
+ Mat petscMatrix;
+
+ /** \brief
+ * Petsc's vector structures for the rhs vector, the solution vector and a
+ * temporary vector for calculating the final residuum.
+ */
+ Vec petscRhsVec, petscSolVec, petscTmpVec;
+
+ /// Arrays definig the non zero pattern of Petsc's matrix.
+ int *d_nnz, *o_nnz;
+ };
+
+ typedef GlobalMatrixSolver ParallelDomainVec;
+
+} //namespace AMDiS
+
+#endif
diff --git a/AMDiS/src/parallel/ParallelDomainBase.cc b/AMDiS/src/parallel/ParallelDomainBase.cc
index 3057c509cf4e42c0b53f7ff82324b99b266f86e2..0f3697fafba27e318a83bab3c7ddf6c741076bd2 100644
--- a/AMDiS/src/parallel/ParallelDomainBase.cc
+++ b/AMDiS/src/parallel/ParallelDomainBase.cc
@@ -25,20 +25,10 @@
#include "ProblemInstat.h"
#include "Debug.h"
-#include "petscksp.h"
-
namespace AMDiS {
using boost::lexical_cast;
- PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
- {
- if (iter % 10 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
- std::cout << "[0] Petsc-Iteration " << iter << ": " << rnorm << std::endl;
-
- return 0;
- }
-
inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
{
return (*dof1 < *dof2);
@@ -55,8 +45,6 @@ namespace AMDiS {
refineManager(problemStat->getRefinementManager(0)),
info(problemStat->getInfo()),
initialPartitionMesh(true),
- d_nnz(NULL),
- o_nnz(NULL),
nRankDofs(0),
rstart(0),
nComponents(problemStat->getNumComponents()),
@@ -280,453 +268,6 @@ namespace AMDiS {
}
- void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult,
- int dispAddRow, int dispAddCol)
- {
- FUNCNAME("ParallelDomainBase::setDofMatrix()");
-
- TEST_EXIT(mat)("No DOFMatrix!\n");
-
- using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
- namespace traits= mtl::traits;
- typedef DOFMatrix::base_matrix_type Matrix;
-
- traits::col<Matrix>::type col(mat->getBaseMatrix());
- traits::const_value<Matrix>::type value(mat->getBaseMatrix());
-
- typedef traits::range_generator<row, Matrix>::type cursor_type;
- typedef traits::range_generator<nz, cursor_type>::type icursor_type;
-
- std::vector<int> cols;
- std::vector<double> values;
- cols.reserve(300);
- values.reserve(300);
-
- // === Traverse all rows of the dof matrix and insert row wise the values ===
- // === to the petsc matrix. ===
-
- for (cursor_type cursor = begin<row>(mat->getBaseMatrix()),
- cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {
-
- cols.clear();
- values.clear();
-
- // Global index of the current row dof.
- int globalRowDof = mapLocalGlobalDofs[*cursor];
- // Test if the current row dof is a periodic dof.
- bool periodicRow = (periodicDof.count(globalRowDof) > 0);
-
-
- // === Traverse all non zero entries of the row and produce vector cols ===
- // === with the column indices of all row entries and vector values ===
- // === with the corresponding values. ===
-
- for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
- icursor != icend; ++icursor) {
-
- // Set only non null values.
- if (value(*icursor) != 0.0) {
- // Global index of the current column index.
- int globalColDof = mapLocalGlobalDofs[col(*icursor)];
- // Calculate the exact position of the column index in the petsc matrix.
- int colIndex = globalColDof * dispMult + dispAddCol;
-
- // If the current row is not periodic, but the current dof index is periodic,
- // we have to duplicate the value to the other corresponding periodic columns.
- if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
- // The value is assign to n matrix entries, therefore, every entry
- // has only 1/n value of the original entry.
- double scalFactor = 1.0 / (periodicDof[globalColDof].size() + 1.0);
-
- // Insert original entry.
- cols.push_back(colIndex);
- values.push_back(value(*icursor) * scalFactor);
-
- // Insert the periodic entries.
- for (std::set<DegreeOfFreedom>::iterator it =
- periodicDof[globalColDof].begin();
- it != periodicDof[globalColDof].end(); ++it) {
- cols.push_back(*it * dispMult + dispAddCol);
- values.push_back(value(*icursor) * scalFactor);
- }
- } else {
- // Neigher row nor column dof index is periodic, simple add entry.
- cols.push_back(colIndex);
- values.push_back(value(*icursor));
- }
- }
- }
-
-
- // === Up to now we have assembled on row. Now, the row must be send to the ===
- // === corresponding rows to the petsc matrix. ===
-
- // Calculate petsc row index.
- int rowIndex = globalRowDof * dispMult + dispAddRow;
-
- if (periodicRow) {
- // The row dof is periodic, so send dof to all the corresponding rows.
-
- double scalFactor = 1.0 / (periodicDof[globalRowDof].size() + 1.0);
-
- int diagIndex = -1;
- for (int i = 0; i < static_cast<int>(values.size()); i++) {
- // Change only the non diagonal values in the col. For the diagonal test
- // we compare the global dof indices of the dof matrix (not of the petsc
- // matrix!).
- if ((cols[i] - dispAddCol) / dispMult != globalRowDof)
- values[i] *= scalFactor;
- else
- diagIndex = i;
- }
-
- // Send the main row to the petsc matrix.
- MatSetValues(petscMatrix, 1, &rowIndex, cols.size(),
- &(cols[0]), &(values[0]), ADD_VALUES);
-
- // Set diagonal element to zero, i.e., the diagonal element of the current
- // row is not send to the periodic row indices.
- if (diagIndex != -1)
- values[diagIndex] = 0.0;
-
- // Send the row to all periodic row indices.
- for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRowDof].begin();
- it != periodicDof[globalRowDof].end(); ++it) {
- int perRowIndex = *it * dispMult + dispAddRow;
- MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(),
- &(cols[0]), &(values[0]), ADD_VALUES);
- }
-
- } else {
- // The row dof is not periodic, simply send the row to the petsc matrix.
- MatSetValues(petscMatrix, 1, &rowIndex, cols.size(),
- &(cols[0]), &(values[0]), ADD_VALUES);
- }
- }
- }
-
-
- void ParallelDomainBase::setDofVector(Vec& petscVec, DOFVector<double>* vec,
- int dispMult, int dispAdd)
- {
- // Traverse all used dofs in the dof vector.
- DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
- for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- // Calculate global row index of the dof.
- int globalRow = mapLocalGlobalDofs[dofIt.getDOFIndex()];
- // Calculate petsc index of the row dof.
- int index = globalRow * dispMult + dispAdd;
-
- if (periodicDof.count(globalRow) > 0) {
- // The dof index is periodic, so devide the value to all dof entries.
-
- double value = *dofIt / (periodicDof[globalRow].size() + 1.0);
- VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
-
- for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRow].begin();
- it != periodicDof[globalRow].end(); ++it) {
- index = *it * dispMult + dispAdd;
- VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
- }
-
- } else {
- // The dof index is not periodic.
- double value = *dofIt;
- VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
- }
- }
- }
-
-
- void ParallelDomainBase::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
- {
- FUNCNAME("ParallelDomainBase::createPetscNnzStructure()");
-
- TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
- TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");
-
- d_nnz = new int[nRankRows];
- o_nnz = new int[nRankRows];
- for (int i = 0; i < nRankRows; i++) {
- d_nnz[i] = 0;
- o_nnz[i] = 0;
- }
-
- using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
- namespace traits = mtl::traits;
- typedef DOFMatrix::base_matrix_type Matrix;
- typedef std::vector<std::pair<int, int> > MatrixNnzEntry;
-
- // Stores to each rank a list of nnz entries (i.e. pairs of row and column index)
- // that this rank will send to. This nnz entries will be assembled on this rank,
- // but because the row DOFs are not DOFs of this rank they will be send to the
- // owner of the row DOFs.
- std::map<int, MatrixNnzEntry> sendMatrixEntry;
-
- for (int i = 0; i < nComponents; i++) {
- for (int j = 0; j < nComponents; j++) {
- if ((*mat)[i][j]) {
- Matrix bmat = (*mat)[i][j]->getBaseMatrix();
-
- traits::col<Matrix>::type col(bmat);
- traits::const_value<Matrix>::type value(bmat);
-
- typedef traits::range_generator<row, Matrix>::type cursor_type;
- typedef traits::range_generator<nz, cursor_type>::type icursor_type;
-
- for (cursor_type cursor = begin<row>(bmat),
- cend = end<row>(bmat); cursor != cend; ++cursor) {
-
- // Map the local row number to the global DOF index and create from it
- // the global PETSc row index of this DOF.
- int petscRowIdx = mapLocalGlobalDofs[*cursor] * nComponents + i;
-
- if (isRankDof[*cursor]) {
-
- // === The current row DOF is a rank dof, so create the corresponding ===
- // === nnz values directly on rank's nnz data. ===
-
- // This is the local row index of the local PETSc matrix.
- int localPetscRowIdx = petscRowIdx - rstart * nComponents;
-
-#if (DEBUG != 0)
- if (localPetscRowIdx < 0 || localPetscRowIdx >= nRankRows) {
- std::cout << "ERROR in rank: " << mpiRank << std::endl;
- std::cout << " Wrong r = " << localPetscRowIdx << " " << *cursor
- << " " << mapLocalGlobalDofs[*cursor] << " "
- << nRankRows << std::endl;
- ERROR_EXIT("Should not happen!\n");
- }
-#endif
-
- // Traverse all non zero entries in this row.
- for (icursor_type icursor = begin<nz>(cursor),
- icend = end<nz>(cursor); icursor != icend; ++icursor) {
- if (value(*icursor) != 0.0) {
- int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
-
- // The row DOF is a rank DOF, if also the column is a rank DOF,
- // increment the d_nnz values for this row, otherwise the o_nnz value.
- if (petscColIdx >= rstart * nComponents &&
- petscColIdx < rstart * nComponents + nRankRows)
- d_nnz[localPetscRowIdx]++;
- else
- o_nnz[localPetscRowIdx]++;
- }
- }
- } else {
-
- // === The current row DOF is not a rank dof, i.e., it will be created ===
- // === on this rank, but after this it will be send to another rank ===
- // === matrix. So we need to send also the corresponding nnz structure ===
- // === of this row to the corresponding rank. ===
-
- // Find out who is the member of this DOF.
- int sendToRank = -1;
- for (RankToDofContainer::iterator it = recvDofs.begin();
- it != recvDofs.end(); ++it) {
- for (DofContainer::iterator dofIt = it->second.begin();
- dofIt != it->second.end(); ++dofIt) {
- if (**dofIt == *cursor) {
-
- if (petscRowIdx == 6717) {
- debug::writeDofMesh(mpiRank, *cursor, feSpace);
- MSG("SEND DOF TO: %d/%d\n", it->first, *cursor);
- }
-
- sendToRank = it->first;
- break;
- }
- }
-
- if (sendToRank != -1)
- break;
- }
-
- TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n");
-
- // Send all non zero entries to the member of the row DOF.
- for (icursor_type icursor = begin<nz>(cursor),
- icend = end<nz>(cursor); icursor != icend; ++icursor) {
- if (value(*icursor) != 0.0) {
- int petscColIdx = mapLocalGlobalDofs[col(*icursor)] * nComponents + j;
-
- sendMatrixEntry[sendToRank].
- push_back(std::make_pair(petscRowIdx, petscColIdx));
- }
- }
-
- } // if (isRankDof[*cursor]) ... else ...
- } // for each row in mat[i][j]
- } // if mat[i][j]
- }
- }
-
- // === Send and recv the nnz row structure to/from other ranks. ===
-
- StdMpi<MatrixNnzEntry> stdMpi(mpiComm, true);
- stdMpi.send(sendMatrixEntry);
- stdMpi.recv(sendDofs);
- stdMpi.startCommunication<int>(MPI_INT);
-
- // === Evaluate the nnz structure this rank got from other ranks and add it to ===
- // === the PETSc nnz data structure. ===
-
- for (std::map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
- it != stdMpi.getRecvData().end(); ++it) {
- if (it->second.size() > 0) {
- for (unsigned int i = 0; i < it->second.size(); i++) {
- int r = it->second[i].first;
- int c = it->second[i].second;
-
- int localRowIdx = r - rstart * nComponents;
-
- TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
- ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
- r, localRowIdx, nRankRows, it->first);
-
- if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
- o_nnz[localRowIdx]++;
- else
- d_nnz[localRowIdx]++;
- }
- }
- }
- }
-
-
- void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
- {
- FUNCNAME("ParallelDomainBase::fillPetscMatrix()");
-
- clock_t first = clock();
-
- // === Create PETSc vector (rhs, solution and a temporary vector). ===
-
- VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
- VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
- VecSetType(petscRhsVec, VECMPI);
-
- VecCreate(PETSC_COMM_WORLD, &petscSolVec);
- VecSetSizes(petscSolVec, nRankRows, nOverallRows);
- VecSetType(petscSolVec, VECMPI);
-
- VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
- VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
- VecSetType(petscTmpVec, VECMPI);
-
- if (!d_nnz)
- createPetscNnzStructure(mat);
-
- // === Create PETSc matrix with the computed nnz data structure. ===
-
- MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
- 0, d_nnz, 0, o_nnz, &petscMatrix);
-
- INFO(info, 8)("Fill petsc matrix 1 needed %.5f seconds\n", TIME_USED(first, clock()));
-
-#if (DEBUG != 0)
- int a, b;
- MatGetOwnershipRange(petscMatrix, &a, &b);
- TEST_EXIT(a == rstart * nComponents)("Wrong matrix ownership range!\n");
- TEST_EXIT(b == rstart * nComponents + nRankRows)("Wrong matrix ownership range!\n");
-#endif
-
- // === Transfer values from DOF matrices to the PETSc matrix. ===
-
- for (int i = 0; i < nComponents; i++)
- for (int j = 0; j < nComponents; j++)
- if ((*mat)[i][j])
- setDofMatrix((*mat)[i][j], nComponents, i, j);
-
- INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", TIME_USED(first, clock()));
-
- MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
-
- // === Transfer values from DOF vector to the PETSc vector. ===
-
- for (int i = 0; i < nComponents; i++)
- setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
-
- VecAssemblyBegin(petscRhsVec);
- VecAssemblyEnd(petscRhsVec);
-
- INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
- }
-
-
- void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
- {
- FUNCNAME("ParallelDomainBase::solvePetscMatrix()");
-
-#if 0
- // Set old solution to be initiual guess for petsc solver.
- for (int i = 0; i < nComponents; i++)
- setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i);
-
- VecAssemblyBegin(petscSolVec);
- VecAssemblyEnd(petscSolVec);
-#endif
-
- // === Init Petsc solver. ===
-
- KSP solver;
- KSPCreate(PETSC_COMM_WORLD, &solver);
- KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN);
- KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
- KSPSetType(solver, KSPBCGS);
- KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
- KSPSetFromOptions(solver);
- // Do not delete the solution vector, use it for the initial guess.
- // KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
-
-
- // === Run Petsc. ===
-
- KSPSolve(solver, petscRhsVec, petscSolVec);
-
- // === Transfere values from Petsc's solution vectors to the dof vectors.
- PetscScalar *vecPointer;
- VecGetArray(petscSolVec, &vecPointer);
-
- for (int i = 0; i < nComponents; i++) {
- DOFVector<double> *dofvec = vec.getDOFVector(i);
- for (int j = 0; j < nRankDofs; j++)
- (*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];
- }
-
- VecRestoreArray(petscSolVec, &vecPointer);
-
-
- // === Synchronize dofs at common dofs, i.e., dofs that correspond to more ===
- // === than one partition. ===
- synchVector(vec);
-
-
- // === Print information about solution process. ===
-
- int iterations = 0;
- KSPGetIterationNumber(solver, &iterations);
- MSG(" Number of iterations: %d\n", iterations);
-
- double norm = 0.0;
- MatMult(petscMatrix, petscSolVec, petscTmpVec);
- VecAXPY(petscTmpVec, -1.0, petscRhsVec);
- VecNorm(petscTmpVec, NORM_2, &norm);
- MSG(" Residual norm: %e\n", norm);
-
-
- // === Destroy Petsc's variables. ===
-
- MatDestroy(petscMatrix);
- VecDestroy(petscRhsVec);
- VecDestroy(petscSolVec);
- VecDestroy(petscTmpVec);
- KSPDestroy(solver);
- }
-
-
void ParallelDomainBase::synchVector(DOFVector<double> &vec)
{
StdMpi<std::vector<double> > stdMpi(mpiComm);
@@ -868,6 +409,7 @@ namespace AMDiS {
// === If there is a non zero pattern computed for Petsc matrix, delete it. So ===
// === it will be recomputed after next assembling. ===
+#if 0
if (d_nnz) {
delete [] d_nnz;
d_nnz = NULL;
@@ -877,6 +419,7 @@ namespace AMDiS {
delete [] o_nnz;
o_nnz = NULL;
}
+#endif
}
@@ -2354,29 +1897,6 @@ namespace AMDiS {
}
- void ParallelDomainBase::solve()
- {
- FUNCNAME("ParallelDomainBase::solve()");
-
-#ifdef _OPENMP
- double wtime = omp_get_wtime();
-#endif
- clock_t first = clock();
-
- fillPetscMatrix(probStat->getSystemMatrix(), probStat->getRHS());
- solvePetscMatrix(*(probStat->getSolution()));
-
-#ifdef _OPENMP
- INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
- TIME_USED(first, clock()),
- omp_get_wtime() - wtime);
-#else
- INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
- TIME_USED(first, clock()));
-#endif
- }
-
-
void ParallelDomainBase::serialize(std::ostream &out)
{
SerUtil::serialize(out, elemWeights);
diff --git a/AMDiS/src/parallel/ParallelDomainBase.h b/AMDiS/src/parallel/ParallelDomainBase.h
index a3136f9fb4313977f30cf10caf36604239f71ed6..e8fbd32c67fd110548c5053d311b6f09508d4730 100644
--- a/AMDiS/src/parallel/ParallelDomainBase.h
+++ b/AMDiS/src/parallel/ParallelDomainBase.h
@@ -38,10 +38,6 @@
#include "BoundaryManager.h"
#include "AMDiS_fwd.h"
-#include "petsc.h"
-#include "petscsys.h"
-#include "petscao.h"
-
namespace AMDiS {
struct DofPtrSortFct {
@@ -56,7 +52,7 @@ namespace AMDiS {
class ParallelDomainBase : public ProblemIterationInterface,
public ProblemTimeInterface
{
- private:
+ protected:
/// Defines a mapping type from DOFs to rank numbers.
typedef std::map<const DegreeOfFreedom*, int> DofToRank;
@@ -163,7 +159,7 @@ namespace AMDiS {
iterationIF->endIteration(adaptInfo);
}
- virtual void solve();
+ virtual void solve() = 0;
virtual int getNumProblems()
{
@@ -181,8 +177,6 @@ namespace AMDiS {
return nRankDofs;
}
- void solvePetscMatrix(SystemVector &vec);
-
virtual ProblemStatBase *getProblem(int number = 0)
{
return NULL;
@@ -267,27 +261,6 @@ namespace AMDiS {
DofToRank& boundaryDofs,
DofToBool& vertexDof);
- /// Creates a new non zero pattern structure for the Petsc matrix.
- void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
-
- /** \brief
- * Create a PETSc matrix and PETSc vectors. The given DOF matrices are used to
- * create the nnz structure of the PETSc matrix and the values are transfered to it.
- * The given DOF vectors are used to the the values of the PETSc rhs vector.
- *
- * \param[in] mat
- * \param[in] vec
- */
- void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
-
- /// Takes a dof matrix and sends the values to the global petsc matrix.
- void setDofMatrix(DOFMatrix* mat, int dispMult = 1,
- int dispAddRow = 0, int dispAddCol = 0);
-
- /// Takes a dof vector and sends its values to a given petsc vector.
- void setDofVector(Vec& petscVec, DOFVector<double>* vec,
- int disMult = 1, int dispAdd = 0);
-
/** \brief
* Checks for all given interior boundaries if the elements fit together on both
* sides of the boundaries. If this is not the case, the mesh is adapted. Because
@@ -541,19 +514,7 @@ namespace AMDiS {
* of the parition it corresponds to is stored.
*/
std::map<int, int> oldPartitionVec;
-
- /// Petsc's matrix structure.
- Mat petscMatrix;
-
- /** \brief
- * Petsc's vector structures for the rhs vector, the solution vector and a
- * temporary vector for calculating the final residuum.
- */
- Vec petscRhsVec, petscSolVec, petscTmpVec;
-
- /// Arrays definig the non zero pattern of Petsc's matrix.
- int *d_nnz, *o_nnz;
-
+
/// Number of DOFs in the rank mesh.
int nRankDofs;
diff --git a/AMDiS/src/parallel/ParallelDomainVec.cc b/AMDiS/src/parallel/ParallelDomainVec.cc
deleted file mode 100644
index 3ba7af68e180f5f13b4aee496d66800e783c3ba4..0000000000000000000000000000000000000000
--- a/AMDiS/src/parallel/ParallelDomainVec.cc
+++ /dev/null
@@ -1,20 +0,0 @@
-#include "boost/lexical_cast.hpp"
-#include "ParallelDomainVec.h"
-#include "ProblemVec.h"
-#include "ProblemInstat.h"
-#include "SystemVector.h"
-#include "Serializer.h"
-#include "BoundaryManager.h"
-
-namespace AMDiS {
-
- using boost::lexical_cast;
-
- ParallelDomainVec::ParallelDomainVec(ProblemVec *problem,
- ProblemInstatVec *problemInstat)
- : ParallelDomainBase(problem, problemInstat)
- {
- FUNCNAME("ParallelDomainVec::ParallelDomainVec()");
-
- }
-}
diff --git a/AMDiS/src/parallel/ParallelDomainVec.h b/AMDiS/src/parallel/ParallelDomainVec.h
deleted file mode 100644
index c234c07990c1d36a788b57cbb714642b0b43a9e8..0000000000000000000000000000000000000000
--- a/AMDiS/src/parallel/ParallelDomainVec.h
+++ /dev/null
@@ -1,40 +0,0 @@
-// ============================================================================
-// == ==
-// == 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 ParallelDomainVec.h */
-
-#ifndef AMDIS_PARALLELDOMAINVEC_H
-#define AMDIS_PARALLELDOMAINVEC_H
-
-#include "ParallelDomainBase.h"
-#include "ProblemVec.h"
-#include "BoundaryManager.h"
-
-namespace AMDiS {
-
- class ParallelDomainVec : public ParallelDomainBase
- {
- public:
- ParallelDomainVec(ProblemVec *problem,
- ProblemInstatVec *problemInstat);
- };
-
-}
-
-#endif // AMDIS_PARALLELDOMAINVEC_H
diff --git a/AMDiS/src/parallel/StdMpi.cc b/AMDiS/src/parallel/StdMpi.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6e1e525fb352bbfb4592ca3452bf6c7f9c3e65f9
--- /dev/null
+++ b/AMDiS/src/parallel/StdMpi.cc
@@ -0,0 +1,220 @@
+#include "StdMpi.h"
+
+namespace AMDiS {
+
+ int intSizeOf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data)
+ {
+ return data.size() * 2;
+ }
+
+ int intSizeOf(std::vector<MeshStructure> &data)
+ {
+ int s = 0;
+ for (unsigned int i = 0; i < data.size(); i++)
+ s += data[i].getCode().size() + 2;
+
+ return s;
+ }
+
+ int intSizeOf(std::vector<WorldVector<double> > &data)
+ {
+ return data.size() * Global::getGeo(WORLD);
+ }
+
+ int intSizeOf(std::vector<int> &data)
+ {
+ return data.size();
+ }
+
+ int intSizeOf(std::vector<double> &data)
+ {
+ return data.size();
+ }
+
+ int intSizeOf(std::vector<std::pair<int, int> > &data)
+ {
+ return data.size() * 2;
+ }
+
+ int intSizeOf(std::vector<AtomicBoundary> &data)
+ {
+ return data.size() * 3;
+ }
+
+ int intSizeOf(std::vector<BoundaryObject> &data)
+ {
+ return data.size() * 5;
+ }
+
+ int intSizeOf(std::vector<const DegreeOfFreedom*> &data)
+ {
+ return data.size();
+ }
+
+ void makeBuf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data, int* buf)
+ {
+ int i = 0;
+ for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator it = data.begin();
+ it != data.end(); ++it) {
+ buf[i++] = *(it->first);
+ buf[i++] = it->second;
+ }
+ }
+
+ void makeFromBuf(std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &data,
+ int *buf, int bufSize)
+ {
+ if (bufSize == 0)
+ return;
+
+ TEST_EXIT(bufSize % 2 == 0)("This should not happen!\n");
+
+ data.clear();
+ data.reserve(bufSize / 2);
+
+ for (int i = 0; i < (bufSize / 2); i++)
+ data.push_back(std::make_pair(buf[i * 2], buf[i * 2 + 1]));
+ }
+
+ void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf)
+ {
+ int pos = 0;
+ for (unsigned int i = 0; i < data.size(); i++) {
+ buf[pos++] = data[i].getCode().size();
+ buf[pos++] = data[i].getNumElements();
+
+ for (unsigned int j = 0; j < data[i].getCode().size(); j++)
+ buf[pos++] = data[i].getCode()[j];
+ }
+ }
+
+ void makeFromBuf(std::vector<MeshStructure> &data,
+ unsigned long int *buf, int bufSize)
+ {
+ int pos = 0;
+
+ while (pos < bufSize) {
+ int codeSize = buf[pos++];
+ int nElements = buf[pos++];
+ std::vector<unsigned long int> code;
+ code.resize(codeSize);
+ for (int i = 0; i < codeSize; i++)
+ code[i] = buf[pos++];
+
+ MeshStructure meshCode;
+ meshCode.init(code, nElements);
+
+ data.push_back(meshCode);
+ }
+ }
+
+ void makeBuf(std::vector<WorldVector<double> > &data, double *buf)
+ {
+ int dimOfWorld = Global::getGeo(WORLD);
+ int pos = 0;
+ for (unsigned int i = 0; i < data.size(); i++)
+ for (int j = 0; j < dimOfWorld; j++)
+ buf[pos++] = data[i][j];
+ }
+
+ void makeFromBuf(std::vector<WorldVector<double> > &data, double *buf, int bufSize)
+ {
+ int dimOfWorld = Global::getGeo(WORLD);
+ TEST_EXIT(bufSize % Global::getGeo(WORLD) == 0)("This should not happen!\n");
+
+ int pos = 0;
+ data.resize(bufSize / Global::getGeo(WORLD));
+ for (unsigned int i = 0; i < data.size(); i++)
+ for (int j = 0; j < dimOfWorld; j++)
+ data[i][j] = buf[pos++];
+ }
+
+ void makeBuf(std::vector<int> &data, int *buf)
+ {
+ for (unsigned int i = 0; i < data.size(); i++)
+ buf[i] = data[i];
+ }
+
+ void makeFromBuf(std::vector<int> &data, int *buf, int bufSize)
+ {
+ data.resize(bufSize);
+
+ for (int i = 0; i < bufSize; i++)
+ data[i] = buf[i];
+ }
+
+ void makeBuf(std::vector<double> &data, double *buf)
+ {
+ for (unsigned int i = 0; i < data.size(); i++)
+ buf[i] = data[i];
+ }
+
+ void makeFromBuf(std::vector<double> &data, double *buf, int bufSize)
+ {
+ data.resize(bufSize);
+
+ for (int i = 0; i < bufSize; i++)
+ data[i] = buf[i];
+ }
+
+ void makeBuf(std::vector<std::pair<int, int> > &data, int *buf)
+ {
+ for (unsigned int i = 0; i < data.size(); i++) {
+ buf[i * 2] = data[i].first;
+ buf[i * 2 + 1] = data[i].second;
+ }
+ }
+
+ void makeBuf(std::vector<AtomicBoundary> &data, int *buf)
+ {
+ for (unsigned int i = 0; i < data.size(); i++) {
+ buf[i * 3] = data[i].rankObj.elIndex;
+ buf[i * 3 + 1] = data[i].rankObj.subObj;
+ buf[i * 3 + 2] = data[i].rankObj.ithObj;
+ }
+ }
+
+ void makeFromBuf(std::vector<AtomicBoundary> &data, int *buf, int bufSize)
+ {
+ if (bufSize == 0)
+ return;
+
+ TEST_EXIT(bufSize % 3 == 0)("This should not happen!\n");
+
+ data.resize(bufSize / 3);
+ for (int i = 0; i < bufSize / 3; i++) {
+ data[i].rankObj.elIndex = buf[i * 3];
+ data[i].rankObj.subObj = static_cast<GeoIndex>(buf[i * 3 + 1]);
+ data[i].rankObj.ithObj = buf[i * 3 + 2];
+ }
+ }
+
+ void makeBuf(std::vector<BoundaryObject> &data, int *buf)
+ {
+ for (unsigned int i = 0; i < data.size(); i++) {
+ buf[i * 5] = data[i].elIndex;
+ buf[i * 5 + 1] = data[i].elType;
+ buf[i * 5 + 2] = data[i].subObj;
+ buf[i * 5 + 3] = data[i].ithObj;
+ buf[i * 5 + 4] = data[i].reverseMode;
+ }
+ }
+
+ void makeFromBuf(std::vector<BoundaryObject> &data, int *buf, int bufSize)
+ {
+ if (bufSize == 0)
+ return;
+
+ TEST_EXIT(bufSize % 5 == 0)("This should not happen!\n");
+
+ data.resize(bufSize / 5);
+ for (int i = 0; i < bufSize / 5; i++) {
+ data[i].elIndex = buf[i * 5];
+ data[i].elType = buf[i * 5 + 1];
+ data[i].subObj = static_cast<GeoIndex>(buf[i * 5 + 2]);
+ data[i].ithObj = buf[i * 5 + 3];
+ data[i].reverseMode = static_cast<bool>(buf[i * 5 + 4]);
+ }
+ }
+
+}
diff --git a/AMDiS/src/parallel/StdMpi.h b/AMDiS/src/parallel/StdMpi.h
index d7eea5fc0a302fca8b980ee34ab823d65272eca8..0c74719f4ae126d45e529fa257f28201f3df0471 100644
--- a/AMDiS/src/parallel/StdMpi.h
+++ b/AMDiS/src/parallel/StdMpi.h
@@ -29,220 +29,57 @@
namespace AMDiS {
- int intSizeOf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data)
- {
- return data.size() * 2;
- }
+ int intSizeOf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data);
- int intSizeOf(std::vector<MeshStructure> &data)
- {
- int s = 0;
- for (unsigned int i = 0; i < data.size(); i++)
- s += data[i].getCode().size() + 2;
+ int intSizeOf(std::vector<MeshStructure> &data);
- return s;
- }
+ int intSizeOf(std::vector<WorldVector<double> > &data);
- int intSizeOf(std::vector<WorldVector<double> > &data)
- {
- return data.size() * Global::getGeo(WORLD);
- }
-
- int intSizeOf(std::vector<int> &data)
- {
- return data.size();
- }
+ int intSizeOf(std::vector<int> &data);
- int intSizeOf(std::vector<double> &data)
- {
- return data.size();
- }
+ int intSizeOf(std::vector<double> &data);
- int intSizeOf(std::vector<std::pair<int, int> > &data)
- {
- return data.size() * 2;
- }
+ int intSizeOf(std::vector<std::pair<int, int> > &data);
- int intSizeOf(std::vector<AtomicBoundary> &data)
- {
- return data.size() * 3;
- }
+ int intSizeOf(std::vector<AtomicBoundary> &data);
- int intSizeOf(std::vector<BoundaryObject> &data)
- {
- return data.size() * 5;
- }
+ int intSizeOf(std::vector<BoundaryObject> &data);
- int intSizeOf(std::vector<const DegreeOfFreedom*> &data)
- {
- return data.size();
- }
+ int intSizeOf(std::vector<const DegreeOfFreedom*> &data);
- void makeBuf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data, int* buf)
- {
- int i = 0;
- for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator it = data.begin();
- it != data.end(); ++it) {
- buf[i++] = *(it->first);
- buf[i++] = it->second;
- }
- }
+ void makeBuf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data, int* buf);
void makeFromBuf(std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &data,
- int *buf, int bufSize)
- {
- if (bufSize == 0)
- return;
-
- TEST_EXIT(bufSize % 2 == 0)("This should not happen!\n");
-
- data.clear();
- data.reserve(bufSize / 2);
-
- for (int i = 0; i < (bufSize / 2); i++)
- data.push_back(std::make_pair(buf[i * 2], buf[i * 2 + 1]));
- }
+ int *buf, int bufSize);
- void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf)
- {
- int pos = 0;
- for (unsigned int i = 0; i < data.size(); i++) {
- buf[pos++] = data[i].getCode().size();
- buf[pos++] = data[i].getNumElements();
-
- for (unsigned int j = 0; j < data[i].getCode().size(); j++)
- buf[pos++] = data[i].getCode()[j];
- }
- }
+ void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf);
void makeFromBuf(std::vector<MeshStructure> &data,
- unsigned long int *buf, int bufSize)
- {
- int pos = 0;
+ unsigned long int *buf, int bufSize);
- while (pos < bufSize) {
- int codeSize = buf[pos++];
- int nElements = buf[pos++];
- std::vector<unsigned long int> code;
- code.resize(codeSize);
- for (int i = 0; i < codeSize; i++)
- code[i] = buf[pos++];
+ void makeBuf(std::vector<WorldVector<double> > &data, double *buf);
- MeshStructure meshCode;
- meshCode.init(code, nElements);
+ void makeFromBuf(std::vector<WorldVector<double> > &data, double *buf, int bufSize);
- data.push_back(meshCode);
- }
- }
+ void makeBuf(std::vector<int> &data, int *buf);
- void makeBuf(std::vector<WorldVector<double> > &data, double *buf)
- {
- int dimOfWorld = Global::getGeo(WORLD);
- int pos = 0;
- for (unsigned int i = 0; i < data.size(); i++)
- for (int j = 0; j < dimOfWorld; j++)
- buf[pos++] = data[i][j];
- }
-
- void makeFromBuf(std::vector<WorldVector<double> > &data, double *buf, int bufSize)
- {
- int dimOfWorld = Global::getGeo(WORLD);
- TEST_EXIT(bufSize % Global::getGeo(WORLD) == 0)("This should not happen!\n");
+ void makeFromBuf(std::vector<int> &data, int *buf, int bufSize);
- int pos = 0;
- data.resize(bufSize / Global::getGeo(WORLD));
- for (unsigned int i = 0; i < data.size(); i++)
- for (int j = 0; j < dimOfWorld; j++)
- data[i][j] = buf[pos++];
- }
+ void makeBuf(std::vector<double> &data, double *buf);
- void makeBuf(std::vector<int> &data, int *buf)
- {
- for (unsigned int i = 0; i < data.size(); i++)
- buf[i] = data[i];
- }
+ void makeFromBuf(std::vector<double> &data, double *buf, int bufSize);
- void makeFromBuf(std::vector<int> &data, int *buf, int bufSize)
- {
- data.resize(bufSize);
+ void makeBuf(std::vector<std::pair<int, int> > &data, int *buf);
- for (int i = 0; i < bufSize; i++)
- data[i] = buf[i];
- }
+ void makeBuf(std::vector<AtomicBoundary> &data, int *buf);
- void makeBuf(std::vector<double> &data, double *buf)
- {
- for (unsigned int i = 0; i < data.size(); i++)
- buf[i] = data[i];
- }
-
- void makeFromBuf(std::vector<double> &data, double *buf, int bufSize)
- {
- data.resize(bufSize);
+ void makeFromBuf(std::vector<AtomicBoundary> &data, int *buf, int bufSize);
- for (int i = 0; i < bufSize; i++)
- data[i] = buf[i];
- }
+ void makeBuf(std::vector<BoundaryObject> &data, int *buf);
- void makeBuf(std::vector<std::pair<int, int> > &data, int *buf)
- {
- for (unsigned int i = 0; i < data.size(); i++) {
- buf[i * 2] = data[i].first;
- buf[i * 2 + 1] = data[i].second;
- }
- }
+ void makeFromBuf(std::vector<BoundaryObject> &data, int *buf, int bufSize);
- void makeBuf(std::vector<AtomicBoundary> &data, int *buf)
- {
- for (unsigned int i = 0; i < data.size(); i++) {
- buf[i * 3] = data[i].rankObj.elIndex;
- buf[i * 3 + 1] = data[i].rankObj.subObj;
- buf[i * 3 + 2] = data[i].rankObj.ithObj;
- }
- }
- void makeFromBuf(std::vector<AtomicBoundary> &data, int *buf, int bufSize)
- {
- if (bufSize == 0)
- return;
-
- TEST_EXIT(bufSize % 3 == 0)("This should not happen!\n");
-
- data.resize(bufSize / 3);
- for (int i = 0; i < bufSize / 3; i++) {
- data[i].rankObj.elIndex = buf[i * 3];
- data[i].rankObj.subObj = static_cast<GeoIndex>(buf[i * 3 + 1]);
- data[i].rankObj.ithObj = buf[i * 3 + 2];
- }
- }
-
- void makeBuf(std::vector<BoundaryObject> &data, int *buf)
- {
- for (unsigned int i = 0; i < data.size(); i++) {
- buf[i * 5] = data[i].elIndex;
- buf[i * 5 + 1] = data[i].elType;
- buf[i * 5 + 2] = data[i].subObj;
- buf[i * 5 + 3] = data[i].ithObj;
- buf[i * 5 + 4] = data[i].reverseMode;
- }
- }
-
- void makeFromBuf(std::vector<BoundaryObject> &data, int *buf, int bufSize)
- {
- if (bufSize == 0)
- return;
-
- TEST_EXIT(bufSize % 5 == 0)("This should not happen!\n");
-
- data.resize(bufSize / 5);
- for (int i = 0; i < bufSize / 5; i++) {
- data[i].elIndex = buf[i * 5];
- data[i].elType = buf[i * 5 + 1];
- data[i].subObj = static_cast<GeoIndex>(buf[i * 5 + 2]);
- data[i].ithObj = buf[i * 5 + 3];
- data[i].reverseMode = static_cast<bool>(buf[i * 5 + 4]);
- }
- }
template<typename SendT, typename RecvT=SendT>
class StdMpi