diff --git a/.dir-locals.el b/.dir-locals.el
new file mode 100644
index 00000000..ba8984af
--- /dev/null
+++ b/.dir-locals.el
@@ -0,0 +1,101 @@
+; Adapted from PETSc
+
+((nil . ((indent-tabs-mode         . nil)
+         (show-trailing-whitespace . t)))
+
+ (c++-mode . ((c-tab-always-indent        . t)
+              (fill-column                . 100)
+              (c-basic-offset		  . 2)
+              (c-comment-only-line-offset . 0)
+              (c-hanging-braces-alist     . ((substatement-open after)
+                                             (brace-list-open after)
+                                             (brace-entry-open)
+                                             (defun-open after)
+                                             (class-open after)
+                                             (inline-open after)
+                                             (block-open after)
+                                             (block-close . c-snug-do-while)
+                                             (statement-case-open after)
+                                             (substatement after)))
+              (c-hanging-colons-alist     . ((member-init-intro before)
+                                             (inher-intro)
+                                             (case-label after)
+                                             (label after)
+                                             (access-label after)))
+              (c-cleanup-list             . (scope-operator
+                                             brace-else-brace
+                                             brace-elseif-brace
+                                             brace-catch-brace
+                                             empty-defun-braces
+                                             list-close-comma
+                                             defun-close-semi))
+              (c-offsets-alist            . ((inexpr-class . +)
+                                             (inexpr-statement . +)
+                                             (lambda-intro-cont . +)
+                                             (inlambda . c-lineup-inexpr-block)
+                                             (template-args-cont c-lineup-template-args +)
+                                             (incomposition . +)
+                                             (inmodule . +)
+                                             (innamespace . +)
+                                             (inextern-lang . +)
+                                             (composition-close . 0)
+                                             (module-close . 0)
+                                             (namespace-close . 0)
+                                             (extern-lang-close . 0)
+                                             (composition-open . 0)
+                                             (module-open . 0)
+                                             (namespace-open . 0)
+                                             (extern-lang-open . 0)
+                                             (friend . 0)
+                                             (cpp-define-intro c-lineup-cpp-define +)
+                                             (cpp-macro-cont . +)
+                                             (cpp-macro . [0])
+                                             (inclass . ++)
+                                             (stream-op . c-lineup-streamop)
+                                             (arglist-cont-nonempty c-lineup-gcc-asm-reg c-lineup-arglist)
+                                             (arglist-cont c-lineup-gcc-asm-reg 0)
+                                             (arglist-intro . +)
+                                             (catch-clause . 0)
+                                             (else-clause . 0)
+                                             (do-while-closure . 0)
+                                             (label . 2)
+                                             (access-label . -)
+                                             (substatement-label . 2)
+                                             (substatement . +)
+                                             (statement-case-open . 0)
+                                             (statement-case-intro . +)
+                                             (statement-block-intro . +)
+                                             (statement-cont . +)
+                                             (statement . 0)
+                                             (brace-entry-open . 0)
+                                             (brace-list-entry . 0)
+                                             (brace-list-intro . +)
+                                             (brace-list-close . 0)
+                                             (brace-list-open . 0)
+                                             (block-close . 0)
+                                             (inher-cont . c-lineup-multi-inher)
+                                             (inher-intro . +)
+                                             (member-init-cont . c-lineup-multi-inher)
+                                             (member-init-intro . +)
+                                             (annotation-var-cont . +)
+                                             (annotation-top-cont . 0)
+                                             (topmost-intro-cont . c-lineup-topmost-intro-cont)
+                                             (topmost-intro . 0)
+                                             (knr-argdecl . 0)
+                                             (func-decl-cont . +)
+                                             (inline-close . 0)
+                                             (inline-open . 0)
+                                             (class-close . 0)
+                                             (class-open . 0)
+                                             (defun-block-intro . +)
+                                             (defun-close . 0)
+                                             (defun-open . 0)
+                                             (string . c-lineup-dont-change)
+                                             (arglist-close . c-lineup-arglist)
+                                             (substatement-open . 0)
+                                             (case-label . 0)
+                                             (block-open . 0)
+                                             (c . 1)
+                                             (comment-intro . 0)
+                                             (knr-argdecl-intro . -)))
+              )))
diff --git a/.gitignore b/.gitignore
index 9e4a6fe9..3ae541f8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,6 +2,8 @@
 /doc/build
 /build*
 /dist
+/coverage
+*.info
 
 # Created by http://www.gitignore.io
 
diff --git a/.travis.yml b/.travis.yml
index 6cd61e9d..e6b26514 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,24 +1,46 @@
 language:
   - cpp
 
-compiler:
-  - gcc
-  - clang
+cache: apt
 
+notifications:
+  email:
+    on_success: change  # default: change
+    on_failure: always  # default: always
+
+# why does Travis not allow for multi-dimensional build matrices?!?
+#  -> https://github.com/travis-ci/travis-ci/issues/1519
 env:
-  - CMAKE_BUILD_TYPE="Release"
-  - CMAKE_BUILD_TYPE="Debug"
+  matrix:
+    - CMAKE_BUILD_TYPE="Release" WITH_MPI=OFF CXX=clang++ BOOST_VERSION=1.53
+    - CMAKE_BUILD_TYPE="Release" WITH_MPI=OFF CXX=g++ BOOST_VERSION=1.53
+    - CMAKE_BUILD_TYPE="Release" WITH_MPI=ON CXX=g++ BOOST_VERSION=1.53
+    - CMAKE_BUILD_TYPE="Debug" WITH_MPI=OFF CXX=clang++ BOOST_VERSION=1.53
+    - CMAKE_BUILD_TYPE="Debug" WITH_MPI=OFF CXX=g++ BOOST_VERSION=1.53
+    - CMAKE_BUILD_TYPE="Debug" WITH_MPI=ON CXX=g++ BOOST_VERSION=1.53
+    - CMAKE_BUILD_TYPE="Release" WITH_MPI=OFF CXX=clang++ BOOST_VERSION=1.55
+    - CMAKE_BUILD_TYPE="Release" WITH_MPI=OFF CXX=g++ BOOST_VERSION=1.55
+    - CMAKE_BUILD_TYPE="Release" WITH_MPI=ON CXX=g++ BOOST_VERSION=1.55
+    - CMAKE_BUILD_TYPE="Debug" WITH_MPI=OFF CXX=clang++ BOOST_VERSION=1.55
+    - CMAKE_BUILD_TYPE="Debug" WITH_MPI=OFF CXX=g++ BOOST_VERSION=1.55
+    - CMAKE_BUILD_TYPE="Debug" WITH_MPI=ON CXX=g++ BOOST_VERSION=1.55
 
 before_install:
-  - sudo add-apt-repository ppa:apokluda/boost1.53 --yes
-  - sudo add-apt-repository ppa:ubuntu-toolchain-r/test --yes # libstdc++-4.8
-  - if [ "${CXX}" == "clang++" ]; then sudo add-apt-repository --yes ppa:h-rayflood/llvm; fi # clang++-3.2
+  # boost packages
+  - if [ "${BOOST_VERSION}" == "1.53" ]; then sudo add-apt-repository ppa:apokluda/boost1.53 --yes; fi
+  - if [ "${BOOST_VERSION}" == "1.55" ]; then sudo add-apt-repository ppa:boost-latest/ppa --yes; fi
+  # libstdc++-4.8
+  - sudo add-apt-repository ppa:ubuntu-toolchain-r/test --yes
+  # clang++-3.2
+  - if [ "${CXX}" == "clang++" ]; then sudo add-apt-repository --yes ppa:h-rayflood/llvm; fi
   - sudo apt-get update
 
 install:
+  # note: clang needs g++-4.8 it for libstdc++ update
   - if [ "${CXX}" == "clang++" ]; then sudo apt-get -qq install clang-3.2; fi
-  - sudo apt-get install g++-4.8 # clang need it for libstdc++ update
-  - sudo apt-get install libboost1.53-all-dev
+  - if [ "${WITH_MPI}" == "ON" ]; then sudo apt-get -qq install mpich2 libmpich2-dev; fi
+  - sudo apt-get install g++-4.8
+  - sudo apt-get install "libboost${BOOST_VERSION}-all-dev" libfftw3-3 libfftw3-dev
 
 before_script:
   # update compilers
@@ -28,7 +50,7 @@ before_script:
 script:
   - mkdir build
   - cd build
-  - cmake -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -Dpfasst_DISABLE_LIBCXX=ON -Dpfasst_BUILD_TESTS=ON -Dpfasst_BUILD_EXAMPLES=ON ..
+  - cmake -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -Dpfasst_WITH_MPI=${WITH_MPI} -Dpfasst_DISABLE_LIBCXX=ON -Dpfasst_BUILD_TESTS=ON -Dpfasst_BUILD_EXAMPLES=ON ..
   - make
-  - make test
-
+  - make test ARGS="-V"
+  - make run_example_all
diff --git a/3rdparty/CMakeLists.txt b/3rdparty/CMakeLists.txt
index 3edacc88..3592cb77 100644
--- a/3rdparty/CMakeLists.txt
+++ b/3rdparty/CMakeLists.txt
@@ -21,31 +21,40 @@ list(APPEND 3rdparty_INCLUDES ${Boost_INCLUDE_DIRS})
 
 
 if(pfasst_BUILD_EXAMPLES)
-    message(STATUS "FFTW3")
-    set(fftw3_SOURCE_DIR "${pfasst_BINARY_DIR}/3rdparty/src/fftw3")
-    set(fftw3_INSTALL_DIR "${pfasst_BINARY_DIR}/3rdparty/src/fftw3-install")
+    find_package(FFTW)
 
-    ExternalProject_Add(
-        fftw3
-        LIST_SEPARATOR " "
-        URL http://fftw.org/fftw-3.3.4.tar.gz
-        URL_MD5 2edab8c06b24feeb3b82bbb3ebf3e7b3
-        TIMEOUT 30
-        UPDATE_COMMAND ""
-        PATCH_COMMAND ""
-        BUILD_IN_SOURCE ON
-        CONFIGURE_COMMAND ${fftw3_SOURCE_DIR}/configure CC=${CMAKE_C_COMPILER} CXX=${CMAKE_CXX_COMPILER} CXXFLAGS=${CMAKE_CXX_FLAGS} --prefix=${fftw3_INSTALL_DIR} --libdir=${fftw3_INSTALL_DIR}/lib
-        BUILD_COMMAND make
-        TEST_COMMAND ""
-        INSTALL_DIR ${fftw3_SOURCE_DIR}-install
-        INSTALL_COMMAND make install
-        LOG_DOWNLOAD ON
-        LOG_CONFIGURE ON
-        LOG_BUILD ON
-    )
-    # Specify include dir
-    set(fftw3_INCLUDES ${fftw3_INSTALL_DIR}/include PARENT_SCOPE)
-    set(fftw3_LIBS ${fftw3_INSTALL_DIR}/lib/libfftw3.a PARENT_SCOPE)
+    if(NOT FFTW_FOUND)
+      message(STATUS "FFTW3 not found on your system.")
+      message(STATUS "Going to download and compile it automatically.")
+      set(fftw3_SOURCE_DIR "${pfasst_BINARY_DIR}/3rdparty/src/fftw3")
+      set(fftw3_INSTALL_DIR "${pfasst_BINARY_DIR}/3rdparty/src/fftw3-install")
+
+      ExternalProject_Add(
+          fftw3
+          LIST_SEPARATOR " "
+          URL http://fftw.org/fftw-3.3.4.tar.gz
+          URL_MD5 2edab8c06b24feeb3b82bbb3ebf3e7b3
+          TIMEOUT 30
+          UPDATE_COMMAND ""
+          PATCH_COMMAND ""
+          BUILD_IN_SOURCE ON
+          CONFIGURE_COMMAND ${fftw3_SOURCE_DIR}/configure CC=${CMAKE_C_COMPILER} CXX=${CMAKE_CXX_COMPILER} CXXFLAGS=${CMAKE_CXX_FLAGS} --prefix=${fftw3_INSTALL_DIR} --libdir=${fftw3_INSTALL_DIR}/lib
+          BUILD_COMMAND make
+          TEST_COMMAND ""
+          INSTALL_DIR ${fftw3_SOURCE_DIR}-install
+          INSTALL_COMMAND make install
+          LOG_DOWNLOAD ON
+          LOG_CONFIGURE ON
+          LOG_BUILD ON
+      )
+      # Specify include dir
+      set(FFTW_INCLUDE_PATH ${fftw3_INSTALL_DIR}/include PARENT_SCOPE)
+      set(FFTW_LIBRARIES ${fftw3_INSTALL_DIR}/lib/libfftw3.a PARENT_SCOPE)
+    else()
+      set(FFTW_INCLUDE_PATH ${FFTW_INCLUDE_PATH} PARENT_SCOPE)
+      set(FFTW_LIBRARIES ${FFTW_LIBRARIES} PARENT_SCOPE)
+    endif()
+    set(FFTW_FOUND ${FFTW_FOUND} PARENT_SCOPE)
 
     if(pfasst_BUILD_TESTS)
         message(STATUS "--------------------------------------------------------------------------------")
@@ -58,10 +67,9 @@ if(pfasst_BUILD_TESTS)
     # Add gmock
     ExternalProject_Add(
       googlemock
-      SVN_REPOSITORY http://googlemock.googlecode.com/svn/trunk/
-      SVN_REVISION -r 449  # release 1.7.0
-      TIMEOUT 10
-      # Disable SVN update
+      URL http://googlemock.googlecode.com/files/gmock-1.7.0.zip
+      URL_MD5 073b984d8798ea1594f5e44d85b20d66
+      TIMEOUT 30
       UPDATE_COMMAND ""
       PATCH_COMMAND ""
       CMAKE_ARGS -DCMAKE_BUILD_TYPE=Release
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 133f05db..7c0d9b8e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,62 @@
 # Changelog                                                                        {#page_changelog}
 
+## v0.2.0 -- MPI PFASST (2014/08/XX)
+
+DOI: [10.5281/zenodo.11047](http://dx.doi.org/10.5281/zenodo.11047)
+
+### Notable Features
+
+* Addition of MPI based PFASST.
+
+### Details
+
+* Addition of MPI based PFASST algorithm using the standard predictor stage in a block mode with 
+  fixed iterations.
+  ([#46][], [#57][], [#59][])
+
+* Addition of a simple scalar example and appropriate tests.
+  ([#61][], [#63][], [#76][], [#81][])
+
+* Further tests on proper calculation of quadrature nodes and weights
+  ([#74][])
+
+* Better handling of 3rd-party dependencies
+  ([#53][], [#55][])
+
+* Various tidying
+  ([#56][], [#60][], [#77][])
+
+* Basic Profiling support and test coverage report
+  ([#78][])
+
+[#46]: https://github.com/Parallel-in-Time/PFASST/pull/46
+[#57]: https://github.com/Parallel-in-Time/PFASST/pull/56
+[#59]: https://github.com/Parallel-in-Time/PFASST/pull/59
+[#53]: https://github.com/Parallel-in-Time/PFASST/pull/53
+[#55]: https://github.com/Parallel-in-Time/PFASST/pull/55
+[#56]: https://github.com/Parallel-in-Time/PFASST/pull/56
+[#60]: https://github.com/Parallel-in-Time/PFASST/pull/60
+[#61]: https://github.com/Parallel-in-Time/PFASST/pull/61
+[#63]: https://github.com/Parallel-in-Time/PFASST/pull/63
+[#74]: https://github.com/Parallel-in-Time/PFASST/pull/74
+[#76]: https://github.com/Parallel-in-Time/PFASST/pull/76
+[#77]: https://github.com/Parallel-in-Time/PFASST/pull/77
+[#78]: https://github.com/Parallel-in-Time/PFASST/pull/78
+[#78]: https://github.com/Parallel-in-Time/PFASST/pull/81
+
+### Contributors
+
+* Matthew Emmett, Lawrence Berkeley National Laboratory ([memmett][])
+* Torbjörn Klatt, Jülich Supercomputing Centre ([torbjoernk][])
+* Daniel Ruprecht, Institute of Computational Science, University of Lugano ([danielru][])
+* Robert Speck, Jülich Supercomputing Centre ([pancetta][])
+
+[memmett]: https://github.com/memmett
+[torbjoernk]: https://github.com/torbjoernk
+[danielru]: https://github.com/danielru
+[pancetta]: https://github.com/pancetta
+
+
 ## v0.1.0 -- First Release (2014/07/25)
 
 DOI: [10.5281/zenodo.11047](http://dx.doi.org/10.5281/zenodo.11047)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7e923660..8b8aed5a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,25 +1,31 @@
 cmake_minimum_required(VERSION 2.8)
-
-set(CMAKE_VERBOSE_MAKEFILE ON)
-
 project(pfasst)
 
+list(APPEND CMAKE_MODULE_PATH ${pfasst_SOURCE_DIR}/cmake)
 include(cmake/utility_functions.cmake)
 include(CheckCXXCompilerFlag)
 include(ExternalProject)
+
 # Set default ExternalProject root directory
 set_directory_properties(PROPERTIES EP_PREFIX ${CMAKE_BINARY_DIR}/3rdparty)
 
-option(pfasst_DISABLE_LIBCXX "Disable use of LLVM's libc++ when compiling with Clang." ON)
-option(pfasst_BUILD_EXAMPLES "Build example programs." ON)
-option(pfasst_BUILD_TESTS "Build test suite for PFASST." ON)
+option(pfasst_DISABLE_LIBCXX "Disable use of LLVM's libc++ when compiling with Clang." ON )
+option(pfasst_BUILD_EXAMPLES "Build example programs."                                 ON )
+option(pfasst_BUILD_TESTS    "Build test suite for PFASST."                            ON )
+option(pfasst_WITH_MPI       "Build with MPI enabled."                                 OFF)
+option(pfasst_WITH_GCC_PROF  "Enable excessive debugging & profiling output with GCC." OFF)
 
-# output directories
+# Set output directories
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${pfasst_SOURCE_DIR}/dist/bin")
 set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${pfasst_SOURCE_DIR}/dist/lib")
 set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY "${pfasst_SOURCE_DIR}/dist/lib")
 
-# check for C++11 support
+if(${pfasst_WITH_MPI})
+    find_package(MPI REQUIRED)
+    message(STATUS "Using MPI C++ Compiler: ${MPI_CXX_COMPILER}")
+endif()
+
+# Check for C++11 support
 if(${CMAKE_CXX_COMPILER_ID} MATCHES GNU)
     check_cxx_compiler_flag(-std=c++11 HAVE_STD11)
     if(HAVE_STD11)
@@ -27,6 +33,9 @@ if(${CMAKE_CXX_COMPILER_ID} MATCHES GNU)
     else()
         message(FATAL_ERROR "No advanced standard C++ support of your GCC (-std=c++11 not defined).")
     endif()
+    if(pfasst_WITH_GCC_PROF)
+        add_to_string_list("${CMAKE_CXX_FLAGS}" CMAKE_CXX_FLAGS "-ggdb3 -pg")
+    endif(pfasst_WITH_GCC_PROF)
 elseif(${CMAKE_CXX_COMPILER_ID} MATCHES Clang)
     check_cxx_compiler_flag(-std=c++11 HAVE_STD11)
     if(HAVE_STD11)
@@ -57,7 +66,7 @@ else()
 endif()
 message(STATUS "Your compiler has C++11 support. Hurray!")
 
-# enable all compiler warnings
+# Enable all compiler warnings
 add_to_string_list("${CMAKE_CXX_FLAGS}" CMAKE_CXX_FLAGS "-Wall -Wextra -Wpedantic")
 
 set(3rdparty_INCLUDES)
@@ -69,13 +78,14 @@ if(pfasst_BUILD_TESTS)
     enable_testing()
 endif(pfasst_BUILD_TESTS)
 
-# adding / including 3rd-party libraries
+# Add / include 3rd-party libraries
 message(STATUS "********************************************************************************")
 message(STATUS "Configuring 3rd party libraries")
 # makes available:
-#  - fftw3_INCLUDES  (if pfasst_BUILD_EXAMPLES)
-#  - fftw3_LIBS      (if pfasst_BUILD_EXAMPLES)
-# and Boost headers in 3rdparty_INCLUDES
+#  - Boost headers in 3rdparty_INCLUDES
+#  - Google test and mock headers in 3rdparty_INCLUDES (if pfasst_BUILD_TESTS)
+#  - FFTW_INCLUDE_PATH (if pfasst_BUILD_EXAMPLES)
+#  - FFTW_LIBRARIES (if pfasst_BUILD_EXAMPLES)
 add_subdirectory(3rdparty)
 
 message(STATUS "********************************************************************************")
@@ -95,6 +105,9 @@ if(pfasst_BUILD_TESTS)
 endif()
 message(STATUS "********************************************************************************")
 
+message(STATUS "C++ Compiler: ${CMAKE_CXX_COMPILER}")
 message(STATUS "C++ Compiler ID: ${CMAKE_CXX_COMPILER_ID}")
+message(STATUS "C++ Compiler Names: ${CMAKE_CXX_COMPILER_NAMES}")
+message(STATUS "C++ Compiler Version: ${CMAKE_CXX_COMPILER_VERSION}")
 message(STATUS "C++ Flags: ${CMAKE_CXX_FLAGS}")
 message(STATUS "C++ link flags: ${CMAKE_CXX_LINK_FLAGS}")
diff --git a/README.md b/README.md
index d5a8dcf3..8b7bb05a 100644
--- a/README.md
+++ b/README.md
@@ -1,18 +1,19 @@
 PFASST                                                                                   {#mainpage}
 ======
 
-The PFASST algorithm is a time-parallel algorithm for solving ODEs and PDEs.
-
-The PFASST project is a C++ implementation of the *parallel full approximation scheme in space and 
-time* (PFASST) algorithm.
-It also contains basic implementations of the *spectral deferred correction* (SDC) and 
-*multi-level spectral deferred correction* (MLSDC) algorithms.
+The PFASST project is a C++ implementation of the *parallel full approximation scheme in space and
+time* (PFASST) algorithm, which in turn is a time-parallel algorithm for solving ODEs and PDEs.  It
+also contains basic implementations of the *spectral deferred correction* (SDC) and *multi-level
+spectral deferred correction* (MLSDC) algorithms.
 
 
 News
 ----
 
-* July 25, 2014: PFASST v0.1.0 released.  Please see the [release notes](#releases) for more 
+* August 29, 2014: PFASST v0.2.0 release. Please sse the [release notes](#releases) for more 
+  information.
+
+* July 25, 2014: PFASST v0.1.0 released.  Please see the [release notes](#releases) for more
   information.
 
 
@@ -32,15 +33,22 @@ Currently, it features the following content:
 * \subpage #page_examples
 * \subpage #page_contributing
 * \subpage #page_style_guide
+* \subpage #page_troubleshooting
 
 
 Releases
 --------
 
+* **v0.2.0** MPI-PFASST Release (2014/08/29)
+
+  Implementation of PFASST with MPI.
+  DOI: [10.5281/zenodo.11047][DOI_v010]
+  See \subpage #page_changelog "the Changelog" for details.
+
 * **v0.1.0** First Release (2014/07/25)
 
-  Initial release with basic implementations of SDC and MLSDC.  
-  DOI: [10.5281/zenodo.11047][DOI_v010]  
+  Initial release with basic implementations of SDC and MLSDC.
+  DOI: [10.5281/zenodo.11047][DOI_v010]
   See \subpage #page_changelog "the Changelog" for details.
 
 [DOI_v010]: http://dx.doi.org/10.5281/zenodo.11047
@@ -60,8 +68,11 @@ Build status
 
 For details see [Travis][pfasst-travis].
 
+A test coverage report is generated and published by [Jenkins][coverage-report].
+
 
-[documentation]: https://pint.fz-juelich.de/ci/view/PFASST/job/PFASST%20%28Docu%29/doxygen/
-[pfasst-travis]: https://travis-ci.org/Parallel-in-Time/PFASST
-[master-status]: https://travis-ci.org/Parallel-in-Time/PFASST.svg?branch=master
-[dev-status]:    https://travis-ci.org/Parallel-in-Time/PFASST.svg?branch=development
+[documentation]:   https://pint.fz-juelich.de/ci/view/PFASST/job/PFASST%20%28Docu%29/doxygen/
+[pfasst-travis]:   https://travis-ci.org/Parallel-in-Time/PFASST
+[master-status]:   https://travis-ci.org/Parallel-in-Time/PFASST.svg?branch=master
+[dev-status]:      https://travis-ci.org/Parallel-in-Time/PFASST.svg?branch=development
+[coverage-report]: https://pint.fz-juelich.de/ci/job/PFASST_GCC_4-8_Coverage/Test_Coverage
diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake
new file mode 100644
index 00000000..03969b53
--- /dev/null
+++ b/cmake/FindFFTW.cmake
@@ -0,0 +1,22 @@
+# - Find FFTW
+# Find the native FFTW3 includes and library
+#
+# FFTW_INCLUDE_PATH - where to find fftw3.h
+# FFTW_LIBRARIES - List of libraries when using FFTW.
+# FFTW_FOUND - True if FFTW found.
+
+if(FFTW_INCLUDE_PATH)
+    # Already in cache, be silent
+    set(FFTW_FIND_QUIETLY TRUE)
+endif(FFTW_INCLUDE_PATH)
+
+find_path(FFTW_INCLUDE_PATH fftw3.h)
+
+find_library(FFTW_LIBRARIES NAMES fftw3)
+
+# handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if
+# all listed variables are TRUE
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(FFTW DEFAULT_MSG FFTW_LIBRARIES FFTW_INCLUDE_PATH)
+
+mark_as_advanced(FFTW_LIBRARIES FFTW_INCLUDE_PATH)
diff --git a/doc/source/examples.md b/doc/source/examples.md
index b7e70c78..f68f5970 100644
--- a/doc/source/examples.md
+++ b/doc/source/examples.md
@@ -3,3 +3,8 @@
 ## Simple Advection Diffusion
 
 See \subpage page_examples_advection_diffusion.
+
+
+## Scalar
+
+See \subpage page_examples_scalar.
diff --git a/doc/source/installing.md b/doc/source/installing.md
index b3c23337..e02d5d52 100644
--- a/doc/source/installing.md
+++ b/doc/source/installing.md
@@ -20,27 +20,27 @@ Just use _Git_ to clone the repository.
 
    * By default everything gets built, i.e. tests and examples.
 
-   * To deactivate the tests, pass `-Dpfasst_BUILD_TESTS=OFF` to the
-     _CMake_ command line.
+   * To deactivate the tests, pass `-Dpfasst_BUILD_TESTS=OFF` to the _CMake_ command line.
 
-   * To deactivate the examples, pass `-Dpfasst_BUILD_EXAMPLES=OFF` to
-     the _CMake_ command line.
+   * To deactivate the examples, pass `-Dpfasst_BUILD_EXAMPLES=OFF` to the _CMake_ command line.
 
-   * If you are not on a _Darvin_ system (i.e. MacOS, OSX, ...) and
-     you are using _Clang_ as the compiler, you want to deactivate the
-     use of LLVM's libc++ by passing `-Dpfasst_DISABLE_LIBCXX=ON` to
-     the _CMake_ command line.
+   * To enable profiling support, you need to specify `-Dpfasst_WITH_GCC_PROF` and have the GNU 
+     compiler selected.
 
-   * For example, to build with _Clang_ on a Linux system without the
-     examples, the call to _CMake_ looks like:
+   * If you are not on a _Darvin_ system (i.e. MacOS, OSX, ...) and you are using _Clang_ as the 
+     compiler, you want to deactivate the use of LLVM's libc++ by passing 
+     `-Dpfasst_DISABLE_LIBCXX=ON` to the _CMake_ command line. (this is the default)
 
-       cmake -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -Dpfasst_DISABLE_LIBCXX=ON -Dpfasst_BUILD_EXAMPLES=OFF ..
+   * For example, to build with _Clang_ on a Linux system without the examples, the call to 
+     _CMake_ looks like:
 
-   The [Google Testing] and [Mocking Framework] are required for the
-   unit tests, and are automatically downloaded and built.
+       cmake -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ [-Dpfasst_DISABLE_LIBCXX=ON] -Dpfasst_BUILD_EXAMPLES=OFF ..
 
-   The [FFTW3] library is required for (some) of the examples, and it
-   automatically downloaded and built as well.
+   The [Google Testing] and [Mocking Framework] are required for the unit tests, and are 
+   automatically downloaded and built.
+
+   The [FFTW3] library is required for (some) of the examples. It will be automatically downloaded 
+   and built as well if it is not found on the system.
 
 3. Run the unit tests (if unit tests have been built):
 
@@ -48,7 +48,9 @@ Just use _Git_ to clone the repository.
 
 4. Run the examples (if they have been built):
 
-       ./examples/advection
+       make run_example_all
+
+   Or go to `./dist/examples` and look for the exetuable you want to run.
 
 [Google Testing]: https://code.google.com/p/googletest/
 [Mocking Framework]: https://code.google.com/p/googlemock/
diff --git a/doc/source/troubleshooting.md b/doc/source/troubleshooting.md
new file mode 100644
index 00000000..ebeb4f64
--- /dev/null
+++ b/doc/source/troubleshooting.md
@@ -0,0 +1,35 @@
+# Trouble Shooting                                                           {#page_troubleshooting}
+
+## I'm getting SEGFAULTs, how do I track them down?
+
+Please try running your program through GDB and/or valgrind.
+
+## How do I run my MPI program through GDB?
+
+To run your MPI program through GDB:
+
+    mpirun -n 4 xterm -e gdb -x my-gdb-commands my-pfasst-program
+
+where `my-gdb-commands` is a plain text file containing GDB commands
+to run your program, and `my-pfasst-program` is your executable.  A
+very simple set of GDB commands for the `my-gdb-commands` is:
+
+    catch throw
+    run
+
+The result is: `mpirun` will launch 4 `xterm` windows, each of which
+will immediately run your program through `gdb`, which subsequently
+runs your program.  Hopefully GDB will catch the SEGFAULT.
+
+## How do I run my MPI program through valgrind?
+
+To run your MPI program through valgrind:
+
+    mpirun -n 4 xterm -e valgrind my-pfasst-program
+
+
+
+
+
+
+
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index b79278b1..5ae8e0d3 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1 +1,10 @@
+set(all_example_runs)
+
 add_subdirectory(advection_diffusion)
+add_subdirectory(scalar)
+
+add_custom_target(run_example_all
+    DEPENDS ${all_example_runs}
+    COMMENT "Running all examples" VERBATIM
+)
+message(STATUS "  - run_example_all")
diff --git a/examples/advection_diffusion/CMakeLists.txt b/examples/advection_diffusion/CMakeLists.txt
index ffa92cbc..f4d2c046 100644
--- a/examples/advection_diffusion/CMakeLists.txt
+++ b/examples/advection_diffusion/CMakeLists.txt
@@ -1,7 +1,10 @@
+set(all_example_runs ${all_example_runs})
+set(all_advec_example_runs)
+
 message(STATUS "  advection_diffusion")
 include_directories(
     ${3rdparty_INCLUDES}
-    ${fftw3_INCLUDES}
+    ${FFTW_INCLUDE_PATH}
     ${pfasst_INCLUDES}
 )
 
@@ -11,13 +14,72 @@ set(advec_examples
     serial_mlsdc_autobuild
 )
 
-foreach(example ${advec_examples})
+set(advec_mpi_examples)
+
+if(${pfasst_WITH_MPI})
+    set(advec_mpi_examples
+        mpi_pfasst
+    )
+endif()
+
+set(all_advec_examples
+    ${advec_examples}
+    ${advec_mpi_examples}
+)
+
+foreach(example ${all_advec_examples})
     add_executable(${example} ${CMAKE_CURRENT_SOURCE_DIR}/${example}.cpp)
-    add_dependencies(${example} fftw3)
+    if(NOT FFTW_FOUND)
+        add_dependencies(${example} fftw3)
+    endif()
     target_link_libraries(${example}
-        ${fftw3_LIBS}
+        ${FFTW_LIBRARIES}
     )
     set_target_properties(${example}
         PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/advection_diffusion
     )
 endforeach(example)
+
+message(STATUS "    creating make targets to run examples")
+foreach(example ${advec_examples})
+    add_custom_target(run_example_advec_${example}
+        COMMAND ${example}
+        DEPENDS ${example}
+        WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/advection_diffusion
+        COMMENT "Running Advection-Diffusion example: ${example}" VERBATIM
+    )
+    list(APPEND all_advec_example_runs run_example_advec_${example})
+    message(STATUS "      - run_example_advec_${example}")
+endforeach(example)
+
+if(${pfasst_WITH_MPI})
+    include_directories(${MPI_CXX_INCLUDE_PATH})
+    foreach(example ${advec_mpi_examples})
+        if(MPI_COMPILE_FLAGS)
+            set_target_properties(${example} PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
+        endif()
+        if(MPI_LINK_FLAGS)
+            set_target_properties(${example} PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+        endif()
+        target_link_libraries(${example} ${MPI_CXX_LIBRARIES})
+    endforeach(example)
+    message(STATUS "    creating make targets to run MPI examples")
+    foreach(example ${advec_mpi_examples})
+        add_custom_target(run_example_advec_${example}
+            COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/advection_diffusion/${example} ${MPIEXEC_POSTFLAGS}
+            DEPENDS ${example}
+            WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/advection_diffusion
+            COMMENT "Running Advection-Diffusion example: ${example}" VERBATIM
+        )
+        list(APPEND all_advec_example_runs run_example_advec_${example})
+        message(STATUS "      - run_example_advec_${example}")
+    endforeach(example)
+endif()
+
+add_custom_target(run_example_advec_all
+    DEPENDS ${all_advec_example_runs}
+    COMMENT "Running all Advection-Diffusion examples" VERBATIM
+)
+message(STATUS "    - run_example_advec_all")
+
+set(all_example_runs ${all_example_runs} ${all_advec_example_runs} PARENT_SCOPE)
diff --git a/examples/advection_diffusion/advection_diffusion_sweeper.hpp b/examples/advection_diffusion/advection_diffusion_sweeper.hpp
index dbdc39a1..f4dc3d97 100644
--- a/examples/advection_diffusion/advection_diffusion_sweeper.hpp
+++ b/examples/advection_diffusion/advection_diffusion_sweeper.hpp
@@ -5,12 +5,15 @@
 #ifndef _ADVECTION_DIFFUSION_SWEEPER_HPP_
 #define _ADVECTION_DIFFUSION_SWEEPER_HPP_
 
+#include <cstdlib>
 #include <cassert>
 #include <complex>
 #include <map>
 #include <ostream>
 #include <vector>
+#include <utility>
 
+#include <pfasst/globals.hpp>
 #include <pfasst/encap/imex_sweeper.hpp>
 
 #include "fft.hpp"
@@ -24,42 +27,150 @@ using pfasst::encap::Encapsulation;
 using pfasst::encap::as_vector;
 
 
+/**
+ * errors at different iterations and time nodes
+ * 
+ * Mapping a pair of step/iteration indices onto the error of the solution.
+ */
 typedef map<pair<size_t, size_t>, double> error_map;
 
 template<typename time = pfasst::time_precision>
 class AdvectionDiffusionSweeper
   : public pfasst::encap::IMEXSweeper<time>
 {
+    //! @{
     FFT fft;
-
     vector<complex<double>> ddx, lap;
+    //! @}
+
+    //! @{
     error_map errors;
+    //! @}
 
+    //! @{
     double v  = 1.0;
     time   t0 = 1.0;
     double nu = 0.02;
     size_t nf1evals = 0;
+    //! @}
 
   public:
+    //! @{
     AdvectionDiffusionSweeper(size_t nvars)
     {
-      ddx.resize(nvars);
-      lap.resize(nvars);
+      this->ddx.resize(nvars);
+      this->lap.resize(nvars);
       for (size_t i = 0; i < nvars; i++) {
         double kx = 2 * PI * ((i <= nvars / 2) ? int(i) : int(i) - int(nvars));
-        ddx[i] = complex<double>(0.0, 1.0) * kx;
-        lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
+        this->ddx[i] = complex<double>(0.0, 1.0) * kx;
+        this->lap[i] = (kx * kx < 1e-13) ? 0.0 : -kx * kx;
       }
     }
 
-    ~AdvectionDiffusionSweeper()
+    virtual ~AdvectionDiffusionSweeper()
+    {
+      cout << "number of f1 evals: " << this->nf1evals << endl;
+    }
+    //! @}
+
+    //! @{
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::predict()
+     */
+    void predict(bool initial) override
+    {
+      pfasst::encap::IMEXSweeper<time>::predict(initial);
+      time t  = this->get_controller()->get_time();
+      time dt = this->get_controller()->get_time_step();
+      this->echo_error(t + dt, true);
+    }
+
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::sweep()
+     */
+    void sweep() override
     {
-      cout << "number of f1 evals: " << nf1evals << endl;
+      pfasst::encap::IMEXSweeper<time>::sweep();
+      time t  = this->get_controller()->get_time();
+      time dt = this->get_controller()->get_time_step();
+      this->echo_error(t + dt);
+    }
+    //! @}
+
+    //! @{
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::f_expl_eval()
+     */
+    void f_expl_eval(shared_ptr<Encapsulation<time>> f_expl_encap,
+                     shared_ptr<Encapsulation<time>> u_encap,
+                     time t) override
+    {
+      UNUSED(t);
+      auto& u = as_vector<double, time>(u_encap);
+      auto& f_expl = as_vector<double, time>(f_expl_encap);
+
+      double c = -v / double(u.size());
+
+      auto* z = this->fft.forward(u);
+      for (size_t i = 0; i < u.size(); i++) {
+        z[i] *= c * this->ddx[i];
+      }
+      this->fft.backward(f_expl);
+
+      this->nf1evals++;
+    }
+
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::f_impl_eval()
+     */
+    void f_impl_eval(shared_ptr<Encapsulation<time>> f_impl_encap,
+                     shared_ptr<Encapsulation<time>> u_encap,
+                     time t) override
+    {
+      UNUSED(t);
+      auto& u = as_vector<double, time>(u_encap);
+      auto& f_impl = as_vector<double, time>(f_impl_encap);
+
+      double c = nu / double(u.size());
+
+      auto* z = this->fft.forward(u);
+      for (size_t i = 0; i < u.size(); i++) {
+        z[i] *= c * this->lap[i];
+      }
+      this->fft.backward(f_impl);
+    }
+
+    /**
+     * @copybrief pfasst::encap::IMEXSweeper::impl_solve()
+     */
+    void impl_solve(shared_ptr<Encapsulation<time>> f_impl_encap,
+                    shared_ptr<Encapsulation<time>> u_encap,
+                    time t, time dt,
+                    shared_ptr<Encapsulation<time>> rhs_encap) override
+    {
+      UNUSED(t);
+      auto& u = as_vector<double, time>(u_encap);
+      auto& f_impl = as_vector<double, time>(f_impl_encap);
+      auto& rhs = as_vector<double, time>(rhs_encap);
+
+      double c = nu * double(dt);
+
+      auto* z = this->fft.forward(rhs);
+      for (size_t i = 0; i < u.size(); i++) {
+        z[i] /= (1.0 - c * this->lap[i]) * double(u.size());
+      }
+      this->fft.backward(u);
+
+      for (size_t i = 0; i < u.size(); i++) {
+        f_impl[i] = (u[i] - rhs[i]) / double(dt);
+      }
     }
+    //! @}
 
+    //! @{
     void exact(shared_ptr<Encapsulation<time>> q, time t)
     {
-      this->exact(as_vector<double,time>(q), t);
+      this->exact(as_vector<double, time>(q), t);
     }
 
     void exact(DVectorT& q, time t)
@@ -81,10 +192,10 @@ class AdvectionDiffusionSweeper
 
     void echo_error(time t, bool predict = false)
     {
-      auto& qend = as_vector<double,time>(this->get_end_state());
+      auto& qend = as_vector<double, time>(this->get_end_state());
       DVectorT qex(qend.size());
 
-      exact(qex, t);
+      this->exact(qex, t);
 
       double max = 0.0;
       for (size_t i = 0; i < qend.size(); i++) {
@@ -98,80 +209,17 @@ class AdvectionDiffusionSweeper
            << " (" << qend.size() << ", " << predict << ")"
            << endl;
 
-      errors.insert(pair<pair<size_t, size_t>, double>
-		    (pair<size_t, size_t>(n, k), max));
+      this->errors.insert(pair<pair<size_t, size_t>, double>(pair<size_t, size_t>(n, k), max));
     }
 
+    /**
+     * retrieve errors at iterations and time nodes
+     */
     error_map get_errors()
     {
-      return errors;
+      return this->errors;
     }
-
-    void predict(bool initial)
-    {
-      pfasst::encap::IMEXSweeper<time>::predict(initial);
-      time t  = this->get_controller()->get_time();
-      time dt = this->get_controller()->get_time_step();
-      echo_error(t + dt, true);
-    }
-
-    void sweep()
-    {
-      pfasst::encap::IMEXSweeper<time>::sweep();
-      time t  = this->get_controller()->get_time();
-      time dt = this->get_controller()->get_time_step();
-      echo_error(t + dt);
-    }
-
-    void f1eval(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/)
-    {
-      auto& q = as_vector<double,time>(_q);
-      auto& f = as_vector<double,time>(_f);
-
-      double c = -v / double(q.size());
-
-      auto* z = fft.forward(q);
-      for (size_t i = 0; i < q.size(); i++) {
-        z[i] *= c * ddx[i];
-      }
-      fft.backward(f);
-
-      nf1evals++;
-    }
-
-    void f2eval(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/)
-    {
-      auto& q = as_vector<double,time>(_q);
-      auto& f = as_vector<double,time>(_f);
-
-      double c = nu / double(q.size());
-
-      auto* z = fft.forward(q);
-      for (size_t i = 0; i < q.size(); i++) {
-        z[i] *= c * lap[i];
-      }
-      fft.backward(f);
-    }
-
-    void f2comp(shared_ptr<Encapsulation<time>> _f, shared_ptr<Encapsulation<time>> _q, time /*t*/, time dt,
-                shared_ptr<Encapsulation<time>> _rhs)
-    {
-      auto& q = as_vector<double,time>(_q);
-      auto& f = as_vector<double,time>(_f);
-      auto& rhs = as_vector<double,time>(_rhs);
-
-      auto* z = fft.forward(rhs);
-      for (size_t i = 0; i < q.size(); i++) {
-        z[i] /= (1.0 - nu * double(dt) * lap[i]) * double(q.size());
-      }
-      fft.backward(q);
-
-      for (size_t i = 0; i < q.size(); i++) {
-        f[i] = (q[i] - rhs[i]) / double(dt);
-      }
-
-    }
-
+    //! @}
 };
 
 #endif
diff --git a/examples/advection_diffusion/fft.hpp b/examples/advection_diffusion/fft.hpp
index cbf53f0f..faf24c9d 100644
--- a/examples/advection_diffusion/fft.hpp
+++ b/examples/advection_diffusion/fft.hpp
@@ -72,7 +72,6 @@ class FFT
         x[i] = real(wk->z[i]);
       }
     }
-
 };
 
 #endif
diff --git a/examples/advection_diffusion/mpi_pfasst.cpp b/examples/advection_diffusion/mpi_pfasst.cpp
new file mode 100644
index 00000000..13adae35
--- /dev/null
+++ b/examples/advection_diffusion/mpi_pfasst.cpp
@@ -0,0 +1,77 @@
+/*
+ * Advection/diffusion example using an encapsulated IMEX sweeper.
+ *
+ * This example uses MPI PFASST.
+ */
+
+#include <memory>
+#include <cassert>
+
+#include <mpi.h>
+#include <fftw3.h>
+
+#include <pfasst.hpp>
+#include <pfasst/pfasst.hpp>
+#include <pfasst/mpi_communicator.hpp>
+#include <pfasst/encap/automagic.hpp>
+#include <pfasst/encap/mpi_vector.hpp>
+
+#include "advection_diffusion_sweeper.hpp"
+#include "spectral_transfer_1d.hpp"
+
+using namespace std;
+using namespace pfasst;
+using namespace pfasst::encap;
+using namespace pfasst::mpi;
+
+error_map run_mpi_pfasst()
+{
+  const size_t nsteps = 4;
+  const double dt     = 0.01;
+  const size_t niters = 4;
+
+  vector<pair<size_t, pfasst::QuadratureType>> nodes = {
+    { 3, pfasst::QuadratureType::GaussLobatto },
+    { 5, pfasst::QuadratureType::GaussLobatto }
+  };
+
+  vector<size_t> ndofs = { 64, 128 };
+
+  auto build_level = [ndofs](size_t level) {
+    auto factory  = make_shared<MPIVectorFactory<double>>(ndofs[level]);
+    auto sweeper  = make_shared<AdvectionDiffusionSweeper<>>(ndofs[level]);
+    auto transfer = make_shared<SpectralTransfer1D<>>();
+
+    return AutoBuildTuple<>(sweeper, transfer, factory);
+  };
+
+  auto initial = [](shared_ptr<EncapSweeper<>> sweeper, shared_ptr<Encapsulation<>> q0) {
+    auto ad = dynamic_pointer_cast<AdvectionDiffusionSweeper<>>(sweeper);
+    assert(ad);
+    ad->exact(q0, 0.0);
+  };
+
+  MPICommunicator comm(MPI_COMM_WORLD);
+  PFASST<> pf;
+
+  auto_build(pf, nodes, build_level);
+  auto_setup(pf, initial);
+
+  pf.set_comm(&comm);
+  pf.set_duration(0.0, nsteps * dt, dt, niters);
+  pf.set_nsweeps({2, 1});
+  pf.run();
+
+  auto fine = pf.get_finest<AdvectionDiffusionSweeper<>>();
+  return fine->get_errors();
+}
+
+#ifndef PFASST_UNIT_TESTING
+int main(int argc, char** argv)
+{
+  MPI_Init(&argc, &argv);
+  run_mpi_pfasst();
+  fftw_cleanup();
+  MPI_Finalize();
+}
+#endif
diff --git a/examples/advection_diffusion/serial_mlsdc.cpp b/examples/advection_diffusion/serial_mlsdc.cpp
index fdf88906..f6c07ca9 100644
--- a/examples/advection_diffusion/serial_mlsdc.cpp
+++ b/examples/advection_diffusion/serial_mlsdc.cpp
@@ -39,7 +39,7 @@ error_map run_serial_mlsdc()
    * (according to 'xrat').
    */
   for (size_t l = 0; l < nlevs; l++) {
-    auto nodes    = compute_nodes<double>(nnodes, "gauss-lobatto");
+    auto nodes    = compute_nodes<double>(nnodes, pfasst::QuadratureType::GaussLobatto);
     auto factory  = make_shared<VectorFactory<double>>(ndofs);
     auto sweeper  = make_shared<AdvectionDiffusionSweeper<>>(ndofs);
     auto transfer = make_shared<SpectralTransfer1D<>>();
diff --git a/examples/advection_diffusion/serial_mlsdc_autobuild.cpp b/examples/advection_diffusion/serial_mlsdc_autobuild.cpp
index ff9ea71e..1a93aa65 100644
--- a/examples/advection_diffusion/serial_mlsdc_autobuild.cpp
+++ b/examples/advection_diffusion/serial_mlsdc_autobuild.cpp
@@ -7,8 +7,11 @@
  * controller.
  */
 
-#include <memory>
+#include <cstdlib>
 #include <cassert>
+#include <memory>
+#include <vector>
+#include <string>
 
 #include <fftw3.h>
 
@@ -24,6 +27,7 @@ using namespace std;
 using namespace pfasst;
 using namespace pfasst::encap;
 
+
 int main(int /*argc*/, char** /*argv*/)
 {
   MLSDC<> mlsdc;
@@ -32,9 +36,9 @@ int main(int /*argc*/, char** /*argv*/)
   const double dt     = 0.01;
   const size_t niters = 4;
 
-  vector<pair<size_t, string>> nodes = {
-    { 3, "gauss-lobatto" },
-    { 5, "gauss-lobatto" }
+  vector<pair<size_t, pfasst::QuadratureType>> nodes = {
+    { 3, pfasst::QuadratureType::GaussLobatto },
+    { 5, pfasst::QuadratureType::GaussLegendre }
   };
 
   vector<size_t> ndofs = { 64, 128 };
diff --git a/examples/advection_diffusion/spectral_transfer_1d.hpp b/examples/advection_diffusion/spectral_transfer_1d.hpp
index bc84c1c6..ae407966 100644
--- a/examples/advection_diffusion/spectral_transfer_1d.hpp
+++ b/examples/advection_diffusion/spectral_transfer_1d.hpp
@@ -6,6 +6,7 @@
 #define _SPECTRAL_TRANSFER_1D_HPP_
 
 #include <cassert>
+#include <cstdlib>
 #include <memory>
 
 #include <pfasst/encap/vector.hpp>
@@ -13,6 +14,9 @@
 
 #include "fft.hpp"
 
+using namespace std;
+
+
 template<typename time = pfasst::time_precision>
 class SpectralTransfer1D
   : public pfasst::encap::PolyInterpMixin<time>
@@ -22,14 +26,13 @@ class SpectralTransfer1D
     FFT fft;
 
   public:
-
-    void interpolate(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src)
+    void interpolate(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src) override
     {
-      auto& fine = as_vector<double,time>(dst);
-      auto& crse = as_vector<double,time>(src);
+      auto& fine = pfasst::encap::as_vector<double, time>(dst);
+      auto& crse = pfasst::encap::as_vector<double, time>(src);
 
-      auto* crse_z = fft.forward(crse);
-      auto* fine_z = fft.get_workspace(fine.size())->z;
+      auto* crse_z = this->fft.forward(crse);
+      auto* fine_z = this->fft.get_workspace(fine.size())->z;
 
       for (size_t i = 0; i < fine.size(); i++) {
         fine_z[i] = 0.0;
@@ -45,13 +48,13 @@ class SpectralTransfer1D
         fine_z[fine.size() - crse.size() / 2 + i] = c * crse_z[crse.size() / 2 + i];
       }
 
-      fft.backward(fine);
+      this->fft.backward(fine);
     }
 
-    void restrict(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src)
+    void restrict(shared_ptr<Encapsulation> dst, shared_ptr<const Encapsulation> src) override
     {
-      auto& fine = as_vector<double,time>(src);
-      auto& crse = as_vector<double,time>(dst);
+      auto& fine = pfasst::encap::as_vector<double, time>(src);
+      auto& crse = pfasst::encap::as_vector<double, time>(dst);
 
       size_t xrat = fine.size() / crse.size();
 
@@ -59,7 +62,6 @@ class SpectralTransfer1D
         crse[i] = fine[xrat*i];
       }
     }
-
 };
 
 #endif
diff --git a/examples/advection_diffusion/vanilla_sdc.cpp b/examples/advection_diffusion/vanilla_sdc.cpp
index de6544a9..3e663368 100644
--- a/examples/advection_diffusion/vanilla_sdc.cpp
+++ b/examples/advection_diffusion/vanilla_sdc.cpp
@@ -4,6 +4,7 @@
  * This example uses a vanilla SDC sweeper.
  */
 
+#include <cstdlib>
 #include <memory>
 
 #include <fftw3.h>
@@ -24,7 +25,7 @@ error_map run_vanilla_sdc()
   const size_t ndofs  = 64;
   const size_t niters = 4;
 
-  auto nodes   = pfasst::compute_nodes(nnodes, "gauss-lobatto");
+  auto nodes   = pfasst::compute_nodes(nnodes, pfasst::QuadratureType::GaussLobatto);
   auto factory = make_shared<pfasst::encap::VectorFactory<double>>(ndofs);
   auto sweeper = make_shared<AdvectionDiffusionSweeper<>>(ndofs);
 
diff --git a/examples/scalar/CMakeLists.txt b/examples/scalar/CMakeLists.txt
new file mode 100644
index 00000000..90df45b9
--- /dev/null
+++ b/examples/scalar/CMakeLists.txt
@@ -0,0 +1,39 @@
+set(all_example_runs ${all_example_runs})
+set(all_scalar_example_runs)
+
+message(STATUS "  scalar")
+include_directories(
+    ${pfasst_INCLUDES}
+    ${3rdparty_INCLUDES}
+)
+
+set(scalar_examples
+    scalar_sdc
+)
+
+foreach(example ${scalar_examples})
+    add_executable(${example} ${CMAKE_CURRENT_SOURCE_DIR}/${example}.cpp)
+    set_target_properties(${example}
+        PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/scalar
+    )
+endforeach(example)
+
+message(STATUS "    creating make targets to run examples")
+foreach(example ${scalar_examples})
+    add_custom_target(run_example_scalar_${example}
+        COMMAND ${example}
+        DEPENDS ${example}
+        WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/scalar
+        COMMENT "Running Scalar example: ${example}" VERBATIM
+    )
+    list(APPEND all_scalar_example_runs run_example_scalar_${example})
+    message(STATUS "      - run_example_scalar_${example}")
+endforeach(example)
+
+add_custom_target(run_example_scalar_all
+  DEPENDS ${all_scalar_example_runs}
+  COMMENT "Running all Scalar examples" VERBATIM
+)
+message(STATUS "    - run_example_scalar_all")
+
+set(all_example_runs ${all_example_runs} ${all_scalar_example_runs} PARENT_SCOPE)
diff --git a/examples/scalar/README.md b/examples/scalar/README.md
new file mode 100644
index 00000000..232ebf42
--- /dev/null
+++ b/examples/scalar/README.md
@@ -0,0 +1,7 @@
+# Scalar                                                                     {#page_examples_scalar}
+
+This directory contains a simple implementations of scalar ODE solver using the PFASST framework.
+
+The SDC sweeper is defined in `scalar_sweeper.hpp`.
+As this simple equation does not require any special handling, all the SDC magic is derived and 
+taken from pfasst::encap::IMEXSweeper::sweep().
diff --git a/examples/scalar/scalar_sdc.cpp b/examples/scalar/scalar_sdc.cpp
new file mode 100644
index 00000000..9dd5b42a
--- /dev/null
+++ b/examples/scalar/scalar_sdc.cpp
@@ -0,0 +1,69 @@
+/**
+ * Scalar test equation example using an encapsulated IMEX sweeper.
+ *
+ * Solves Dahlquist's test equation with single level SDC
+ *
+ * \\( u' = \\lambda*u \\quad\\text{ , } u(0) = u_0 \\)
+ *
+ * with complex lambda, treating the real part implicitly and the imaginary part 
+ * explicitly.
+ *
+ */
+
+#include<complex>
+
+#include<pfasst.hpp>
+#include<pfasst/sdc.hpp>
+#include<pfasst/encap/vector.hpp>
+
+#include "scalar_sweeper.hpp"
+
+double run_scalar_sdc(const size_t nsteps, const double dt, const size_t nnodes,
+                      const size_t niters, const complex<double> lambda, 
+                      const pfasst::QuadratureType nodetype)
+{
+  pfasst::SDC<> sdc;
+
+  // For test equation, set initial value to \\( 1+0i \\)
+  const complex<double> y0 = complex<double>(1.0, 0.0);
+
+  auto nodes = pfasst::compute_nodes(nnodes, nodetype);
+
+  // This is a scalar example, so we use the encap::VectorFactory with fixed length of 1 and 
+  //  complex type.
+  auto factory = make_shared<pfasst::encap::VectorFactory<complex<double>>>(1);
+  auto sweeper = make_shared<ScalarSweeper<>>(lambda, y0);
+
+  sweeper->set_nodes(nodes);
+  sweeper->set_factory(factory);
+
+  sdc.add_level(sweeper);
+
+  // Final time Tend = dt*nsteps
+  sdc.set_duration(0.0, dt*nsteps, dt, niters);
+  sdc.setup();
+
+  auto q0 = sweeper->get_state(0);
+  sweeper->exact(q0, 0.0);
+
+  sdc.run();
+
+  return sweeper->get_errors();
+}
+
+#ifndef PFASST_UNIT_TESTING
+/**
+ * Main routine running the scalar example with a preset parameters
+ */
+int main(int /*argc*/, char** /*argv*/)
+{
+  const size_t nsteps = 2;
+  const double dt     = 1.0;
+  const size_t nnodes = 4;
+  const size_t niters = 6;
+  const complex<double> lambda = complex<double>(-1.0, 1.0);
+  const pfasst::QuadratureType nodetype = pfasst::QuadratureType::GaussLobatto;
+
+  run_scalar_sdc(nsteps, dt, nnodes, niters, lambda, nodetype);
+}
+#endif
diff --git a/examples/scalar/scalar_sweeper.hpp b/examples/scalar/scalar_sweeper.hpp
new file mode 100644
index 00000000..af0a38b9
--- /dev/null
+++ b/examples/scalar/scalar_sweeper.hpp
@@ -0,0 +1,188 @@
+#ifndef _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
+#define _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
+
+#include <complex>
+#include <vector>
+
+#include <pfasst/encap/imex_sweeper.hpp>
+#include <pfasst/encap/vector.hpp>
+
+using namespace std;
+
+/**
+ * Sweeper for scalar test equation
+ *
+ * \\( u' = \\lambda*u \\quad\\text{ , } u(0) = u_0 \\)
+ *
+ * with complex lambda using an IMEX scheme. Derived from the generic imex_sweeper.
+ *
+ */
+template<typename time = pfasst::time_precision>
+class ScalarSweeper
+  : public pfasst::encap::IMEXSweeper<time>
+{
+  private:
+    typedef pfasst::encap::Encapsulation<time> encap_type;
+
+    //! Define a type for a complex PFASST vector encapsulation.
+    typedef pfasst::encap::VectorEncapsulation<complex<double>> complex_vector_type;
+
+     //! Parameter lambda and initial value \\( u_0 \\)
+    complex<double> lambda, u0;
+
+    //! The complex unit \\( i = \\sqrt{-1} \\)
+    const complex<double> i_complex = complex<double>(0, 1);
+
+     //! Error at the final time. For the scalar example, an analytical solution is known.
+    double error;
+
+     //! Counters for how often `f_expl_eval`, `f_impl_eval` and `impl_solve` are called.
+    size_t n_f_expl_eval, n_f_impl_eval, n_impl_solve;
+
+  public:
+    /**
+     * Generic constructor; initialize all function call counters with zero.
+     * @param[in] lambda coefficient in test equation
+     * @param[in] u0 initial value at \\( t=0 \\)
+     */ 
+    ScalarSweeper(const complex<double>& lambda, const complex<double>& u0)
+      :   lambda(lambda)
+        , u0(u0)
+        , error(0.0)
+        , n_f_expl_eval(0)
+        , n_f_impl_eval(0)
+        , n_impl_solve(0)
+    {}
+
+    /**
+     * Upon destruction, report final error and number of function calls
+     */
+    virtual ~ScalarSweeper()
+    {
+      cout << "Final error:                    " << scientific << this->error << endl;
+      cout << "Number of explicit evaluations: " << this->n_f_expl_eval << endl;
+      cout << "Number of implicit evaluations: " << this->n_f_impl_eval << endl;
+      cout << "Number of implicit solves:      " << this->n_impl_solve << endl;
+    }
+
+    /**
+     * Compute error between last state and exact solution at time tand print it to cout
+     *
+     * @param[in] t Time
+     */
+    void echo_error(time t)
+    {
+      auto& qend = pfasst::encap::as_vector<complex<double>, time>(this->get_end_state());
+
+      complex_vector_type qex(qend.size());
+
+      this->exact(qex, t);
+      double max_err = abs(qend[0] - qex[0]) / abs(qex[0]);
+      cout << "err: " << scientific << max_err << endl;
+      this->error = max_err;
+    }
+
+    /**
+     * Returns error, but does not update it!
+     */
+    double get_errors()
+    {
+      return this->error;
+    }
+
+    /**
+     * Prediction step and update of error. Uses predictor as provided by IMEXSweeper.
+     */
+    void predict(bool initial) override
+    {
+      pfasst::encap::IMEXSweeper<time>::predict(initial);
+      time t  = this->get_controller()->get_time();
+      time dt = this->get_controller()->get_time_step();
+      this->echo_error(t + dt);
+    }
+
+    /**
+     * Perform a sweep and update error. Uses sweep as provided by IMEXSweeper.
+     */
+    void sweep() override
+    {
+      pfasst::encap::IMEXSweeper<time>::sweep();
+      time t  = this->get_controller()->get_time();
+      time dt = this->get_controller()->get_time_step();
+      this->echo_error(t + dt);
+    }
+
+    /**
+     * Computes the exact solution \\( u_0 \\exp \\left( \\lambda*t \\right) \\) at a given time t.
+     *
+     * @param[in] t Time
+     */
+    void exact(complex_vector_type& q, time t)
+    {
+      q[0] = this->u0 * exp(this->lambda * t);
+    }
+
+    void exact(shared_ptr<encap_type> q_encap, time t)
+    {
+      auto& q = pfasst::encap::as_vector<complex<double>, time>(q_encap);
+      this->exact(q, t);
+    }
+
+    /**
+     * Evaluate the explicit part of the right hand side: Multiply with \\( \\text{imag}(\\lambda) \\)
+     */
+    void f_expl_eval(shared_ptr<encap_type> f_encap,
+                     shared_ptr<encap_type> q_encap, time t) override
+    {
+      UNUSED(t);
+      auto& f = pfasst::encap::as_vector<complex<double>, time>(f_encap);
+      auto& q = pfasst::encap::as_vector<complex<double>, time>(q_encap);
+
+      // f_expl = multiply with imaginary part of lambda
+      f[0] = this->i_complex * imag(this->lambda) * q[0];
+
+      this->n_f_expl_eval++;
+    }
+
+    /**
+     * Evaluate the implicit part of the right hand side: Multiply with \\( \\text{real}(\\lambda) \\)
+     */
+    void f_impl_eval(shared_ptr<encap_type> f_encap,
+                     shared_ptr<encap_type> q_encap, time t) override
+    {
+      UNUSED(t);
+      auto& f = pfasst::encap::as_vector<complex<double>, time>(f_encap);
+      auto& q = pfasst::encap::as_vector<complex<double>, time>(q_encap);
+
+      // f_impl = multiply with real part of lambda
+      f[0] = real(this->lambda) * q[0];
+
+      this->n_f_impl_eval++;
+    }
+
+    /**
+     * For given \\( b \\), solve 
+     * \\( \\left( \\mathbb{I}_d - \\Delta t \\text{real}(\\lambda) \\right) u = b \\) 
+     * for \\( u \\) and set f_encap to \\( \\text{real}(\\lambda) u \\)
+     */
+    void impl_solve(shared_ptr<encap_type> f_encap,
+                    shared_ptr<encap_type> q_encap, time t, time dt,
+                    shared_ptr<encap_type> rhs_encap) override
+    {
+      UNUSED(t);
+      auto& f = pfasst::encap::as_vector<complex<double>, time>(f_encap);
+      auto& q = pfasst::encap::as_vector<complex<double>, time>(q_encap);
+      auto& rhs = pfasst::encap::as_vector<complex<double>, time>(rhs_encap);
+
+      // invert f_impl = multiply with inverse of real part of lambda
+      double inv = 1.0 / (1.0 - double(dt) * real(this->lambda));
+      q[0] = inv * rhs[0];
+
+      // compute f_impl_eval of q[0], i.e. multiply with real(lambda)
+      f[0] = real(this->lambda) * q[0];
+
+      this->n_impl_solve++;
+    }
+};
+
+#endif  // _EXAMPLES__SCALAR__SCALAR_SWEEPER_HPP_
diff --git a/generate_coverage.sh b/generate_coverage.sh
new file mode 100755
index 00000000..00c35b2e
--- /dev/null
+++ b/generate_coverage.sh
@@ -0,0 +1,98 @@
+#!/bin/sh
+basepath="`pwd`"
+
+function print_help {
+  echo "#######################################################################"
+  echo "###               Generation of Test Coverage Report                ###"
+  echo "#                                                                     #"
+  echo "# This only works for builds made with GCC and the following CMake    #"
+  echo "# variables:                                                          #"
+  echo "#   -Dpfasst_WITH_GCC_PROF=ON -Dpfasst_BUILD_TESTS=ON                 #"
+  echo "#                                                                     #"
+  echo "# First (and only) parameter must be the name of the build directory  #"
+  echo "#                                                                     #"
+  echo "# Example:                                                            #"
+  echo "#   ./generate_coverage.sh build_gcc                                  #"
+  echo "#                                                                     #"
+  echo "#######################################################################"
+  return 0
+}
+
+
+if [[ $# -ne 1 ]]
+then
+  print_help
+  echo "ERROR: Please name the build directory as the first parameter."
+  exit -1
+fi
+
+builddir="${basepath}/${1}"
+if [[ ! -d "${builddir}" ]]
+then
+  echo "ERROR: Given build directory '${bilddir}' does not exist."
+  exit -1
+fi
+cd "${builddir}"
+
+coveragedir="${basepath}/coverage"
+rm -rf "${coveragedir}"
+mkdir -p "${coveragedir}"
+
+for testdir in `find ${builddir} -type d | grep -o 'tests/.*/.*\dir'`
+do
+  testname=`expr "${testdir}" : '.*\(test_[a-zA-Z\-_]*\)\.dir'`
+  echo "# ${testname}"
+  cd "${testdir}"
+  echo "    Deleting old tracing data"
+  lcov --zerocounters  --directory . \
+       > "${coveragedir}/log_${testname}_zero.log" 2&>1
+
+  cd "${builddir}"
+  echo "    Running test and generating tracing data"
+  ctest -R ${testname} \
+        > "${coveragedir}/log_${testname}_test.log"
+  if [[ $? -ne 0 ]]
+  then
+    echo "ERROR: test '${testname}' failed."
+    echo "       See '${coveragedir}/log_${testname}_test.log' for details."
+    cd ${basepath}
+    exit -1
+  fi
+
+  cd "${testdir}"
+  echo "    Capturing all tracing data"
+  lcov --directory . --capture --output-file "${testname}.info.tmp" \
+       > "${coveragedir}/log_${testname}_capture.log"
+
+  echo "    Extracting interesting tracing data"
+  lcov --extract "${testname}.info.tmp" "*${basepath}/include/**/*" \
+       --output-file "${testname}.info" \
+       > "${coveragedir}/log_${testname}_extract.log"
+  rm "${testname}.info.tmp"
+
+  cd "${builddir}"
+  echo "    Aggegrating coverage data"
+  if [[ -e "${coveragedir}/all_tests.info" ]]
+  then
+    lcov --add-tracefile "${coveragedir}/all_tests.info" \
+         --add-tracefile "${testdir}/${testname}.info" \
+         --output-file "${coveragedir}/all_tests.info" \
+         > "${coveragedir}/log_${testname}_aggregate.log"
+  else
+    lcov --add-tracefile "${testdir}/${testname}.info" \
+         --output-file "${coveragedir}/all_tests.info" \
+         > "${coveragedir}/log_${testname}_aggregate.log"
+  fi
+done
+
+cd ${basepath}
+echo "# Generating HTML report"
+genhtml --output-directory "${coveragedir}" \
+        --demangle-cpp --num-spaces 2 --sort \
+        --title "PFASST++ Test Coverage" --prefix "${basepath}/include" \
+        --function-coverage --branch-coverage --legend \
+        "${coveragedir}/all_tests.info" \
+        > "${coveragedir}/log_generate.log"
+
+echo "# Coverage report can be found in ${coveragedir}/index.html"
+echo "# Everything done."
diff --git a/include/pfasst/controller.hpp b/include/pfasst/controller.hpp
index b01bf678..6d410e44 100644
--- a/include/pfasst/controller.hpp
+++ b/include/pfasst/controller.hpp
@@ -9,6 +9,7 @@
 #include <memory>
 #include <cassert>
 #include <iterator>
+#include <iostream>
 
 #include "interfaces.hpp"
 
@@ -23,18 +24,22 @@ namespace pfasst
   class Controller
   {
     protected:
+      //! @{
       deque<shared_ptr<ISweeper<time>>>  levels;
       deque<shared_ptr<ITransfer<time>>> transfer;
+      //! @}
 
+      //! @{
       size_t step, iteration, max_iterations;
       time t, dt, tend;
+      //! @}
 
     public:
       //! @{
-      void setup()
+      virtual void setup()
       {
         for (auto l = coarsest(); l <= finest(); ++l) {
-	  l.current()->set_controller(this);
+          l.current()->set_controller(this);
           l.current()->setup();
         }
       }
@@ -74,13 +79,13 @@ namespace pfasst
       template<typename R = ISweeper<time>>
       shared_ptr<R> get_finest()
       {
-	return get_level<R>(nlevels()-1);
+        return get_level<R>(nlevels()-1);
       }
 
       template<typename R = ISweeper<time>>
       shared_ptr<R> get_coarsest()
       {
-	return get_level<R>(0);
+        return get_level<R>(0);
       }
 
       template<typename R = ITransfer<time>>
@@ -108,24 +113,28 @@ namespace pfasst
        * implementing a `RandomAccessIterator`.
        */
       class LevelIter
-        : iterator<random_access_iterator_tag, shared_ptr<ISweeper<time>>, size_t,
+        : iterator<random_access_iterator_tag, shared_ptr<ISweeper<time>>, int,
                    ISweeper<time>*, ISweeper<time>>
       {
           Controller* ts;
 
         public:
-          typedef size_t                     difference_type;
+          //! @{
+          typedef int                        difference_type;
           typedef shared_ptr<ISweeper<time>> value_type;
           typedef ISweeper<time>*            pointer;
           typedef ISweeper<time>             reference;
           typedef random_access_iterator_tag iterator_category;
+          //! @}
 
-          size_t level;
+          //! @{
+          int level;
+          //! @}
 
           //! @{
-          LevelIter(size_t level, Controller* ts)
+          LevelIter(int level, Controller* ts)
             : ts(ts), level(level)
-          {}
+          { }
           //! @}
 
           //! @{
@@ -178,7 +187,7 @@ namespace pfasst
       LevelIter coarsest() { return LevelIter(0, this); }
       //! @}
 
-
+      //! @{
       /**
        * Get current time step number.
        */
@@ -187,6 +196,12 @@ namespace pfasst
         return step;
       }
 
+      void set_step(size_t n)
+      {
+        t += (n - step)*dt;
+        step = n;
+      }
+
       time get_time_step()
       {
         return dt;
@@ -227,8 +242,9 @@ namespace pfasst
       {
         return max_iterations;
       }
-
+      //! @}
   };
-}
+
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/encap/automagic.hpp b/include/pfasst/encap/automagic.hpp
index daf5b8d0..2e07508c 100644
--- a/include/pfasst/encap/automagic.hpp
+++ b/include/pfasst/encap/automagic.hpp
@@ -16,12 +16,14 @@ namespace pfasst
   {
 
     template<typename time = time_precision>
-    using AutoBuildTuple = tuple<shared_ptr<pfasst::encap::EncapSweeper<time>>,
-                                 shared_ptr<pfasst::ITransfer<time>>,
-                                 shared_ptr<pfasst::encap::EncapFactory<time>>>;
+    using AutoBuildTuple = tuple<
+                             shared_ptr<pfasst::encap::EncapSweeper<time>>,
+                             shared_ptr<pfasst::ITransfer<time>>,
+                             shared_ptr<pfasst::encap::EncapFactory<time>>
+                           >;
 
     template<typename ControllerT, typename BuildT, typename time = time_precision>
-    void auto_build(ControllerT& c, vector<pair<size_t, string>> nodes, BuildT build)
+    void auto_build(ControllerT& c, vector<pair<size_t, pfasst::QuadratureType>> nodes, BuildT build)
     {
       for (size_t l = 0; l < nodes.size(); l++) {
         auto nds = pfasst::compute_nodes<time>(get<0>(nodes[l]), get<1>(nodes[l]));
@@ -48,7 +50,7 @@ namespace pfasst
       }
     }
 
-  }
-}
+  }  // ::pfasst::encap
+} // ::pfasst
 
 #endif
diff --git a/include/pfasst/encap/encap_sweeper.hpp b/include/pfasst/encap/encap_sweeper.hpp
index 1a73ad34..21b905de 100644
--- a/include/pfasst/encap/encap_sweeper.hpp
+++ b/include/pfasst/encap/encap_sweeper.hpp
@@ -5,8 +5,11 @@
 #ifndef _PFASST_ENCAP_ENCAP_SWEEPER_HPP_
 #define _PFASST_ENCAP_ENCAP_SWEEPER_HPP_
 
+#include <cstdlib>
 #include <vector>
+#include <memory>
 
+#include "../globals.hpp"
 #include "../interfaces.hpp"
 #include "../quadrature.hpp"
 #include "encapsulation.hpp"
@@ -20,50 +23,105 @@ namespace pfasst
     class EncapSweeper
       : public ISweeper<time>
     {
+      public:
+        //! @{
+        typedef Encapsulation<time> encap_type;
+        typedef EncapFactory<time> factory_type;
+        //! @}
+
+      private:
+        //! @{
         vector<time> nodes;
-        shared_ptr<EncapFactory<time>> factory;
+        vector<bool> is_proper;
+        shared_ptr<factory_type> factory;
+        //! @}
 
       public:
+        //! @{
+        virtual ~EncapSweeper()
+        {}
+        //! @}
 
-        void set_nodes(vector<time> nodes)
+        //! @{
+        virtual void set_nodes(vector<time> nodes)
         {
-          this->nodes = nodes;
+          auto augmented = pfasst::augment_nodes(nodes);
+          this->nodes = get<0>(augmented);
+          this->is_proper = get<1>(augmented);
         }
 
-        const vector<time> get_nodes() const
+        virtual const vector<time> get_nodes() const
         {
           return nodes;
         }
 
-        void set_factory(shared_ptr<EncapFactory<time>> factory)
+        const vector<bool> get_is_proper() const
         {
-          this->factory = shared_ptr<EncapFactory<time>>(factory);
+          return is_proper;
         }
 
-        shared_ptr<EncapFactory<time>> get_factory() const
+        virtual void set_factory(shared_ptr<factory_type> factory)
+        {
+          this->factory = factory;
+        }
+
+        virtual shared_ptr<factory_type> get_factory() const
         {
           return factory;
         }
 
-        virtual void set_state(shared_ptr<const Encapsulation<time>> /*q0*/, size_t /*m*/)
+        /**
+         * sets solution values at time node index `m`
+         *
+         * @param[in] u0 values to set
+         * @param[in] m 0-based node index
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void set_state(shared_ptr<const encap_type> u0, size_t m)
         {
+          UNUSED(u0); UNUSED(m);
           throw NotImplementedYet("sweeper");
         }
 
-        virtual shared_ptr<Encapsulation<time>> get_state(size_t /*m*/) const
+        /**
+         * retrieve solution values of current iteration at time node index `m`
+         *
+         * @param[in] m 0-based index of time node
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual shared_ptr<Encapsulation<time>> get_state(size_t m) const
         {
+          UNUSED(m);
           throw NotImplementedYet("sweeper");
           return NULL;
         }
 
-        virtual shared_ptr<Encapsulation<time>> get_tau(size_t /*m*/) const
+        /**
+         * retrieves FAS correction of current iteration at time node index `m`
+         *
+         * @param[in] m 0-based index of time node
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual shared_ptr<Encapsulation<time>> get_tau(size_t m) const
         {
+          UNUSED(m);
           throw NotImplementedYet("sweeper");
           return NULL;
         }
 
-        virtual shared_ptr<Encapsulation<time>> get_saved_state(size_t /*m*/) const
+        /**
+         * retrieves solution values of previous iteration at time node index `m`
+         *
+         * @param[in] m 0-based index of time node
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual shared_ptr<Encapsulation<time>> get_saved_state(size_t m) const
         {
+          UNUSED(m);
           throw NotImplementedYet("sweeper");
           return NULL;
         }
@@ -72,25 +130,103 @@ namespace pfasst
         {
           return this->get_state(this->get_nodes().size() - 1);
         }
-
-      virtual void evaluate(size_t /*m*/)
+        //! @}
+
+        //! @{
+        /**
+         * @copydoc ISweeper::advance()
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void advance() override
         {
           throw NotImplementedYet("sweeper");
         }
 
-        virtual void advance()
+        virtual void spread() override
         {
+          for (size_t m = 1; m < nodes.size(); m++) {
+            this->get_state(m)->copy(this->get_state(0));
+          }
+        }
+
+        /**
+         * evaluates the right hand side at given time node
+         *
+         * This evaluates the right hand side at the given time node with index `m` as returned by
+         * pfasst::encap::EncapSweeper::get_nodes:
+         *
+         * @param[in] m index of the time node to evaluate at
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void evaluate(size_t m)
+        {
+          UNUSED(m);
           throw NotImplementedYet("sweeper");
         }
 
-        virtual void integrate(time /*dt*/, vector<shared_ptr<Encapsulation<time>>> /*dst*/) const
+        /**
+         * integrates values of right hand side at all time nodes \\( t \\in [0,M-1] \\)
+         * simultaneously
+         *
+         * @param[in] dt width of the time interval to integrate
+         * @param[in,out] dst integrated values
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void integrate(time dt, vector<shared_ptr<encap_type>> dst) const
         {
+          UNUSED(dt); UNUSED(dst);
           throw NotImplementedYet("sweeper");
         }
+        //! @}
+
+        //! @{
+        virtual void post(ICommunicator* comm, int tag) override
+        {
+          this->get_state(0)->post(comm, tag);
+        }
+
+        virtual void send(ICommunicator* comm, int tag, bool blocking) override
+        {
+          this->get_state(this->get_nodes().size() - 1)->send(comm, tag, blocking);
+        }
+
+        virtual void recv(ICommunicator* comm, int tag, bool blocking) override
+        {
+          this->get_state(0)->recv(comm, tag, blocking);
+        }
+
+        virtual void broadcast(ICommunicator* comm) override
+        {
+          if (comm->rank() == comm->size() - 1) {
+            this->get_state(0)->copy(this->get_state(this->get_nodes().size() - 1));
+          }
+          this->get_state(0)->broadcast(comm);
+        }
+        //! @}
     };
 
-  }
 
-}
+    template<typename time>
+    EncapSweeper<time>& as_encap_sweeper(shared_ptr<ISweeper<time>> x)
+    {
+      shared_ptr<EncapSweeper<time>> y = dynamic_pointer_cast<EncapSweeper<time>>(x);
+      assert(y);
+      return *y.get();
+    }
+
+
+    template<typename time>
+    const EncapSweeper<time>& as_encap_sweeper(shared_ptr<const ISweeper<time>> x)
+    {
+      shared_ptr<const EncapSweeper<time>> y = dynamic_pointer_cast<const EncapSweeper<time>>(x);
+      assert(y);
+      return *y.get();
+    }
+
+  }  // ::pfasst::encap
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/encap/encapsulation.hpp b/include/pfasst/encap/encapsulation.hpp
index abf4df93..84419a68 100644
--- a/include/pfasst/encap/encapsulation.hpp
+++ b/include/pfasst/encap/encapsulation.hpp
@@ -8,6 +8,7 @@
 #include <vector>
 #include <memory>
 
+#include "../globals.hpp"
 #include "../interfaces.hpp"
 #include "../quadrature.hpp"
 
@@ -21,7 +22,7 @@ namespace pfasst
     typedef enum EncapType { solution, function } EncapType;
 
     /**
-     * basic encapsulation.
+     * Data/solution encapsulation.
      *
      * An Encapsulation provides basic functionality of the user data used by PFASST such as
      * mathematical operation @em axpy \\(y=ax+y\\) and packing/unpacking for message passing.
@@ -33,17 +34,32 @@ namespace pfasst
     {
       public:
         //! @{
-        virtual ~Encapsulation() { }
+        virtual ~Encapsulation()
+        {}
         //! @}
 
         //! @{
         // required for time-parallel communications
-        virtual void send()
+        virtual void post(ICommunicator* comm, int tag)
         {
+          UNUSED(comm); UNUSED(tag);
+        }
+
+        virtual void send(ICommunicator* comm, int tag, bool blocking)
+        {
+          UNUSED(comm); UNUSED(tag); UNUSED(blocking);
+          throw NotImplementedYet("pfasst");
+        }
+
+        virtual void recv(ICommunicator* comm, int tag, bool blocking)
+        {
+          UNUSED(comm); UNUSED(tag); UNUSED(blocking);
           throw NotImplementedYet("pfasst");
         }
-        virtual void recv()
+
+        virtual void broadcast(ICommunicator* comm)
         {
+          UNUSED(comm);
           throw NotImplementedYet("pfasst");
         }
         //! @}
@@ -68,18 +84,33 @@ namespace pfasst
          * Here, \\(a\\) is a time point and \\(x\\) another data structure (usually of the
          * same type).
          */
-        virtual void saxpy(time /*a*/, shared_ptr<const Encapsulation<time>> /*x*/)
+        virtual void saxpy(time a, shared_ptr<const Encapsulation<time>> x)
         {
+          UNUSED(a); UNUSED(x);
           throw NotImplementedYet("encap");
         }
+
         /**
          * defines matrix-vector multiplication for this data type.
          */
-        virtual void mat_apply(vector<shared_ptr<Encapsulation<time>>> /*dst*/, time /*a*/, matrix<time> /*m*/,
-			       vector<shared_ptr<Encapsulation<time>>> /*src*/, bool zero = true)
+        virtual void mat_apply(vector<shared_ptr<Encapsulation<time>>> dst, time a, matrix<time> mat,
+                               vector<shared_ptr<Encapsulation<time>>> src, bool zero = true)
         {
-	  (void) zero;
-          throw NotImplementedYet("encap");
+          size_t ndst = dst.size();
+          size_t nsrc = src.size();
+
+          if (zero) {
+            for (auto elem : dst) { elem->zero(); }
+          }
+
+          for (size_t n = 0; n < ndst; n++) {
+            for (size_t m = 0; m < nsrc; m++) {
+              auto s = mat(n, m);
+              if (s != 0.0) {
+                dst[n]->saxpy(a*s, src[m]);
+              }
+            }
+          }
         }
         //! @}
     };
@@ -91,7 +122,7 @@ namespace pfasst
         virtual shared_ptr<Encapsulation<time>> create(const EncapType) = 0;
     };
 
-  }
-}
+  }  // ::pfasst::encap
+} // ::pfasst
 
 #endif
diff --git a/include/pfasst/encap/imex_sweeper.hpp b/include/pfasst/encap/imex_sweeper.hpp
index 89df6447..1514cfc2 100644
--- a/include/pfasst/encap/imex_sweeper.hpp
+++ b/include/pfasst/encap/imex_sweeper.hpp
@@ -2,11 +2,14 @@
 #ifndef _PFASST_ENCAP_IMEX_SWEEPER_HPP_
 #define _PFASST_ENCAP_IMEX_SWEEPER_HPP_
 
+#include <cstdlib>
+#include <cassert>
 #include <vector>
 #include <memory>
 
 #include <boost/numeric/ublas/matrix.hpp>
 
+#include "../globals.hpp"
 #include "encapsulation.hpp"
 #include "encap_sweeper.hpp"
 
@@ -18,179 +21,407 @@ namespace pfasst
   {
     using pfasst::encap::Encapsulation;
 
+    /**
+     * interface for an semi-implicit sweeper
+     *
+     * Given an ODE \\( \\frac{\\partial}{\\partial t}u(t) = F(t,u) \\) where the function of the
+     * right hand side \\( F(t,u) \\) can be split into a non-stiff and a stiff part.
+     * To reduce complexity and computational effort one would want to solve the non-stiff part
+     * explicitly and the stiff part implicitly.
+     * Therefore, we define the splitting \\( F(t,u) = F_{expl}(t,u) + F_{impl}(t,u) \\).
+     *
+     * This sweeper requires three interfaces to be implement for such ODEs: two routines to
+     * evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a
+     * given state, and one that solves (perhaps with an external solver) the backward-Euler
+     * equation \\( U^{n+1} - \\Delta t F_{\\rm impl}(U^{n+1}) = RHS \\) for \\( U^{n+1} \\).
+     *
+     * @tparam time precision type of the time dimension
+     */
     template<typename time = time_precision>
     class IMEXSweeper
       : public pfasst::encap::EncapSweeper<time>
     {
-        vector<shared_ptr<Encapsulation<time>>> Q, pQ, S, T, Fe, Fi;
-        matrix<time> Smat, SEmat, SImat;
-
-      public:
-        ~IMEXSweeper()
-        {}
-
-        void set_state(shared_ptr<const Encapsulation<time>> q0, size_t m)
+      protected:
+        //! @{
+        /**
+         * solution values \\( u(\\tau) \\) at all time nodes \\( \\tau \\in [0, M-1] \\) of the
+         * current iteration
+         */
+        vector<shared_ptr<Encapsulation<time>>> u_state;
+
+        /**
+         * solution values \\( u(t) \\) at all time nodes \\( t \\in [0, M-1] \\) of the
+         * previous iteration
+         */
+        vector<shared_ptr<Encapsulation<time>>> previous_u_state;
+
+        /**
+         * node-to-node integrated values of \\( F(t,u) \\) at all time nodes \\( t \\in
+         * [0, M-1] \\) of the current iteration
+         */
+        vector<shared_ptr<Encapsulation<time>>> s_integrals;
+
+        /**
+         * FAS corrections \\( \\tau_t \\) at all time nodes \\( t \\in [0, M-1] \\) of the current
+         * iteration
+         */
+        vector<shared_ptr<Encapsulation<time>>> fas_corrections;
+
+        /**
+         * values of the explicit part of the right hand side \\( F_{expl}(t,u) \\) at all time
+         * nodes \\( t \\in [0, M-1] \\) of the current iteration
+         */
+        vector<shared_ptr<Encapsulation<time>>> fs_expl;
+
+        /**
+         * values of the explicit part of the right hand side \\( F_{impl}(t,u) \\) at all time
+         * nodes \\( t \\in [0, M-1] \\) of the current iteration
+         */
+        vector<shared_ptr<Encapsulation<time>>> fs_impl;
+        //! @}
+
+        //! @{
+        /**
+         * quadrature matrix containing weights for node-to-node integration
+         */
+        matrix<time> s_mat;
+
+        /**
+         * quadrature matrix containing weights for 0-to-last node integration
+         */
+        matrix<time> b_mat;
+
+        /**
+         * quadrature matrix containing weights for node-to-node integration of explicit part
+         *
+         * @see IMEXSweeper::setup(bool) for a short description
+         */
+        matrix<time> s_mat_expl;
+
+        /**
+         * quadrature matrix containing weights for node-to-node integration of implicit part
+         *
+         * @see IMEXSweeper::setup(bool) for a short description
+         */
+        matrix<time> s_mat_impl;
+        //! @}
+
+
+        /**
+        * Set end state to \\( U_0 + \\int F_{expl} + F_{expl} \\).
+        */
+        void integrate_end_state(time dt)
         {
-          Q[m]->copy(q0);
+          vector<shared_ptr<Encapsulation<time>>> dst = { this->u_state.back() };
+          dst[0]->copy(this->u_state.front());
+          dst[0]->mat_apply(dst, dt, this->b_mat, this->fs_expl, false);
+          dst[0]->mat_apply(dst, dt, this->b_mat, this->fs_impl, false);
         }
 
-        shared_ptr<Encapsulation<time>> get_state(size_t m) const
+        inline bool last_node_is_virtual()
         {
-          return Q[m];
+          return !this->get_is_proper().back();
         }
 
-        shared_ptr<Encapsulation<time>> get_tau(size_t m) const
+
+      public:
+        //! @{
+        virtual ~IMEXSweeper()
+        {}
+        //! @}
+
+        //! @{
+        virtual void set_state(shared_ptr<const Encapsulation<time>> u0, size_t m) override
         {
-          return T[m];
+          this->u_state[m]->copy(u0);
         }
 
-        shared_ptr<Encapsulation<time>> get_saved_state(size_t m) const
+        virtual shared_ptr<Encapsulation<time>> get_state(size_t m) const override
         {
-          return pQ[m];
+          return this->u_state[m];
         }
 
-        virtual void advance()
+        virtual shared_ptr<Encapsulation<time>> get_tau(size_t m) const override
         {
-          Q[0]->copy(Q.back());
-          Fe[0]->copy(Fe.back());
-          Fi[0]->copy(Fi.back());
+          return this->fas_corrections[m];
         }
 
-        virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const
+        virtual shared_ptr<Encapsulation<time>> get_saved_state(size_t m) const override
         {
-          dst[0]->mat_apply(dst, dt, Smat, Fe, true);
-          dst[0]->mat_apply(dst, dt, Smat, Fi, false);
+          return this->previous_u_state[m];
         }
-
-        void setup(bool coarse)
+        //! @}
+
+        //! @{
+        /**
+         * @copydoc ISweeper::setup(bool)
+         *
+         * To reduce computational overhead, we precompute the partial integration matrices
+         * IMEXSweeper::s_mat_expl \\( (\\tilde{s}^{expl})_{m,j} \\) and IMEXSweeper::s_mat_impl
+         * \\( (\\tilde{s}^{impl})_{m,j} \\) by incorporating known values of the SDC sweep
+         * equation.
+         *
+         * Let \\( F = F_{impl} + F_{expl} \\), \\( f = F_{impl} \\), \\( g = F_{expl} \\) and
+         * \\( \\Delta t_m = t_{m+1} - t_m \\).
+         * @f{eqnarray*}{
+         *   u_{m+1}^{k+1} &=& u_m^k + \Delta t_m \left( f_{m+1}^{k+1} - f_{m+1}^k \right)
+         *                     + \Delta t_m \left( g_m^{k+1} - g_m^k \right)
+         *                     + \sum_{j=1}^M s_{m,j} F_j^k \\
+         *                 &=& u_m^k + \Delta t_m f_{m+1}^{k+1} + \Delta t_m g_m^{k+1}
+         *                     + \sum_{j=1}^M s_{m,j} F_j^k
+         *                     - \Delta t_m \left( f_{m+1}^k - g_m^k \right) \\
+         *                 &=& u_m^k + \Delta t_m f_{m+1}^{k+1} + \Delta t_m g_m^{k+1}
+         *                     + \sum_{j=1}^M s_{m,j} f_j^k - \Delta t_m f_{m+1}^k
+         *                     + \sum_{j=1}^M s_{m,j} g_j^k - \Delta t_m g_m^k \\
+         *                 &=& u_m^k + \Delta t_m f_{m+1}^{k+1} + \Delta t_m g_m^{k+1}
+         *                     + \sum_{j=1}^M \tilde{s}_{m,j}^{impl} f_j^k
+         *                     + \sum_{j=1}^M \tilde{s}_{m,j}^{expl} g_j^k
+         * @f}
+         * with
+         * @f[
+         *  \tilde{s}_{m,j}^{impl} =
+         *    \begin{cases}
+         *      s_{m,j} - \Delta t_m &\mbox{if } j \equiv m+1 \\
+         *      s_{m,j} &\mbox{else}
+         *    \end{cases} \\
+         *  \tilde{s}_{m,j}^{expl} =
+         *    \begin{cases}
+         *      s_{m,j} - \Delta t_m &\mbox{if } j \equiv m \\
+         *      s_{m,j} &\mbox{else}
+         *    \end{cases} \\
+         * @f]
+         *
+         * This procedure can be derived from the so called \\( Q_{\\Delta} \\) notation.
+         */
+        virtual void setup(bool coarse) override
         {
           auto nodes = this->get_nodes();
+          auto is_proper = this->get_is_proper();
           assert(nodes.size() >= 1);
 
-          Smat = compute_quadrature(nodes, nodes, 's');
+          this->s_mat = compute_quadrature(nodes, nodes, is_proper, QuadratureMatrix::S);
 
-          SEmat = Smat;
-          SImat = Smat;
+          auto q_mat = compute_quadrature(nodes, nodes, is_proper, QuadratureMatrix::Q);
+          this->b_mat = matrix<time>(1, nodes.size());
+          for (size_t m = 0; m < nodes.size(); m++) {
+            this->b_mat(0, m) = q_mat(nodes.size() - 2, m);
+          }
+
+          this->s_mat_expl = this->s_mat;
+          this->s_mat_impl = this->s_mat;
           for (size_t m = 0; m < nodes.size() - 1; m++) {
             time ds = nodes[m + 1] - nodes[m];
-            SEmat(m, m)     -= ds;
-            SImat(m, m + 1) -= ds;
+            this->s_mat_expl(m, m)     -= ds;
+            this->s_mat_impl(m, m + 1) -= ds;
           }
 
           for (size_t m = 0; m < nodes.size(); m++) {
-            Q.push_back(this->get_factory()->create(pfasst::encap::solution));
+            this->u_state.push_back(this->get_factory()->create(pfasst::encap::solution));
             if (coarse) {
-              pQ.push_back(this->get_factory()->create(pfasst::encap::solution));
+              this->previous_u_state.push_back(this->get_factory()->create(pfasst::encap::solution));
             }
-            Fe.push_back(this->get_factory()->create(pfasst::encap::function));
-            Fi.push_back(this->get_factory()->create(pfasst::encap::function));
+            // XXX: can prolly save some f's here for virtual nodes...
+            this->fs_expl.push_back(this->get_factory()->create(pfasst::encap::function));
+            this->fs_impl.push_back(this->get_factory()->create(pfasst::encap::function));
           }
 
           for (size_t m = 0; m < nodes.size() - 1; m++) {
-            S.push_back(this->get_factory()->create(pfasst::encap::solution));
+            this->s_integrals.push_back(this->get_factory()->create(pfasst::encap::solution));
             if (coarse) {
-              T.push_back(this->get_factory()->create(pfasst::encap::solution));
+              this->fas_corrections.push_back(this->get_factory()->create(pfasst::encap::solution));
             }
           }
         }
 
-        virtual void sweep()
+        /**
+         * Compute low-order provisional solution.
+         *
+         * This does not simply copy the initial value to all time nodes but carries out a few
+         * forward/backward IMEX Euler steps between the nodes.
+         *
+         * @param[in] initial if `true` the explicit and implicit part of the right hand side of the
+         *     ODE get evaluated with the initial value
+         */
+        virtual void predict(bool initial) override
         {
           const auto   nodes  = this->get_nodes();
           const size_t nnodes = nodes.size();
           assert(nnodes >= 1);
 
-	  time dt = this->get_controller()->get_time_step();
+          time dt = this->get_controller()->get_time_step();
           time t  = this->get_controller()->get_time();
 
-          // integrate
-          S[0]->mat_apply(S, dt, SEmat, Fe, true);
-          S[0]->mat_apply(S, dt, SImat, Fi, false);
-          if (T.size() > 0) {
-            for (size_t m = 0; m < nnodes - 1; m++) {
-              S[m]->saxpy(1.0, T[m]);
-            }
-          }
+          if (initial) { this->evaluate(0); }
 
-          // sweep
           shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
 
-          for (size_t m = 0; m < nnodes - 1; m++) {
+          for (size_t m = 0; m < nnodes - (this->last_node_is_virtual() ? 2 : 1); m++) {
             time ds = dt * (nodes[m + 1] - nodes[m]);
-
-            rhs->copy(Q[m]);
-            rhs->saxpy(ds, Fe[m]);
-            rhs->saxpy(1.0, S[m]);
-            f2comp(Fi[m + 1], Q[m + 1], t, ds, rhs);
-            f1eval(Fe[m + 1], Q[m + 1], t + ds);
+            rhs->copy(this->u_state[m]);
+            rhs->saxpy(ds, this->fs_expl[m]);
+            this->impl_solve(this->fs_impl[m + 1], this->u_state[m + 1], t, ds, rhs);
+            this->f_expl_eval(this->fs_expl[m + 1], this->u_state[m + 1], t + ds);
 
             t += ds;
           }
+
+          if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
         }
 
-        virtual void predict(bool initial)
+        virtual void sweep() override
         {
           const auto   nodes  = this->get_nodes();
           const size_t nnodes = nodes.size();
           assert(nnodes >= 1);
 
-	  time dt = this->get_controller()->get_time_step();
+          time dt = this->get_controller()->get_time_step();
           time t  = this->get_controller()->get_time();
 
-          if (initial) {
-            f1eval(Fe[0], Q[0], t);
-            f2eval(Fi[0], Q[0], t);
+          // integrate
+          this->s_integrals[0]->mat_apply(this->s_integrals, dt, this->s_mat_expl, this->fs_expl, true);
+          this->s_integrals[0]->mat_apply(this->s_integrals, dt, this->s_mat_impl, this->fs_impl, false);
+          if (this->fas_corrections.size() > 0) {
+            for (size_t m = 0; m < nnodes - 1; m++) {
+              this->s_integrals[m]->saxpy(1.0, this->fas_corrections[m]);
+            }
           }
 
+          // sweep
           shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
 
-          for (size_t m = 0; m < nnodes - 1; m++) {
+          for (size_t m = 0; m < nnodes - (this->last_node_is_virtual() ? 2 : 1); m++) {
             time ds = dt * (nodes[m + 1] - nodes[m]);
-            rhs->copy(Q[m]);
-            rhs->saxpy(ds, Fe[m]);
-            f2comp(Fi[m + 1], Q[m + 1], t, ds, rhs);
-            f1eval(Fe[m + 1], Q[m + 1], t + ds);
+
+            rhs->copy(this->u_state[m]);
+            rhs->saxpy(ds, this->fs_expl[m]);
+            rhs->saxpy(1.0, this->s_integrals[m]);
+            this->impl_solve(this->fs_impl[m + 1], this->u_state[m + 1], t, ds, rhs);
+            this->f_expl_eval(this->fs_expl[m + 1], this->u_state[m + 1], t + ds);
 
             t += ds;
           }
+
+          if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
         }
 
-        virtual void save()
+        virtual void advance() override
         {
-          for (size_t m = 0; m < pQ.size(); m++) {
-            pQ[m]->copy(Q[m]);
+          this->u_state[0]->copy(this->u_state.back());
+          this->fs_expl[0]->copy(this->fs_expl.back());
+          this->fs_impl[0]->copy(this->fs_impl.back());
+        }
+
+        virtual void save(bool initial_only) override
+        {
+          if (initial_only) {
+            this->previous_u_state[0]->copy(u_state[0]);
+          } else {
+            for (size_t m = 0; m < this->previous_u_state.size(); m++) {
+              this->previous_u_state[m]->copy(u_state[m]);
+            }
           }
         }
 
-        virtual void evaluate(size_t m)
+        /**
+         * @copybrief EncapSweeper::evaluate()
+         *
+        * If the node `m` is virtual (not proper), we can save some
+        * implicit/explicit evaluations.
+        */
+        virtual void evaluate(size_t m) override
         {
-          time t = this->get_nodes()[m]; // XXX
-          f1eval(Fe[m], Q[m], t);
-          f2eval(Fi[m], Q[m], t);
+          time t0 = this->get_controller()->get_time();
+          time dt = this->get_controller()->get_time_step();
+          time t =  t0 + dt * this->get_nodes()[m];
+          if ((m == 0) || this->get_is_proper()[m]) {
+            this->f_expl_eval(this->fs_expl[m], this->u_state[m], t);
+          }
+          if (this->get_is_proper()[m]) {
+            this->f_impl_eval(this->fs_impl[m], this->u_state[m], t);
+          }
         }
 
-        virtual void f1eval(shared_ptr<Encapsulation<time>> /*f*/, shared_ptr<Encapsulation<time>> /*q*/,
-			    time /*t*/)
+        /**
+         * @copybrief EncapSweeper::integrate()
+         *
+         * @param[in] dt width of time interval to integrate over
+         * @param[in,out] dst integrated values; will get zeroed out beforehand
+         */
+        virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const override
+        {
+          dst[0]->mat_apply(dst, dt, this->s_mat, this->fs_expl, true);
+          dst[0]->mat_apply(dst, dt, this->s_mat, this->fs_impl, false);
+        }
+        //! @}
+
+        //! @{
+        /**
+         * Evaluates the explicit part of the right hand side of the ODE at the given time.
+         *
+         * @param[in,out] f_expl_encap Encapsulation to store the evaluated right hand side
+         * @param[in] u_encap Encapsulation storing the solution values to use for computing the
+         *     explicit part of the right hand side of the ODE
+         * @param[in] t time point of the evaluation
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void f_expl_eval(shared_ptr<Encapsulation<time>> f_expl_encap,
+                                 shared_ptr<Encapsulation<time>> u_encap,
+                                 time t)
         {
-          throw NotImplementedYet("imex (f1eval)");
+          UNUSED(f_expl_encap); UNUSED(u_encap); UNUSED(t);
+          throw NotImplementedYet("imex (f_expl_eval)");
         }
 
-        virtual void f2eval(shared_ptr<Encapsulation<time>> /*f*/, shared_ptr<Encapsulation<time>> /*q*/,
-			    time /*t*/)
+        /**
+         * Evaluates the implicit part of the right hand side of the ODE at the given time.
+         *
+         * This is typically called to compute the implicit part of the right hand side at the first
+         * collocation node, and on all nodes after restriction or interpolation.
+         *
+         * @param[in,out] f_impl_encap Encapsulation to store the evaluated right hand side
+         * @param[in] u_encap Encapsulation storing the solution values to use for computing the
+         *     implicit part of the right hand side of the ODE
+         * @param[in] t time point of the evaluation
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void f_impl_eval(shared_ptr<Encapsulation<time>> f_impl_encap,
+                                 shared_ptr<Encapsulation<time>> u_encap,
+                                 time t)
         {
-          throw NotImplementedYet("imex (f2eval)");
+          UNUSED(f_impl_encap); UNUSED(u_encap); UNUSED(t);
+          throw NotImplementedYet("imex (f_impl_eval)");
         }
 
-        virtual void f2comp(shared_ptr<Encapsulation<time>> /*f*/, shared_ptr<Encapsulation<time>> /*q*/,
-                            time /*t*/, time /*dt*/,
-                            shared_ptr<Encapsulation<time>> /*rhs*/)
+        /**
+         * Solves \\( U - \\Delta t f_{\\rm impl}(U) = RHS \\) for \\( U \\).
+         *
+         * During an IMEX SDC sweep, the correction equation is evolved using a forward-Euler
+         * stepper for the explicit piece, and a backward-Euler stepper for the implicit piece.
+         * This routine (implemented by the user) performs the solve required to perform one
+         * backward-Euler sub-step, and also returns \\( f_{\\rm impl}(U) \\).
+         *
+         * @param[in,out] f_encap Encapsulation to store the evaluated right hand side
+         * @param[in,out] u_encap Encapsulation to store the solution of the backward-Euler sub-step
+         * @param[in] t time point (of \\( RHS \\))
+         * @param[in] dt sub-step size to the previous time point (\\( \\Delta t \\))
+         * @param[in] rhs_encap Encapsulation storing \\( RHS \\)
+         *
+         * @note This method must be implemented in derived sweepers.
+         */
+        virtual void impl_solve(shared_ptr<Encapsulation<time>> f_encap,
+                                shared_ptr<Encapsulation<time>> u_encap,
+                                time t, time dt,
+                                shared_ptr<Encapsulation<time>> rhs_encap)
         {
+          UNUSED(f_encap); UNUSED(u_encap); UNUSED(t); UNUSED(dt); UNUSED(rhs_encap);
           throw NotImplementedYet("imex (f2comp)");
         }
-
+        //! @}
     };
 
-  }
-}
+  }  // ::pfasst::encap
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/encap/mpi_vector.hpp b/include/pfasst/encap/mpi_vector.hpp
new file mode 100644
index 00000000..93bfae18
--- /dev/null
+++ b/include/pfasst/encap/mpi_vector.hpp
@@ -0,0 +1,143 @@
+/*
+ * MPI enabled vector.
+ */
+
+#ifndef _PFASST_MPI_VECTOR_HPP_
+#define _PFASST_MPI_VECTOR_HPP_
+
+#include <cassert>
+
+#include <mpi.h>
+
+#include "vector.hpp"
+#include "../mpi_communicator.hpp"
+
+using namespace std;
+
+namespace pfasst
+{
+  namespace encap
+  {
+
+    using namespace pfasst::mpi;
+
+    template<typename scalar, typename time = time_precision>
+    class MPIVectorEncapsulation
+      : public VectorEncapsulation<scalar, time>
+    {
+        //! @{
+        MPI_Request recv_request = MPI_REQUEST_NULL;
+        MPI_Request send_request = MPI_REQUEST_NULL;
+        //! @}
+
+        //! @{
+        inline MPICommunicator& as_mpi(ICommunicator* comm)
+        {
+          auto mpi = dynamic_cast<MPICommunicator*>(comm);
+          assert(mpi);
+          return *mpi;
+        }
+        //! @}
+
+      public:
+        //! @{
+        MPIVectorEncapsulation(int size)
+          : VectorEncapsulation<scalar, time>(size)
+        {}
+
+        virtual ~MPIVectorEncapsulation()
+        {}
+        //! @}
+
+        //! @{
+        void post(ICommunicator* comm, int tag) override
+        {
+          auto& mpi = as_mpi(comm);
+          if (mpi.size() == 1) { return; }
+          if (mpi.rank() == 0) { return; }
+
+          int err = MPI_Irecv(this->data(), sizeof(scalar) * this->size(), MPI_CHAR,
+                              (mpi.rank() - 1) % mpi.size(), tag, mpi.comm, &recv_request);
+          if (err != MPI_SUCCESS) {
+            throw MPIError();
+          }
+        }
+
+        virtual void recv(ICommunicator* comm, int tag, bool blocking) override
+        {
+          auto& mpi = as_mpi(comm);
+          if (mpi.size() == 1) { return; }
+          if (mpi.rank() == 0) { return; }
+
+          int err;
+          if (blocking) {
+            MPI_Status stat;
+            err = MPI_Recv(this->data(), sizeof(scalar) * this->size(), MPI_CHAR,
+                           (mpi.rank() - 1) % mpi.size(), tag, mpi.comm, &stat);
+          } else {
+            MPI_Status stat;
+            err = MPI_Wait(&recv_request, &stat);
+          }
+
+          if (err != MPI_SUCCESS) {
+            throw MPIError();
+          }
+        }
+
+        virtual void send(ICommunicator* comm, int tag, bool blocking) override
+        {
+          auto& mpi = as_mpi(comm);
+          if (mpi.size() == 1) { return; }
+          if (mpi.rank() == mpi.size() - 1) { return; }
+
+          int err = MPI_SUCCESS;
+          if (blocking) {
+            err = MPI_Send(this->data(), sizeof(scalar) * this->size(), MPI_CHAR,
+                           (mpi.rank() + 1) % mpi.size(), tag, mpi.comm);
+          } else {
+            MPI_Status stat;
+            err = MPI_Wait(&send_request, &stat);
+            if (err != MPI_SUCCESS) {
+              throw MPIError();
+            }
+
+            err = MPI_Isend(this->data(), sizeof(scalar) * this->size(), MPI_CHAR,
+                            (mpi.rank() + 1) % mpi.size(), tag, mpi.comm, &send_request);
+          }
+
+          if (err != MPI_SUCCESS) {
+            throw MPIError();
+          }
+        }
+
+        virtual void broadcast(ICommunicator* comm) override
+        {
+          auto& mpi = as_mpi(comm);
+          int err = MPI_Bcast(this->data(), sizeof(scalar) * this->size(), MPI_CHAR,
+                              comm->size()-1, mpi.comm);
+
+          if (err != MPI_SUCCESS) {
+            throw MPIError();
+          }
+        }
+        //! @}
+    };
+
+    template<typename scalar, typename time = time_precision>
+    class MPIVectorFactory
+      : public EncapFactory<time>
+    {
+        int size;
+      public:
+        int dofs() { return size; }
+        MPIVectorFactory(const int size) : size(size) { }
+        shared_ptr<Encapsulation<time>> create(const EncapType)
+        {
+          return make_shared<MPIVectorEncapsulation<scalar, time>>(size);
+        }
+    };
+
+  }  // ::pfasst::encap
+}  // ::pfasst
+
+#endif
diff --git a/include/pfasst/encap/poly_interp.hpp b/include/pfasst/encap/poly_interp.hpp
index 343bf5d9..acb92238 100644
--- a/include/pfasst/encap/poly_interp.hpp
+++ b/include/pfasst/encap/poly_interp.hpp
@@ -9,6 +9,7 @@
 #include <cassert>
 #include <memory>
 
+#include "../globals.hpp"
 #include "../interfaces.hpp"
 #include "encap_sweeper.hpp"
 
@@ -21,60 +22,67 @@ namespace pfasst
     class PolyInterpMixin
       : public pfasst::ITransfer<time>
     {
-        using EncapVecT = vector<shared_ptr<Encapsulation<time>>>;
-        matrix<time> tmat, fmat;
+        //! @{
+        typedef vector<shared_ptr<Encapsulation<time>>> EncapVecT;
+        matrix<time> tmat;
+        matrix<time> fmat;
+        //! @}
 
       public:
-        virtual ~PolyInterpMixin() { }
-
-        virtual void interpolate(shared_ptr<ISweeper<time>> dst,
-                                 shared_ptr<const ISweeper<time>> src,
-                                 bool interp_delta_from_initial,
-                                 bool interp_initial)
+        //! @{
+        virtual ~PolyInterpMixin()
+        {}
+        //! @}
+
+        //! @{
+        virtual void interpolate_initial(shared_ptr<ISweeper<time>> dst,
+                                         shared_ptr<const ISweeper<time>> src) override
         {
-          shared_ptr<EncapSweeper<time>> fine = dynamic_pointer_cast<EncapSweeper<time>>(dst);
-          assert(fine);
-          shared_ptr<const EncapSweeper<time>> crse = dynamic_pointer_cast<const EncapSweeper<time>>(src);
-          assert(crse);
+          auto& fine = as_encap_sweeper(dst);
+          auto& crse = as_encap_sweeper(src);
+
+          auto crse_factory = crse.get_factory();
+          auto fine_factory = fine.get_factory();
+
+          auto crse_delta = crse_factory->create(solution);
+          this->restrict(crse_delta, fine.get_state(0));
+          crse_delta->saxpy(-1.0, crse.get_state(0));
+
+          auto fine_delta = fine_factory->create(solution);
+          this->interpolate(fine_delta, crse_delta);
+          fine.get_state(0)->saxpy(-1.0, fine_delta);
 
-          this->interpolate(fine, crse, interp_delta_from_initial, interp_initial);
+          fine.evaluate(0);
         }
 
-        virtual void interpolate(shared_ptr<EncapSweeper<time>> fine,
-                                 shared_ptr<const EncapSweeper<time>> crse,
-                                 bool interp_delta_from_initial,
-                                 bool interp_initial)
+
+        virtual void interpolate(shared_ptr<ISweeper<time>> dst,
+                                 shared_ptr<const ISweeper<time>> src,
+                                 bool interp_initial) override
         {
+          auto& fine = as_encap_sweeper(dst);
+          auto& crse = as_encap_sweeper(src);
+
           if (tmat.size1() == 0) {
-            tmat = pfasst::compute_interp<time>(fine->get_nodes(), crse->get_nodes());
+            tmat = pfasst::compute_interp<time>(fine.get_nodes(), crse.get_nodes());
           }
 
-          size_t nfine = fine->get_nodes().size();
-          size_t ncrse = crse->get_nodes().size();
+          size_t nfine = fine.get_nodes().size();
+          size_t ncrse = crse.get_nodes().size();
 
-          auto crse_factory = crse->get_factory();
-          auto fine_factory = fine->get_factory();
+          auto crse_factory = crse.get_factory();
+          auto fine_factory = fine.get_factory();
 
           EncapVecT fine_state(nfine), fine_delta(ncrse);
 
-          for (size_t m = 0; m < nfine; m++) { fine_state[m] = fine->get_state(m); }
+          for (size_t m = 0; m < nfine; m++) { fine_state[m] = fine.get_state(m); }
           for (size_t m = 0; m < ncrse; m++) { fine_delta[m] = fine_factory->create(solution); }
 
-          if (interp_delta_from_initial) {
-            for (size_t m = 1; m < nfine; m++) {
-              fine_state[m]->copy(fine_state[0]);
-            }
-          }
-
           auto crse_delta = crse_factory->create(solution);
           size_t m0 = interp_initial ? 0 : 1;
           for (size_t m = m0; m < ncrse; m++) {
-            crse_delta->copy(crse->get_state(m));
-            if (interp_delta_from_initial) {
-              crse_delta->saxpy(-1.0, crse->get_state(0));
-            } else {
-              crse_delta->saxpy(-1.0, crse->get_saved_state(m));
-            }
+            crse_delta->copy(crse.get_state(m));
+            crse_delta->saxpy(-1.0, crse.get_saved_state(m));
             interpolate(fine_delta[m], crse_delta);
           }
 
@@ -82,33 +90,44 @@ namespace pfasst
             fine_delta[0]->zero();
           }
 
-          fine->get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
+          fine.get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
 
-          for (size_t m = m0; m < nfine; m++) { fine->evaluate(m); }
+          for (size_t m = m0; m < nfine; m++) { fine.evaluate(m); }
         }
 
-        virtual void restrict(shared_ptr<ISweeper<time>> dst,
-                              shared_ptr<const ISweeper<time>> src,
-                              bool restrict_initial)
+
+        virtual void interpolate(shared_ptr<Encapsulation<time>> dst,
+                                 shared_ptr<const Encapsulation<time>> src)
         {
-          shared_ptr<EncapSweeper<time>> crse = dynamic_pointer_cast<EncapSweeper<time>>(dst);
-          assert(crse);
-          shared_ptr<const EncapSweeper<time>> fine = dynamic_pointer_cast<const EncapSweeper<time>>(src);
-          assert(fine);
+          UNUSED(dst); UNUSED(src);
+          throw NotImplementedYet("mlsdc/pfasst");
+        }
+        //! @}
 
-          this->restrict(crse, fine, restrict_initial);
+
+        //! @{
+        virtual void restrict_initial(shared_ptr<ISweeper<time>> dst,
+                                      shared_ptr<const ISweeper<time>> src) override
+        {
+          auto& crse = as_encap_sweeper(dst);
+          auto& fine = as_encap_sweeper(src);
+          this->restrict(crse.get_state(0), fine.get_state(0));
         }
 
-        virtual void restrict(shared_ptr<EncapSweeper<time>> crse,
-                              shared_ptr<const EncapSweeper<time>> fine,
-                              bool restrict_initial)
+
+        virtual void restrict(shared_ptr<ISweeper<time>> dst,
+                              shared_ptr<const ISweeper<time>> src,
+                              bool restrict_initial) override
         {
-          auto dnodes = crse->get_nodes();
-          auto snodes = fine->get_nodes();
+          auto& crse = as_encap_sweeper(dst);
+          auto& fine = as_encap_sweeper(src);
 
-          size_t ncrse = crse->get_nodes().size();
+          auto dnodes = crse.get_nodes();
+          auto snodes = fine.get_nodes();
+
+          size_t ncrse = crse.get_nodes().size();
           assert(ncrse > 1);
-          size_t nfine = fine->get_nodes().size();
+          size_t nfine = fine.get_nodes().size();
 
           int trat = (int(nfine) - 1) / (int(ncrse) - 1);
 
@@ -117,31 +136,32 @@ namespace pfasst
             if (dnodes[m] != snodes[m * trat]) {
               throw NotImplementedYet("coarse nodes must be nested");
             }
-            this->restrict(crse->get_state(m), fine->get_state(m * trat));
+            this->restrict(crse.get_state(m), fine.get_state(m * trat));
           }
 
-          for (size_t m = m0; m < ncrse; m++) { crse->evaluate(m); }
+          for (size_t m = m0; m < ncrse; m++) { crse.evaluate(m); }
         }
 
-        virtual void fas(time dt, shared_ptr<ISweeper<time>> dst,
-                                  shared_ptr<const ISweeper<time>> src)
-        {
-          shared_ptr<EncapSweeper<time>> crse = dynamic_pointer_cast<EncapSweeper<time>>(dst);
-          assert(crse);
-          shared_ptr<const EncapSweeper<time>> fine = dynamic_pointer_cast<const EncapSweeper<time>>(src);
-          assert(fine);
 
-          this->fas(dt, crse, fine);
+        virtual void restrict(shared_ptr<Encapsulation<time>> dst,
+                              shared_ptr<const Encapsulation<time>> src)
+        {
+          UNUSED(dst); UNUSED(src);
+          throw NotImplementedYet("mlsdc/pfasst");
         }
+        //! @}
 
-        virtual void fas(time dt, shared_ptr<EncapSweeper<time>> crse,
-                         shared_ptr<const EncapSweeper<time>> fine)
+        virtual void fas(time dt, shared_ptr<ISweeper<time>> dst,
+                         shared_ptr<const ISweeper<time>> src) override
         {
-          size_t ncrse = crse->get_nodes().size(); assert(ncrse >= 1);
-          size_t nfine = fine->get_nodes().size(); assert(nfine >= 1);
+          auto& crse = pfasst::encap::as_encap_sweeper(dst);
+          auto& fine = pfasst::encap::as_encap_sweeper(src);
+
+          size_t ncrse = crse.get_nodes().size(); assert(ncrse >= 1);
+          size_t nfine = fine.get_nodes().size(); assert(nfine >= 1);
 
-          auto crse_factory = crse->get_factory();
-          auto fine_factory = fine->get_factory();
+          auto crse_factory = crse.get_factory();
+          auto fine_factory = fine.get_factory();
 
           EncapVecT crse_z2n(ncrse - 1), fine_z2n(nfine - 1), rstr_z2n(ncrse - 1);
 
@@ -150,13 +170,13 @@ namespace pfasst
           for (size_t m = 0; m < nfine - 1; m++) { fine_z2n[m] = fine_factory->create(solution); }
 
           // compute '0 to node' integral on the coarse level
-          crse->integrate(dt, crse_z2n);
+          crse.integrate(dt, crse_z2n);
           for (size_t m = 1; m < ncrse - 1; m++) {
             crse_z2n[m]->saxpy(1.0, crse_z2n[m - 1]);
           }
 
           // compute '0 to node' integral on the fine level
-          fine->integrate(dt, fine_z2n);
+          fine.integrate(dt, fine_z2n);
           for (size_t m = 1; m < nfine - 1; m++) {
             fine_z2n[m]->saxpy(1.0, fine_z2n[m - 1]);
           }
@@ -169,7 +189,7 @@ namespace pfasst
 
           // compute 'node to node' tau correction
           EncapVecT tau(ncrse - 1), rstr_and_crse(2 * (ncrse - 1));
-          for (size_t m = 0; m < ncrse - 1; m++) { tau[m] = crse->get_tau(m); }
+          for (size_t m = 0; m < ncrse - 1; m++) { tau[m] = crse.get_tau(m); }
           for (size_t m = 0; m < ncrse - 1; m++) { rstr_and_crse[m] = rstr_z2n[m]; }
           for (size_t m = 0; m < ncrse - 1; m++) { rstr_and_crse[ncrse - 1 + m] = crse_z2n[m]; }
 
@@ -191,22 +211,9 @@ namespace pfasst
           tau[0]->mat_apply(tau, 1.0, fmat, rstr_and_crse, true);
         }
 
-        // required for interp/restrict helpers
-        virtual void interpolate(shared_ptr<Encapsulation<time>> /*dst*/,
-                                 shared_ptr<const Encapsulation<time>> /*src*/)
-        {
-          throw NotImplementedYet("mlsdc/pfasst");
-        }
-
-        virtual void restrict(shared_ptr<Encapsulation<time>> /*dst*/,
-                              shared_ptr<const Encapsulation<time>> /*src*/)
-        {
-          throw NotImplementedYet("mlsdc/pfasst");
-        }
-
     };
 
-  }
-}
+  }  // ::pfasst::encap
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/encap/vector.hpp b/include/pfasst/encap/vector.hpp
index 41b6ee2d..e2739c98 100644
--- a/include/pfasst/encap/vector.hpp
+++ b/include/pfasst/encap/vector.hpp
@@ -5,6 +5,7 @@
 #ifndef _PFASST_VECTOR_HPP_
 #define _PFASST_VECTOR_HPP_
 
+#include <cstdlib>
 #include <algorithm>
 #include <memory>
 #include <vector>
@@ -32,8 +33,6 @@ namespace pfasst
     {
       public:
         //! @{
-        /**
-         */
         VectorEncapsulation(int size)
           : vector<scalar>(size)
         {
@@ -74,14 +73,15 @@ namespace pfasst
         //! @}
 
         //! @{
-        void zero()
+        void zero() override
         {
           this->assign(this->size(), scalar(0.0));
         }
 
-        void copy(shared_ptr<const Encapsulation<time>> x)
+        void copy(shared_ptr<const Encapsulation<time>> x) override
         {
-          shared_ptr<const VectorEncapsulation<scalar, time>> x_cast = dynamic_pointer_cast<const VectorEncapsulation<scalar, time>>(x);
+          shared_ptr<const VectorEncapsulation<scalar, time>> x_cast = \
+            dynamic_pointer_cast<const VectorEncapsulation<scalar, time>>(x);
           assert(x_cast);
           this->copy(x_cast);
         }
@@ -93,9 +93,10 @@ namespace pfasst
         //! @}
 
         //! @{
-        void saxpy(time a, shared_ptr<const Encapsulation<time>> x)
+        void saxpy(time a, shared_ptr<const Encapsulation<time>> x) override
         {
-          shared_ptr<const VectorEncapsulation<scalar, time>> x_cast = dynamic_pointer_cast<const VectorEncapsulation<scalar, time>>(x);
+          shared_ptr<const VectorEncapsulation<scalar, time>> x_cast = \
+            dynamic_pointer_cast<const VectorEncapsulation<scalar, time>>(x);
           assert(x_cast);
 
           this->saxpy(a, x_cast);
@@ -113,7 +114,7 @@ namespace pfasst
          *     `dynamic_cast` into pfasst::encap::VectorEncapsulation std::abort is called.
          */
         void mat_apply(vector<shared_ptr<Encapsulation<time>>> dst, time a, matrix<time> mat,
-                       vector<shared_ptr<Encapsulation<time>>> src, bool zero = true)
+                       vector<shared_ptr<Encapsulation<time>>> src, bool zero = true) override
         {
           size_t ndst = dst.size();
           size_t nsrc = src.size();
@@ -131,7 +132,8 @@ namespace pfasst
           dst_cast[0]->mat_apply(dst_cast, a, mat, src_cast, zero);
         }
 
-        void mat_apply(vector<shared_ptr<VectorEncapsulation<scalar, time>>> dst, time a, matrix<time> mat,
+        void mat_apply(vector<shared_ptr<VectorEncapsulation<scalar, time>>> dst,
+                       time a, matrix<time> mat,
                        vector<shared_ptr<VectorEncapsulation<scalar, time>>> src, bool zero = true)
         {
           size_t ndst = dst.size();
@@ -171,7 +173,8 @@ namespace pfasst
      *     precision of the time points; defaults to pfasst::time_precision
      */
     template<typename scalar, typename time = time_precision>
-    class VectorFactory : public EncapFactory<time>
+    class VectorFactory
+      : public EncapFactory<time>
     {
         int size;
       public:
@@ -201,7 +204,7 @@ namespace pfasst
       return *y.get();
     }
 
-  }
-}
+  }  // ::pfasst::encap
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/globals.hpp b/include/pfasst/globals.hpp
new file mode 100644
index 00000000..80d2d3f0
--- /dev/null
+++ b/include/pfasst/globals.hpp
@@ -0,0 +1,30 @@
+#ifndef _GLOBALS__HPP_
+#define _GLOBALS__HPP_
+
+/**
+ * denoting unused function parameters for omitting compiler warnings
+ *
+ * To denote an unused function parameter just use this makro on it in the function body:
+ *
+ * ~~~{.cpp}
+ * void func(auto p) {
+ *   UNUSED(p);
+ * }
+ * ~~~
+ *
+ * which renders to
+ *
+ * ~~~{.cpp}
+ * void func(auto p) {
+ *   (void)(p);
+ * }
+ * ~~~
+ *
+ * which is the standard and compiler independet way of omitting warnings on unused parameters while
+ * still being able to fully document the parameter with Doxygen.
+ *
+ * @param[in] expr parameter to be denoted as unused
+ */
+#define UNUSED(expr) (void)(expr)
+
+#endif  // _GLOBALS__HPP_
diff --git a/include/pfasst/interfaces.hpp b/include/pfasst/interfaces.hpp
index fa0d1778..92234105 100644
--- a/include/pfasst/interfaces.hpp
+++ b/include/pfasst/interfaces.hpp
@@ -12,6 +12,8 @@
 #include <memory>
 #include <string>
 
+#include "globals.hpp"
+
 using namespace std;
 
 namespace pfasst
@@ -19,6 +21,7 @@ namespace pfasst
 
   using time_precision = double;
 
+  // forward declare for ISweeper
   template<typename time>
   class Controller;
 
@@ -28,33 +31,50 @@ namespace pfasst
    * Used by PFASST to mark methods that are required for a particular algorithm (SDC/MLSDC/PFASST)
    * that may not be necessary for all others.
    */
-  class NotImplementedYet : public exception
+  class NotImplementedYet
+    : public exception
   {
       string msg;
     public:
-      NotImplementedYet(string msg) : msg(msg) { }
+      NotImplementedYet(string msg)
+        : msg(msg)
+      {}
+
       const char* what() const throw()
       {
         return (string("Not implemented/supported yet, required for: ") + this->msg).c_str();
       }
   };
 
+
   /**
    * value exception.
    *
    * Thrown when a PFASST routine is passed an invalid value.
    */
-  class ValueError : public exception
+  class ValueError
+    : public exception
   {
       string msg;
     public:
-      ValueError(string msg) : msg(msg) { }
+      ValueError(string msg)
+        : msg(msg)
+      {}
+
       const char* what() const throw()
       {
         return (string("ValueError: ") + this->msg).c_str();
       }
   };
 
+  class ICommunicator
+  {
+    public:
+      virtual ~ICommunicator() { }
+      virtual int size() = 0;
+      virtual int rank() = 0;
+  };
+
   /**
    * abstract SDC sweeper.
    * @tparam time time precision
@@ -63,32 +83,42 @@ namespace pfasst
   template<typename time = time_precision>
   class ISweeper
   {
+    protected:
       Controller<time>* controller;
 
     public:
-      virtual ~ISweeper() { }
+      //! @{
+      virtual ~ISweeper()
+      {}
+      //! @}
 
+      //! @{
       /**
        * set the sweepers controller.
        */
       void set_controller(Controller<time>* ctrl)
       {
-	controller = ctrl;
+        this->controller = ctrl;
       }
 
       Controller<time>* get_controller()
       {
-	assert(controller);
-	return controller;
+        assert(this->controller);
+        return this->controller;
       }
+      //! @}
 
+      //! @{
       /**
        * setup (allocate etc) the sweeper.
        * @param[in] coarse
        *     `true` if this sweeper exists on a coarsened MLSDC or PFASST level.
        *     This implies that space for an FAS correction and "saved" solutions are necessary.
        */
-      virtual void setup(bool coarse = false) { (void) coarse; }
+      virtual void setup(bool coarse = false)
+      {
+        UNUSED(coarse);
+      }
 
       /**
        * perform a predictor sweep.
@@ -104,15 +134,16 @@ namespace pfasst
       virtual void predict(bool initial) = 0;
 
       /**
-       * perform one SDC sweep/iteration.
-       * Compute a correction and update solution values.
-       * Note that this function can assume that valid function values exist from a previous
-       * pfasst::ISweeper::sweep() or pfasst::ISweeper::predict().
+       * Perform one SDC sweep/iteration.
+       *
+       * Compute a correction and update solution values.  Note that this function can assume that
+       * valid function values exist from a previous pfasst::ISweeper::sweep() or
+       * pfasst::ISweeper::predict().
        */
       virtual void sweep() = 0;
 
       /**
-       * advance from one time step to the next.
+       * Advance from one time step to the next.
        *
        * Essentially this means copying the solution and function values from the last node to the
        * first node.
@@ -120,12 +151,49 @@ namespace pfasst
       virtual void advance() = 0;
 
       /**
-       * save solutions (and/or function values) at all nodes.
+       * Save states (and/or function values) at all nodes.
        *
        * This is typically done in MLSDC/PFASST immediately after a call to restrict.
        * The saved states are used to compute deltas during interpolation.
+       *
+       * @note This method must be implemented in derived sweepers.
        */
-      virtual void save() { NotImplementedYet("mlsdc/pfasst"); }
+      virtual void save(bool initial_only=false)
+      {
+        UNUSED(initial_only);
+        throw NotImplementedYet("mlsdc/pfasst");
+      }
+
+      virtual void spread()
+      {
+        throw NotImplementedYet("pfasst");
+      }
+      //! @}
+
+      //! @{
+      virtual void post(ICommunicator* comm, int tag)
+      {
+        UNUSED(comm); UNUSED(tag);
+      };
+
+      virtual void send(ICommunicator* comm, int tag, bool blocking)
+      {
+        UNUSED(comm); UNUSED(tag); UNUSED(blocking);
+        throw NotImplementedYet("pfasst");
+      }
+
+      virtual void recv(ICommunicator* comm, int tag, bool blocking)
+      {
+        UNUSED(comm); UNUSED(tag); UNUSED(blocking);
+        throw NotImplementedYet("pfasst");
+      }
+
+      virtual void broadcast(ICommunicator* comm)
+      {
+        UNUSED(comm);
+        throw NotImplementedYet("pfasst");
+      }
+      //! @}
 
   };
 
@@ -138,23 +206,46 @@ namespace pfasst
   class ITransfer
   {
     public:
-      // XXX: pass level iterator to these routines as well
-      virtual ~ITransfer() { }
+      //! @{
+      virtual ~ITransfer()
+      {}
+      //! @}
 
+      //! @{
       /**
-       * interpolate, in time and space, from the coarse sweeper to the fine sweeper.
-       * @param[in] interp_delta_from_initial
-       *     `true` if the delta computed at each node should be relative to the initial condition.
+       * Interpolate initial condition from the coarse sweeper to the fine sweeper.
+       */
+      virtual void interpolate_initial(shared_ptr<ISweeper<time>> dst,
+                                       shared_ptr<const ISweeper<time>> src)
+      {
+        UNUSED(dst); UNUSED(src);
+        throw NotImplementedYet("pfasst");
+      }
+
+      /**
+       * Interpolate, in time and space, from the coarse sweeper to the fine sweeper.
        * @param[in] interp_initial
        *     `true` if a delta for the initial condtion should also be computed (PFASST).
        */
       virtual void interpolate(shared_ptr<ISweeper<time>> dst,
                                shared_ptr<const ISweeper<time>> src,
-                               bool interp_delta_from_initial = false,
                                bool interp_initial = false) = 0;
 
       /**
-       * restrict, in time and space, from the fine sweeper to the coarse sweeper.
+       * Restrict initial condition from the fine sweeper to the coarse sweeper.
+       * @param[in] restrict_initial
+       *     `true` if the initial condition should also be restricted.
+       */
+      virtual void restrict_initial(shared_ptr<ISweeper<time>> dst,
+                                    shared_ptr<const ISweeper<time>> src)
+      {
+        UNUSED(dst); UNUSED(src);
+        throw NotImplementedYet("pfasst");
+      }
+
+
+      /**
+       * Restrict, in time and space, from the fine sweeper to the coarse sweeper.
        * @param[in] restrict_initial
        *     `true` if the initial condition should also be restricted.
        */
@@ -162,25 +253,15 @@ namespace pfasst
                             shared_ptr<const ISweeper<time>> src,
                             bool restrict_initial = false) = 0;
 
+
       /**
-       * compute FAS correction between the coarse and fine sweepers.
+       * Compute FAS correction between the coarse and fine sweepers.
        */
       virtual void fas(time dt, shared_ptr<ISweeper<time>> dst,
                        shared_ptr<const ISweeper<time>> src) = 0;
-
-  };
-
-  /**
-   * abstract time communicator.
-   */
-  class ICommunicator
-  {
-    public:
-      virtual void post() { }
-      virtual void send() = 0;
-      virtual void recv() = 0;
+      //! @}
   };
 
-}
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/mlsdc.hpp b/include/pfasst/mlsdc.hpp
index af315e0d..d1109522 100644
--- a/include/pfasst/mlsdc.hpp
+++ b/include/pfasst/mlsdc.hpp
@@ -21,8 +21,10 @@ namespace pfasst
   class MLSDC
     : public Controller<time>
   {
+    protected:
       vector<size_t> nsweeps;
-      bool predict, initial;
+      bool predict;
+      bool initial;
 
       typedef typename pfasst::Controller<time>::LevelIter LevelIter;
 
@@ -40,16 +42,21 @@ namespace pfasst
       }
 
     public:
-      void setup()
+      void setup() override
       {
         nsweeps.resize(this->nlevels());
         fill(nsweeps.begin(), nsweeps.end(), 1);
         for (auto leviter = this->coarsest(); leviter <= this->finest(); ++leviter) {
-	  leviter.current()->set_controller(this);
+          leviter.current()->set_controller(this);
           leviter.current()->setup(leviter != this->finest());
         }
       }
 
+      void set_nsweeps(vector<size_t> nsweeps)
+      {
+        this->nsweeps = nsweeps;
+      }
+
       /**
        * evolve ODE using MLSDC.
        *
@@ -59,11 +66,13 @@ namespace pfasst
       void run()
       {
         for (; this->get_time() < this->get_end_time(); this->advance_time()) {
-          predict = true;   // use predictor for first fine sweep of each step
-          initial = this->get_step() == 0; // only evaluate node 0 functions on first step
+          predict = true;  // use predictor for first fine sweep of each step
+          initial = this->get_step() == 0;  // only evaluate node 0 functions on first step
 
           // iterate by performing v-cycles
-	  for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations(); this->advance_iteration()) {
+          for (this->set_iteration(0);
+               this->get_iteration() < this->get_max_iterations();
+               this->advance_iteration()) {
             cycle_v(this->finest());
           }
 
@@ -140,6 +149,6 @@ namespace pfasst
 
   };
 
-}
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/mpi_communicator.hpp b/include/pfasst/mpi_communicator.hpp
new file mode 100644
index 00000000..d250c377
--- /dev/null
+++ b/include/pfasst/mpi_communicator.hpp
@@ -0,0 +1,74 @@
+/*
+ * Interfaces for SDC/MLSDC/PFASST algorithms.
+ */
+
+#ifndef _PFASST_MPI_COMMUNICATOR_HPP_
+#define _PFASST_MPI_COMMUNICATOR_HPP_
+
+#include <exception>
+
+#include <mpi.h>
+
+#include "interfaces.hpp"
+
+using namespace std;
+
+namespace pfasst
+{
+  namespace mpi
+  {
+
+    class MPIError
+      : public exception
+    {
+      public:
+        const char* what() const throw()
+        {
+          return "mpi error";
+        }
+    };
+
+
+    class MPICommunicator
+      : public ICommunicator
+    {
+        //! @{
+        int _rank;
+        int _size;
+        //! @}
+
+      public:
+        //! @{
+        MPI_Comm comm;
+        //! @}
+
+        //! @{
+        MPICommunicator()
+        {}
+
+        MPICommunicator(MPI_Comm comm)
+        {
+          set_comm(comm);
+        }
+
+        virtual ~MPICommunicator()
+        {}
+        //! @}
+
+        //! @{
+        void set_comm(MPI_Comm comm)
+        {
+          this->comm = comm;
+          MPI_Comm_size(this->comm, &(this->_size));
+          MPI_Comm_rank(this->comm, &(this->_rank));
+        }
+
+        int size() { return this->_size; }
+        int rank() { return this->_rank; }
+        //! @
+    };
+
+  }  // ::pfasst::mpi
+}  // ::pfasst
+
+#endif
diff --git a/include/pfasst/pfasst.hpp b/include/pfasst/pfasst.hpp
index 00dd06af..f5c0a91a 100644
--- a/include/pfasst/pfasst.hpp
+++ b/include/pfasst/pfasst.hpp
@@ -5,30 +5,224 @@
 #ifndef _PFASST_PFASST_HPP_
 #define _PFASST_PFASST_HPP_
 
-#include "controller.hpp"
+#include "mlsdc.hpp"
 
 using namespace std;
 
 namespace pfasst
 {
 
-  template<typename timeT>
-  class PFASST 
-    : public Controller<timeT>
+  template<typename time = pfasst::time_precision>
+  class PFASST
+    : public MLSDC<time>
   {
+      ICommunicator* comm;
 
-      void run(timeT dt, size_t nsteps)
+      typedef typename pfasst::Controller<time>::LevelIter LevelIter;
+
+      bool predict; //<! whether to use a 'predict' sweep
+      bool initial; //<! whether we're sweeping from a new initial condition
+
+      void perform_sweeps(size_t level)
+      {
+        auto sweeper = this->get_level(level);
+        for (size_t s = 0; s < this->nsweeps[level]; s++) {
+          if (predict) {
+            sweeper->predict(initial & predict);
+            predict = false;
+          } else {
+            sweeper->sweep();
+          }
+        }
+      }
+
+    public:
+      void set_comm(ICommunicator* comm)
+      {
+        this->comm = comm;
+      }
+
+      /**
+       * Predictor: restrict initial down, preform coarse sweeps, return to finest.
+       */
+      void predictor()
+      {
+        this->get_finest()->spread();
+
+        // restrict fine initial condition
+        for (auto l = this->finest() - 1; l >= this->coarsest(); --l) {
+          auto crse = l.current();
+          auto fine = l.fine();
+          auto trns = l.transfer();
+          trns->restrict_initial(crse, fine);
+          crse->spread();
+          crse->save();
+        }
+
+        // perform sweeps on the coarse level based on rank
+        predict = true;
+        auto crse = this->coarsest().current();
+        for (int nstep = 0; nstep < comm->rank() + 1; nstep++) {
+          // XXX: set iteration and step?
+          perform_sweeps(0);
+          if (nstep < comm->rank()) {
+            crse->advance();
+          }
+        }
+
+        // return to finest level, sweeping as we go
+        for (auto l = this->coarsest() + 1; l <= this->finest(); ++l) {
+          auto crse = l.coarse();
+          auto fine = l.current();
+          auto trns = l.transfer();
+
+          trns->interpolate(fine, crse, true);
+          if (l < this->finest()) {
+            perform_sweeps(l.level);
+          }
+        }
+      }
+
+      void broadcast()
+      {
+        this->get_finest()->broadcast(comm);
+        initial = true;
+      }
+
+      int tag (int level)
+      {
+        return level * 10000 + this->get_iteration() + 10;
+      }
+
+      int tag(LevelIter l)
+      {
+        return tag(l.level);
+      }
+
+      void post()
+      {
+        for (auto l = this->coarsest() + 1; l <= this->finest(); ++l) {
+          l.current()->post(comm, tag(l));
+        }
+      }
+
+
+      /**
+       * Evolve ODE using PFASST.
+       *
+       * This assumes that the user has set initial conditions on the
+       * finest level.
+       *
+       * Currently uses "block mode" PFASST with the standard predictor.
+       */
+      void run()
+      {
+        int nblocks = int(this->get_end_time() / this->get_time_step()) / comm->size();
+
+        if (nblocks == 0) {
+          throw ValueError("invalid duration: there are more time processors than time steps");
+        }
+
+        initial = true;
+        for (int nblock = 0; nblock < nblocks; nblock++) {
+          this->set_step(nblock * comm->size() + comm->rank());
+
+          predictor();
+
+          for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations();
+               this->advance_iteration()) {
+            post();
+            cycle_v(this->finest());
+          }
+
+          if (nblock < nblocks - 1) {
+            broadcast();
+          }
+        }
+      }
+
+      /**
+       * Cycle down: sweep on current (fine), restrict to coarse.
+       */
+      LevelIter cycle_down(LevelIter l)
       {
+        auto fine = l.current();
+        auto crse = l.coarse();
+        auto trns = l.transfer();
+
+        perform_sweeps(l.level);
+
+        if (l == this->finest()) {
+          // note: convergence tests belong here
+        }
 
+        fine->send(comm, tag(l), false);
+
+        auto dt = this->get_time_step();
+        trns->restrict(crse, fine, true);
+        trns->fas(dt, crse, fine);
+        crse->save();
+
+        return l - 1;
       }
 
-      // leveliter cycle_down(leveliter levels, double t, double dt);
-      // leveliter cycle_up(leveliter levels, double t, double dt);
-      // leveliter cycle_bottom(leveliter levels, double t, double dt);
-      // leveliter cycle_top(leveliter levels, double t, double dt);
-      // leveliter cycle_v(leveliter levels, double t, double dt);
+      /**
+       * Cycle up: interpolate coarse correction to fine, sweep on
+       * current (fine).
+       *
+       * Note that if the fine level corresponds to the finest MLSDC
+       * level, we don't perform a sweep.  In this case the only
+       * operation that is performed here is interpolation.
+       */
+      LevelIter cycle_up(LevelIter l)
+      {
+        auto fine = l.current();
+        auto crse = l.coarse();
+        auto trns = l.transfer();
+
+        trns->interpolate(fine, crse, true);
+
+        fine->recv(comm, tag(l), false);
+        trns->interpolate_initial(fine, crse);
+
+        if (l < this->finest()) {
+          perform_sweeps(l.level);
+        }
+
+        return l + 1;
+      }
+
+      /**
+       * Cycle bottom: sweep on the current (coarsest) level.
+       */
+      LevelIter cycle_bottom(LevelIter l)
+      {
+        auto crse = l.current();
+
+        crse->recv(comm, tag(l), true);
+        this->perform_sweeps(l.level);
+        crse->send(comm, tag(l), true);
+        return l + 1;
+      }
+
+      /**
+       * Perform an MLSDC V-cycle.
+       */
+      LevelIter cycle_v(LevelIter l)
+      {
+        // this v-cycle is a bit different than in mlsdc
+
+        if (l.level == 0) {
+          l = cycle_bottom(l);
+        } else {
+          l = cycle_down(l);
+          l = cycle_v(l);
+          l = cycle_up(l);
+        }
+        return l;
+      }
   };
 
-}
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/quadrature.hpp b/include/pfasst/quadrature.hpp
index 301ca72b..8ad78378 100644
--- a/include/pfasst/quadrature.hpp
+++ b/include/pfasst/quadrature.hpp
@@ -29,7 +29,8 @@ namespace pfasst
       vector<CoeffT> c;
 
     public:
-      Polynomial(size_t n) : c(n)
+      Polynomial(size_t n)
+        : c(n)
       {
         fill(c.begin(), c.end(), 0.0);
       }
@@ -160,20 +161,28 @@ namespace pfasst
       }
   };
 
+
+  enum class QuadratureType {
+      GaussLegendre
+    , GaussLobatto
+    , GaussRadau
+    , ClenshawCurtis
+    , Uniform
+  };
+
+
   template<typename node = time_precision>
-  vector<node> compute_nodes(size_t nnodes, string qtype)
+  vector<node> compute_nodes(size_t nnodes, QuadratureType qtype)
   {
     vector<node> nodes(nnodes);
 
-    if (qtype == "gauss-legendre") {
-
+    if (qtype == QuadratureType::GaussLegendre) {
       auto roots = Polynomial<node>::legendre(nnodes).roots();
       for (size_t j = 0; j < nnodes; j++) {
         nodes[j] = 0.5 * (1.0 + roots[j]);
       }
 
-    } else if (qtype == "gauss-lobatto") {
-
+    } else if (qtype == QuadratureType::GaussLobatto) {
       auto roots = Polynomial<node>::legendre(nnodes - 1).differentiate().roots();
       assert(nnodes >= 2);
       for (size_t j = 0; j < nnodes - 2; j++) {
@@ -182,8 +191,7 @@ namespace pfasst
       nodes.front() = 0.0;
       nodes.back() = 1.0;
 
-    } else if (qtype == "gauss-radau") {
-
+    } else if (qtype == QuadratureType::GaussRadau) {
       auto l   = Polynomial<node>::legendre(nnodes);
       auto lm1 = Polynomial<node>::legendre(nnodes - 1);
       for (size_t i = 0; i < nnodes; i++) {
@@ -195,45 +203,63 @@ namespace pfasst
       }
       nodes.back() = 1.0;
 
-    } else if (qtype == "clenshaw-curtis") {
-
+    } else if (qtype == QuadratureType::ClenshawCurtis) {
       for (size_t j = 0; j < nnodes; j++) {
-        nodes[j] = 0.5 * (1.0 - cos(j * PI / (nnodes-1)));
+        nodes[j] = 0.5 * (1.0 - cos(j * PI / (nnodes - 1)));
       }
 
-    } else if (qtype == "uniform") {
-
+    } else if (qtype == QuadratureType::Uniform) {
       for (size_t j = 0; j < nnodes; j++) {
-        nodes[j] = node(j) / (nnodes-1);
+        nodes[j] = node(j) / (nnodes - 1);
       }
 
     } else {
-
       throw ValueError("invalid node type passed to compute_nodes.");
-
     }
 
     return nodes;
   }
 
+  template<typename node>
+  auto augment_nodes(vector<node> const orig) -> pair<vector<node>, vector<bool>> {
+    vector<node> nodes = orig;
+
+    bool left = nodes.front() == node(0.0);
+    bool right = nodes.back() == node(1.0);
+
+    if (!left)  { nodes.insert(nodes.begin(), node(0.0)); }
+    if (!right) { nodes.insert(nodes.end(),   node(1.0)); }
+
+    vector<bool> is_proper(nodes.size(), true);
+    is_proper.front() = left;
+    is_proper.back() = right;
+
+    return pair<vector<node>, vector<bool>>(nodes, is_proper);
+  }
+
+//  enum class QuadratureMatrix { S, Q, QQ }; // returning QQ might be cool for 2nd-order stuff
+  enum class QuadratureMatrix { S, Q };
+  
   template<typename node = time_precision>
-  matrix<node> compute_quadrature(vector<node> dst, vector<node> src, char type)
+  matrix<node> compute_quadrature(vector<node> dst, vector<node> src, vector<bool> is_proper,
+                                  QuadratureMatrix type)
   {
     const size_t ndst = dst.size();
     const size_t nsrc = src.size();
 
     assert(ndst >= 1);
-    matrix<node> mat(ndst - 1, nsrc);
+    matrix<node> mat(ndst - 1, nsrc, node(0.0));
 
     Polynomial<node> p(nsrc + 1), p1(nsrc + 1);
 
     for (size_t i = 0; i < nsrc; i++) {
+      if (!is_proper[i]) { continue; }
 
       // construct interpolating polynomial coefficients
       p[0] = 1.0;
       for (size_t j = 1; j < nsrc + 1; j++) { p[j] = 0.0; }
       for (size_t m = 0; m < nsrc; m++) {
-        if (m == i) { continue; }
+        if ((!is_proper[m]) || (m == i)) { continue; }
 
         // p_{m+1}(x) = (x - x_j) * p_m(x)
         p1[0] = 0.0;
@@ -247,12 +273,14 @@ namespace pfasst
       auto P = p.integrate();
       for (size_t j = 1; j < ndst; j++) {
         node q = 0.0;
-        if (type == 's') {
+        if (type == QuadratureMatrix::S) {
           q = P.evaluate(dst[j]) - P.evaluate(dst[j - 1]);
-        } else {
+        } else if (type == QuadratureMatrix::Q) {
           q = P.evaluate(dst[j]) - P.evaluate(0.0);
+        } else {
+          throw ValueError("Further matrix types are not implemented yet");
         }
-
+          
         mat(j - 1, i) = q / den;
       }
     }
@@ -282,14 +310,14 @@ namespace pfasst
         if (abs(num) > 1e-32) {
           mat(i, j) = num / den;
         } else {
-	  mat(i, j) = 0.0;
-	}
+          mat(i, j) = 0.0;
+        }
       }
     }
 
     return mat;
   }
 
-}
+}  // ::pfasst
 
 #endif
diff --git a/include/pfasst/sdc.hpp b/include/pfasst/sdc.hpp
index 2eb2be16..d40d87dd 100644
--- a/include/pfasst/sdc.hpp
+++ b/include/pfasst/sdc.hpp
@@ -18,27 +18,27 @@ namespace pfasst
     : public Controller<time>
   {
     public:
-
       void run()
       {
         auto sweeper = this->get_level(0);
 
         for (; this->get_time() < this->get_end_time(); this->advance_time()) {
           bool initial = this->get_step() == 0;
-	  for (this->set_iteration(0); this->get_iteration() < this->get_max_iterations(); this->advance_iteration()) {
-	    bool predict = this->get_iteration() == 0;
-	    if (predict) {
-	      sweeper->predict(initial);
-	    } else {
-	      sweeper->sweep();
-	    }
+          for (this->set_iteration(0);
+               this->get_iteration() < this->get_max_iterations();
+               this->advance_iteration()) {
+            bool predict = this->get_iteration() == 0;
+            if (predict) {
+              sweeper->predict(initial);
+            } else {
+              sweeper->sweep();
+            }
           }
           sweeper->advance();
         }
       }
-
   };
 
-}
+}  // ::pfasst
 
 #endif
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b4c29d0d..56471aaf 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,28 +1,33 @@
 # Building and Running Tests
 include_directories(
     ${3rdparty_INCLUDES}
-    ${fftw3_INCLUDES}
     ${pfasst_INCLUDES}
 )
 
 set(TESTS
-    test-quadrature
-    test-advection-diffusion
+    test_quadrature
 )
 
 foreach(test ${TESTS})
     message(STATUS "  ${test}")
     add_executable(${test} ${test}.cpp)
-    add_dependencies(${test} googlemock fftw3)
+    add_dependencies(${test} googlemock)
     target_link_libraries(${test}
         ${3rdparty_DEPENDEND_LIBS}
-        ${fftw3_LIBS}
         ${pfasst_DEPENDED_LIBS}
     )
     set_target_properties(${test}
         PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests
     )
+    if(pfasst_WITH_GCC_PROF AND ${CMAKE_CXX_COMPILER_ID} MATCHES GNU)
+        set_target_properties(${test}
+            PROPERTIES COMPILE_FLAGS "${CMAKE_CXX_FLAGS} -ftest-coverage -fprofile-arcs"
+                       LINK_FLAGS "-fprofile-arcs"
+        )
+    endif()
     add_test(NAME ${test}
         COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/${test} --gtest_output=xml:${test}_out.xml
     )
 endforeach(test)
+
+add_subdirectory(examples)
diff --git a/tests/examples/CMakeLists.txt b/tests/examples/CMakeLists.txt
new file mode 100644
index 00000000..ae247470
--- /dev/null
+++ b/tests/examples/CMakeLists.txt
@@ -0,0 +1,3 @@
+
+add_subdirectory(advection_diffusion)
+add_subdirectory(scalar)
diff --git a/tests/examples/advection_diffusion/CMakeLists.txt b/tests/examples/advection_diffusion/CMakeLists.txt
new file mode 100644
index 00000000..0ccaa798
--- /dev/null
+++ b/tests/examples/advection_diffusion/CMakeLists.txt
@@ -0,0 +1,79 @@
+# Building and Running Tests
+include_directories(
+    ${3rdparty_INCLUDES}
+    ${FFTW_INCLUDE_PATH}
+    ${pfasst_INCLUDES}
+)
+
+set(TESTS
+    test_advection_diffusion
+)
+
+if(${pfasst_WITH_MPI})
+    set(MPI_TESTS
+        test_mpi_advection_diffusion
+    )
+endif()
+
+foreach(test ${TESTS})
+    message(STATUS "  ${test}")
+    add_executable(${test} ${test}.cpp)
+    add_dependencies(${test} googlemock)
+    if(NOT FFTW_FOUND)
+      add_dependencies(${test} googlemock fftw3)
+    endif()
+    target_link_libraries(${test}
+        ${3rdparty_DEPENDEND_LIBS}
+        ${FFTW_LIBRARIES}
+        ${pfasst_DEPENDED_LIBS}
+    )
+    set_target_properties(${test}
+        PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/advection_diffusion
+    )
+    if(pfasst_WITH_GCC_PROF AND ${CMAKE_CXX_COMPILER_ID} MATCHES GNU)
+        set_target_properties(${test}
+            PROPERTIES COMPILE_FLAGS "${CMAKE_CXX_FLAGS} -ftest-coverage -fprofile-arcs"
+                       LINK_FLAGS "-fprofile-arcs"
+        )
+    endif()
+    add_test(NAME ${test}
+        COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/advection_diffusion/${test} --gtest_output=xml:${test}_out.xml
+    )
+endforeach(test)
+
+if(${pfasst_WITH_MPI})
+    include_directories(${MPI_CXX_INCLUDE_PATH})
+    foreach(test ${MPI_TESTS})
+        message(STATUS "  ${test}")
+        add_executable(${test} ${test}.cpp)
+        add_dependencies(${test} googlemock)
+        if(NOT FFTW_FOUND)
+          add_dependencies(${test} googlemock fftw3)
+        endif()
+        if(MPI_COMPILE_FLAGS)
+            if(pfasst_WITH_GCC_PROF AND ${CMAKE_CXX_COMPILER_ID} MATCHES GNU)
+                set_target_properties(${test}
+                    PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS} -ftest-coverage -fprofile-arcs"
+                               LINK_FLAGS "-fprofile-arcs"
+                )
+            else()
+                set_target_properties(${test} PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
+            endif()
+        endif()
+        if(MPI_LINK_FLAGS)
+            set_target_properties(${test} PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+        endif()
+        target_link_libraries(${test}
+            ${3rdparty_DEPENDEND_LIBS}
+            ${FFTW_LIBRARIES}
+            ${pfasst_DEPENDED_LIBS}
+            ${MPI_CXX_LIBRARIES}
+        )
+        set_target_properties(${test}
+            PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/advection_diffusion
+        )
+        add_test(NAME ${test}
+            COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/advection_diffusion/${test} ${MPIEXEC_POSTFLAGS}
+        )
+    endforeach(test)
+endif()
diff --git a/tests/test-advection-diffusion.cpp b/tests/examples/advection_diffusion/test_advection_diffusion.cpp
similarity index 100%
rename from tests/test-advection-diffusion.cpp
rename to tests/examples/advection_diffusion/test_advection_diffusion.cpp
diff --git a/tests/examples/advection_diffusion/test_mpi_advection_diffusion.cpp b/tests/examples/advection_diffusion/test_mpi_advection_diffusion.cpp
new file mode 100644
index 00000000..751b382d
--- /dev/null
+++ b/tests/examples/advection_diffusion/test_mpi_advection_diffusion.cpp
@@ -0,0 +1,52 @@
+/*
+ * Tests for the advection-diffusion examples.
+ */
+
+#include <algorithm>
+#include <iostream>
+#include <tuple>
+#include <vector>
+
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+using namespace ::testing;
+
+#include <mpi.h>
+
+#define PFASST_UNIT_TESTING
+#include "../examples/advection_diffusion/mpi_pfasst.cpp"
+#undef PFASST_UNIT_TESTING
+
+using namespace std;
+
+TEST(ErrorTest, MPIPFASST)
+{
+  typedef error_map::value_type vtype;
+
+  auto errors = run_mpi_pfasst();
+  auto get_step  = [](const vtype x) { return get<0>(get<0>(x)); };
+  auto get_iter  = [](const vtype x) { return get<1>(get<0>(x)); };
+  auto get_error = [](const vtype x) { return get<1>(x); };
+
+  auto max_iter = get_iter(*std::max_element(errors.begin(), errors.end(),
+             [get_iter](const vtype p1, const vtype p2) { return get_iter(p1) < get_iter(p2); }));
+
+  vector<double> ex = { 1.1168e-13, 4.8849e-13, 5.3268e-13, 2.3059e-12 };
+  for (auto& x: errors) {
+    if (get_iter(x) == max_iter) {
+      EXPECT_NEAR(get_error(x), ex[get_step(x)], *max_element(ex.cbegin(), ex.cend()));
+    }
+  }
+}
+
+int main(int argc, char** argv)
+{
+  testing::InitGoogleTest(&argc, argv);
+  MPI_Init(&argc, &argv);
+  int result = 1, max_result;  // GTest return value 1 (failure), 0 (success)
+  result = RUN_ALL_TESTS();
+  MPI_Allreduce(&result, &max_result, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Finalize();
+  return max_result;
+}
diff --git a/tests/examples/scalar/CMakeLists.txt b/tests/examples/scalar/CMakeLists.txt
new file mode 100644
index 00000000..550616d9
--- /dev/null
+++ b/tests/examples/scalar/CMakeLists.txt
@@ -0,0 +1,32 @@
+# Building and Running Tests
+include_directories(
+    ${3rdparty_INCLUDES}
+    ${pfasst_INCLUDES}
+)
+
+set(TESTS
+    test_scalar_conv
+    test_scalar_highprec
+)
+
+foreach(test ${TESTS})
+    message(STATUS "  ${test}")
+    add_executable(${test} ${test}.cpp)
+    add_dependencies(${test} googlemock)
+    target_link_libraries(${test}
+        ${3rdparty_DEPENDEND_LIBS}
+        ${pfasst_DEPENDED_LIBS}
+    )
+    set_target_properties(${test}
+        PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/scalar
+    )
+    if(pfasst_WITH_GCC_PROF AND ${CMAKE_CXX_COMPILER_ID} MATCHES GNU)
+        set_target_properties(${test}
+            PROPERTIES COMPILE_FLAGS "${CMAKE_CXX_FLAGS} -ftest-coverage -fprofile-arcs"
+                       LINK_FLAGS "-fprofile-arcs"
+        )
+    endif()
+    add_test(NAME ${test}
+        COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tests/examples/scalar/${test} --gtest_output=xml:${test}_out.xml
+    )
+endforeach(test)
diff --git a/tests/examples/scalar/test_scalar_conv.cpp b/tests/examples/scalar/test_scalar_conv.cpp
new file mode 100644
index 00000000..9392bf04
--- /dev/null
+++ b/tests/examples/scalar/test_scalar_conv.cpp
@@ -0,0 +1,203 @@
+/*
+ * Tests for the scalar example solving the test equation
+ */
+#include <cmath>
+
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+using namespace ::testing;
+
+#include <pfasst/quadrature.hpp>
+
+#define PFASST_UNIT_TESTING
+#include "../examples/scalar/scalar_sdc.cpp"
+#undef PFASST_UNIT_TESTING
+    
+/*
+ * parameterized test fixture with number of nodes as parameter
+ */
+class ConvergenceTest
+  : public TestWithParam<tuple<size_t, pfasst::QuadratureType>>
+{
+  protected:
+    size_t nnodes; // parameter
+
+    complex<double> lambda;
+    double Tend;
+    vector<size_t> nsteps;
+    size_t nsteps_l;
+    vector<double> err;
+    vector<double> convrate;
+    double dt;
+    size_t niters;
+    pfasst::QuadratureType nodetype;
+    size_t nnodes_in_call;
+
+    // set parameters base on node type
+    void set_parameters()
+    {
+      switch (this->nodetype)
+      {
+        case pfasst::QuadratureType::GaussLobatto:
+          this->niters = 2 * this->nnodes - 2;
+          this->Tend = 4.0;
+          this->lambda = complex<double>(-1.0, 1.0);
+          this->nsteps = { 2, 5, 10, 15, 20 };
+          this->nnodes_in_call = this->nnodes;
+          break;
+
+        case pfasst::QuadratureType::GaussLegendre:
+          this->niters = 2 * this->nnodes;
+          this->Tend = 6.0;
+          this->lambda = complex<double>(-1.0, 2.0);
+          this->nsteps = { 2, 4, 6, 8, 10 };
+          this->nnodes_in_call = this->nnodes + 2;
+          break;
+        
+        case pfasst::QuadratureType::GaussRadau:
+          this->niters = 2 * this->nnodes - 1;
+          this->Tend = 5.0;
+          this->lambda = complex<double>(-1.0, 2.0);
+          this->nsteps = { 4, 6, 8, 10, 12 };
+          this->nnodes_in_call = this->nnodes + 1;
+          break;
+         
+        // NOTE: At the moment, both Clenshaw Curtis and equidistant nodes do not
+        // reproduce the expected convergence rate... something is wrong, either with
+        // the test or with the nodes
+        
+        // Also: What is the ACTUAL number of quadrature nodes in both cases?
+        case pfasst::QuadratureType::ClenshawCurtis:
+          this->niters = this->nnodes;
+          this->Tend = 1.0;
+          this->lambda = complex<double>(-1.0, 1.0);
+          this->nsteps = {3, 5, 7, 9, 11};
+          this->nnodes_in_call = this->nnodes + 1;
+          break;        
+              
+        case pfasst::QuadratureType::Uniform:
+          this->niters = this->nnodes;
+          this->Tend = 5.0;
+          this->lambda = complex<double>(-1.0, 1.0);
+          this->nsteps = {3, 5, 7, 9, 11};
+          this->nnodes_in_call = this->nnodes;
+          break;
+         
+
+        default:
+          break;
+      }
+    }
+
+  public:
+    virtual void SetUp()
+    {
+      this->nnodes = get<0>(GetParam());
+      this->nodetype = get<1>(GetParam());
+      this->set_parameters();
+      this->nsteps_l = this->nsteps.size();
+      this->err.resize(this->nsteps.size());
+      this->convrate.resize(this->nsteps.size());
+
+      // run to compute errors
+      for (size_t i = 0; i <= this->nsteps_l - 1; ++i) {
+        this->dt = this->Tend / double(this->nsteps[i]);
+        this->err[i] = run_scalar_sdc(this->nsteps[i], this->dt, this->nnodes_in_call, 
+                                      this->niters, this->lambda, this->nodetype);
+      }
+
+      // compute convergence rates
+      for (size_t i = 0; i <= this->nsteps_l - 2; ++i) {
+        this->convrate[i] = log10(this->err[i+1] / this->err[i]) / 
+                              log10(double(this->nsteps[i]) / double(this->nsteps[i + 1]));
+      }
+    }
+
+    virtual void TearDown()
+    {}
+};
+
+/*
+ * The test below verifies that the code approximately (up to a safety factor) reproduces
+ * the theoretically expected rate of convergence
+ */
+TEST_P(ConvergenceTest, AllNodes)
+{
+  for (size_t i = 0; i <= nsteps_l - 2; ++i) {
+    switch (this->nodetype)
+    {
+      case pfasst::QuadratureType::GaussLobatto:
+        // Expect convergence rate of 2*nodes-2 from collocation formula, doing an identical number
+        // of iteration should suffice to reach this as each iteration should increase order by one
+        EXPECT_THAT(convrate[i],
+                    DoubleNear(double(2 * nnodes - 2), 0.99)) << "Convergence rate for " 
+                                                              << this->nnodes
+                                                              << " Gauss-Lobatto nodes "
+                                                              << " at node " << i 
+                                                              << " not within expected range.";
+        break;
+
+      case pfasst::QuadratureType::GaussLegendre:
+        // convergence rates for Legendre nodes should be 2*nodes but are actually better, so
+        // use Ge here
+        EXPECT_THAT(convrate[i], Ge<double>(2 * this->nnodes)) << "Convergence rate for " 
+                                                               << this->nnodes
+                                                               << " Gauss-Legendre nodes "
+                                                               << " at node " << i 
+                                                               << " not within expected range.";
+        break;
+
+      
+      case pfasst::QuadratureType::GaussRadau:
+        // convergence rate for Radau nodes should be 2*nodes-1
+        // For some case, the convergence rate is ALMOST that value, hence put in the 0.99
+        EXPECT_THAT(convrate[i], Ge<double>(0.99* 2 * this->nnodes - 1)) << "Convergence rate for " 
+                                                                         << this->nnodes
+                                                                         << " Gauss-Radu nodes "
+                                                                         << " at node " << i 
+                                                                         << " not within expected range.";
+        break;
+      
+      case pfasst::QuadratureType::ClenshawCurtis:
+        // Clenshaw Curtis should be of order nnodes
+        EXPECT_THAT(convrate[i], Ge<double>(this->nnodes))  << "Convergence rate for " 
+                                                            << this->nnodes
+                                                            << " Clenshaw-Curtis nodes "
+                                                            << " at node " << i 
+                                                            << " not within expected range.";
+        break;
+       
+      case pfasst::QuadratureType::Uniform:
+        // Equidistant nodes should be of order nnodes
+        EXPECT_THAT(convrate[i], Ge<double>(this->nnodes)) << "Convergence rate for " 
+                                                           << this->nnodes
+                                                           << " equidistant nodes "
+                                                           << " at node " << i 
+                                                           << " not within expected range.";
+        break;
+       
+
+      default:
+        EXPECT_TRUE(false);
+        break;
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P(ScalarSDC, ConvergenceTest,
+                        Combine(Range<size_t>(2, 7),
+                                Values(pfasst::QuadratureType::GaussLobatto,
+                                       pfasst::QuadratureType::GaussLegendre,
+                                       pfasst::QuadratureType::GaussRadau/*,
+                                       pfasst::QuadratureType::ClenshawCurtis,
+                                       pfasst::QuadratureType::Uniform*/))
+);
+
+int main(int argc, char** argv)
+{
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}
+
+
diff --git a/tests/examples/scalar/test_scalar_highprec.cpp b/tests/examples/scalar/test_scalar_highprec.cpp
new file mode 100644
index 00000000..653a0924
--- /dev/null
+++ b/tests/examples/scalar/test_scalar_highprec.cpp
@@ -0,0 +1,95 @@
+/*
+ * Tests for the scalar example solving the test equation
+ */
+#include <cmath>
+
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+using namespace ::testing;
+
+#include <pfasst/quadrature.hpp>
+
+#define PFASST_UNIT_TESTING
+#include "../examples/scalar/scalar_sdc.cpp"
+#undef PFASST_UNIT_TESTING
+
+/*
+ * parameterized test that verifies that given sufficiently many nodes and iterations,
+ * SDC can reproduce the analytic solution with very high precision
+ */
+class HighPrecisionTest
+  : public TestWithParam<pfasst::QuadratureType>
+{
+  protected:
+    pfasst::QuadratureType nodetype; // parameter
+    
+    const complex<double> lambda = complex<double>(-1.0,1.0);
+    const double dt = 0.2; // = Tend for single step
+    const size_t nsteps = 1;
+    const size_t niters = 30;
+    const size_t nnodes = 8; 
+    size_t nnodes_in_call;
+    double err;
+    
+    void set_parameters()
+    {
+      switch (this->nodetype)
+      {
+        case pfasst::QuadratureType::GaussLobatto:
+          this->nnodes_in_call = this->nnodes;
+          break;
+        
+        case pfasst::QuadratureType::GaussLegendre:
+          this->nnodes_in_call = this->nnodes + 2;
+          break;
+          
+        case pfasst::QuadratureType::GaussRadau:
+          this->nnodes_in_call = this->nnodes + 1;
+          break;
+        
+        case pfasst::QuadratureType::ClenshawCurtis:
+          this->nnodes_in_call = this->nnodes + 1;
+          break;
+        
+        case pfasst::QuadratureType::Uniform:
+          this->nnodes_in_call = this->nnodes + 1;
+          break;  
+        
+        default:
+          break;
+      }
+      
+    }
+            
+  public:
+    virtual void SetUp()
+    {
+      this->nodetype = GetParam();
+      this->set_parameters();
+      this->err = run_scalar_sdc(this->nsteps, this->dt, this->nnodes_in_call,
+                                 this->niters, this->lambda, this->nodetype);
+    }
+    
+    virtual void TearDown()
+    {}
+};
+
+TEST_P(HighPrecisionTest, AllNodes)
+{
+  EXPECT_THAT(err, Le<double>(5e-12)) << "Failed to bring relative error below 5e-12";
+}
+
+INSTANTIATE_TEST_CASE_P(ScalarSDC, HighPrecisionTest,
+                                Values(pfasst::QuadratureType::GaussLobatto,
+                                       pfasst::QuadratureType::GaussLegendre,
+                                       pfasst::QuadratureType::GaussRadau,
+                                       pfasst::QuadratureType::ClenshawCurtis,
+                                       pfasst::QuadratureType::Uniform)
+);
+
+int main(int argc, char** argv)
+{
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}
\ No newline at end of file
diff --git a/tests/test-quadrature.cpp b/tests/test_quadrature.cpp
similarity index 61%
rename from tests/test-quadrature.cpp
rename to tests/test_quadrature.cpp
index 49948cf4..475a41e2 100644
--- a/tests/test-quadrature.cpp
+++ b/tests/test_quadrature.cpp
@@ -8,6 +8,8 @@
 #include <gtest/gtest.h>
 #include <gmock/gmock.h>
 
+using namespace ::testing;
+
 #include <pfasst/quadrature.hpp>
 
 using namespace std;
@@ -73,14 +75,14 @@ TEST(NodesTest, GaussLegendreNodes)
                                0.9745539561713793
                              };
 
-  auto l3 = pfasst::compute_nodes<long double>(3, "gauss-legendre");
-  EXPECT_THAT(l3, testing::Pointwise(DoubleNear(), l3e));
+  auto l3 = pfasst::compute_nodes<long double>(3, pfasst::QuadratureType::GaussLegendre);
+  EXPECT_THAT(l3, Pointwise(DoubleNear(), l3e));
 
-  auto l5 = pfasst::compute_nodes<long double>(5, "gauss-legendre");
-  EXPECT_THAT(l5, testing::Pointwise(DoubleNear(), l5e));
+  auto l5 = pfasst::compute_nodes<long double>(5, pfasst::QuadratureType::GaussLegendre);
+  EXPECT_THAT(l5, Pointwise(DoubleNear(), l5e));
 
-  auto l7 = pfasst::compute_nodes<long double>(7, "gauss-legendre");
-  EXPECT_THAT(l7, testing::Pointwise(DoubleNear(), l7e));
+  auto l7 = pfasst::compute_nodes<long double>(7, pfasst::QuadratureType::GaussLegendre);
+  EXPECT_THAT(l7, Pointwise(DoubleNear(), l7e));
 }
 
 TEST(NodesTest, GaussLobattoNodes)
@@ -121,20 +123,20 @@ TEST(NodesTest, GaussLobattoNodes)
                                1.0
                              };
 
-  auto l2 = pfasst::compute_nodes<long double>(2, "gauss-lobatto");
-  EXPECT_THAT(l2, testing::Pointwise(DoubleNear(), l2e));
+  auto l2 = pfasst::compute_nodes<long double>(2, pfasst::QuadratureType::GaussLobatto);
+  EXPECT_THAT(l2, Pointwise(DoubleNear(), l2e));
 
-  auto l3 = pfasst::compute_nodes<long double>(3, "gauss-lobatto");
-  EXPECT_THAT(l3, testing::Pointwise(DoubleNear(), l3e));
+  auto l3 = pfasst::compute_nodes<long double>(3, pfasst::QuadratureType::GaussLobatto);
+  EXPECT_THAT(l3, Pointwise(DoubleNear(), l3e));
 
-  auto l5 = pfasst::compute_nodes<long double>(5, "gauss-lobatto");
-  EXPECT_THAT(l5, testing::Pointwise(DoubleNear(), l5e));
+  auto l5 = pfasst::compute_nodes<long double>(5, pfasst::QuadratureType::GaussLobatto);
+  EXPECT_THAT(l5, Pointwise(DoubleNear(), l5e));
 
-  auto l7 = pfasst::compute_nodes<long double>(7, "gauss-lobatto");
-  EXPECT_THAT(l7, testing::Pointwise(DoubleNear(), l7e));
+  auto l7 = pfasst::compute_nodes<long double>(7, pfasst::QuadratureType::GaussLobatto);
+  EXPECT_THAT(l7, Pointwise(DoubleNear(), l7e));
 
-  auto l9 = pfasst::compute_nodes<long double>(9, "gauss-lobatto");
-  EXPECT_THAT(l9, testing::Pointwise(DoubleNear(), l9e));
+  auto l9 = pfasst::compute_nodes<long double>(9, pfasst::QuadratureType::GaussLobatto);
+  EXPECT_THAT(l9, Pointwise(DoubleNear(), l9e));
 }
 
 TEST(NodesTest, ClenshawCurtisNodes)
@@ -175,54 +177,55 @@ TEST(NodesTest, ClenshawCurtisNodes)
                                 1.0
                               };
 
-  auto cc2 = pfasst::compute_nodes<long double>(2, "clenshaw-curtis");
-  EXPECT_THAT(cc2, testing::Pointwise(DoubleNear(), cc2e));
+  auto cc2 = pfasst::compute_nodes<long double>(2, pfasst::QuadratureType::ClenshawCurtis);
+  EXPECT_THAT(cc2, Pointwise(DoubleNear(), cc2e));
 
-  auto cc3 = pfasst::compute_nodes<long double>(3, "clenshaw-curtis");
-  EXPECT_THAT(cc3, testing::Pointwise(DoubleNear(), cc3e));
+  auto cc3 = pfasst::compute_nodes<long double>(3, pfasst::QuadratureType::ClenshawCurtis);
+  EXPECT_THAT(cc3, Pointwise(DoubleNear(), cc3e));
 
-  auto cc5 = pfasst::compute_nodes<long double>(5, "clenshaw-curtis");
-  EXPECT_THAT(cc5, testing::Pointwise(DoubleNear(), cc5e));
+  auto cc5 = pfasst::compute_nodes<long double>(5, pfasst::QuadratureType::ClenshawCurtis);
+  EXPECT_THAT(cc5, Pointwise(DoubleNear(), cc5e));
 
-  auto cc7 = pfasst::compute_nodes<long double>(7, "clenshaw-curtis");
-  EXPECT_THAT(cc7, testing::Pointwise(DoubleNear(), cc7e));
+  auto cc7 = pfasst::compute_nodes<long double>(7, pfasst::QuadratureType::ClenshawCurtis);
+  EXPECT_THAT(cc7, Pointwise(DoubleNear(), cc7e));
 
-  auto cc9 = pfasst::compute_nodes<long double>(9, "clenshaw-curtis");
-  EXPECT_THAT(cc9, testing::Pointwise(DoubleNear(), cc9e));
+  auto cc9 = pfasst::compute_nodes<long double>(9, pfasst::QuadratureType::ClenshawCurtis);
+  EXPECT_THAT(cc9, Pointwise(DoubleNear(), cc9e));
 }
 
 TEST(NodesTest, UniformNodes)
 {
   const long double u2e[2] = { 0.0,
                                1.0
-                              };
+                             };
 
   const long double u3e[3] = { 0.0,
                                0.5,
                                1.0
-                              };
+                             };
 
   const long double u5e[5] = { 0.0,
                                0.25,
                                0.5,
                                0.75,
                                1.0
-                              };
+                             };
 
-  auto u2 = pfasst::compute_nodes<long double>(2, "uniform");
-  EXPECT_THAT(u2, testing::Pointwise(DoubleNear(), u2e));
+  auto u2 = pfasst::compute_nodes<long double>(2, pfasst::QuadratureType::Uniform);
+  EXPECT_THAT(u2, Pointwise(DoubleNear(), u2e));
 
-  auto u3 = pfasst::compute_nodes<long double>(3, "uniform");
-  EXPECT_THAT(u3, testing::Pointwise(DoubleNear(), u3e));
+  auto u3 = pfasst::compute_nodes<long double>(3, pfasst::QuadratureType::Uniform);
+  EXPECT_THAT(u3, Pointwise(DoubleNear(), u3e));
 
-  auto u5 = pfasst::compute_nodes<long double>(5, "uniform");
-  EXPECT_THAT(u5, testing::Pointwise(DoubleNear(), u5e));
+  auto u5 = pfasst::compute_nodes<long double>(5, pfasst::QuadratureType::Uniform);
+  EXPECT_THAT(u5, Pointwise(DoubleNear(), u5e));
 }
 
 TEST(QuadratureTest, GaussLobattoNodes)
 {
-  auto l3 = pfasst::compute_nodes<long double>(3, "gauss-lobatto");
-  auto s3 = pfasst::compute_quadrature(l3, l3, 's');
+  auto l3 = pfasst::compute_nodes<long double>(3, pfasst::QuadratureType::GaussLobatto);
+  auto a3 = pfasst::augment_nodes(l3);
+  auto s3 = pfasst::compute_quadrature(get<0>(a3), get<0>(a3), get<1>(a3), pfasst::QuadratureMatrix::S);
   const long double s3e[6] = { 0.20833333333333333,
                                0.33333333333333333,
                                -0.04166666666666666,
@@ -230,11 +233,12 @@ TEST(QuadratureTest, GaussLobattoNodes)
                                0.33333333333333333,
                                0.20833333333333333
                              };
+    
+  EXPECT_THAT(s3.data(), Pointwise(DoubleNear(), s3e));
 
-  EXPECT_THAT(s3.data(), testing::Pointwise(DoubleNear(), s3e));
-
-  auto l5 = pfasst::compute_nodes<long double>(5, "gauss-lobatto");
-  auto s5 = pfasst::compute_quadrature(l5, l5, 's');
+  auto l5 = pfasst::compute_nodes<long double>(5, pfasst::QuadratureType::GaussLobatto);
+  auto a5 = pfasst::augment_nodes(l5);
+  auto s5 = pfasst::compute_quadrature(get<0>(a5), get<0>(a5), get<1>(a5), pfasst::QuadratureMatrix::S);
   const long double s5e[] = { 0.067728432186156897969267419174073482,
                               0.11974476934341168251615379970493965,
                               -0.021735721866558113665511351745074292,
@@ -255,16 +259,17 @@ TEST(QuadratureTest, GaussLobattoNodes)
                               -0.021735721866558113665511351745074289,
                               0.11974476934341168251615379970493965,
                               0.067728432186156897969267419174073482
-                             };
-  EXPECT_THAT(s5.data(), testing::Pointwise(DoubleNear(), s5e));
+                            };
+  EXPECT_THAT(s5.data(), Pointwise(DoubleNear(), s5e));
 }
 
 TEST(QuadratureTest, ClenshawCurtisNodes)
 {
-  auto c4 = pfasst::compute_nodes<long double>(4, "clenshaw-curtis");
-  auto s4 = pfasst::compute_quadrature(c4, c4, 's');
+  auto c4 = pfasst::compute_nodes<long double>(4, pfasst::QuadratureType::ClenshawCurtis);
+  auto a4 = pfasst::augment_nodes(c4);
+  auto s4 = pfasst::compute_quadrature(get<0>(a4), get<0>(a4), get<1>(a4), pfasst::QuadratureMatrix::S);
   const long double s4e[] = { 0.10243055555555555555555555555555556,
-                               0.16319444444444444444444444444444444,
+                              0.16319444444444444444444444444444444,
                               -0.024305555555555555555555555555555556,
                               0.0086805555555555555555555555555555557,
                               -0.055555555555555555555555555555555556,
@@ -275,12 +280,61 @@ TEST(QuadratureTest, ClenshawCurtisNodes)
                               -0.024305555555555555555555555555555554,
                               0.16319444444444444444444444444444444,
                               0.10243055555555555555555555555555556
-                             };
-  EXPECT_THAT(s4.data(), testing::Pointwise(DoubleNear(), s4e));
+                            };
+  EXPECT_THAT(s4.data(), Pointwise(DoubleNear(), s4e));
+}
+
+class QmatTest
+  : public ::TestWithParam<tuple<size_t, pfasst::QuadratureType>>
+{
+  protected:
+    size_t nnodes;
+    pfasst::QuadratureType qtype;
+    vector<long double> nodes;
+    matrix<long double> q;
+
+  public:
+    virtual void SetUp()
+    {
+      nnodes = get<0>(GetParam());
+      qtype = get<1>(GetParam());
+      
+      auto cnodes = pfasst::compute_nodes<long double>(nnodes, qtype);
+      auto n1 = pfasst::augment_nodes(cnodes);
+      nodes = get<0>(n1);
+      q = pfasst::compute_quadrature(nodes, nodes, get<1>(n1), pfasst::QuadratureMatrix::Q);
+      
+      cout << "Quadrature type no. " << int(qtype) << " -- Number of nodes " << nnodes << endl;
+      
+    }
+  
+    virtual void TearDown()
+    {}
+};
+
+TEST_P(QmatTest, AllNodes)
+{
+  long double qsum;
+  for (size_t m = 0; m < q.size1(); ++m) {
+    qsum = 0;
+    for (size_t j = 0; j < q.size2(); ++j) {
+      qsum += q(m,j);
+    }
+    EXPECT_NEAR(qsum, nodes[m+1], (long double)(3E-12));
+  }
 }
 
+INSTANTIATE_TEST_CASE_P(Quadrature, QmatTest,
+                        ::Combine(::Range<size_t>(2, 14),
+                                  Values<pfasst::QuadratureType>(pfasst::QuadratureType::GaussLegendre,
+                                                                 pfasst::QuadratureType::GaussLobatto,
+                                                                 pfasst::QuadratureType::GaussRadau,
+                                                                 pfasst::QuadratureType::ClenshawCurtis,
+                                                                 pfasst::QuadratureType::Uniform)));
+
+
 int main(int argc, char** argv)
 {
-  testing::InitGoogleTest(&argc, argv);
+  InitGoogleTest(&argc, argv);
   return RUN_ALL_TESTS();
 }