From c4dedb88840349237f61c01f696e77588d3747d0 Mon Sep 17 00:00:00 2001
From: David Blom <d.s.blom@tudelft.nl>
Date: Thu, 1 Sep 2016 16:32:59 +0200
Subject: [PATCH] Use El::VR, El::STAR distribution for elemental matrices
 (#341)

---
 src/RBFMeshMotionSolver/AdaptiveCoarsening.C  | 42 +++++++++----------
 src/RBFMeshMotionSolver/AdaptiveCoarsening.H  | 18 ++++----
 src/RBFMeshMotionSolver/Coarsener.H           |  6 +--
 src/RBFMeshMotionSolver/ElRBFInterpolation.C  | 12 +++---
 src/RBFMeshMotionSolver/ElRBFInterpolation.H  | 12 +++---
 .../ElRBFMeshMotionSolver.C                   |  8 ++--
 src/RBFMeshMotionSolver/NoCoarsening.C        |  6 +--
 src/RBFMeshMotionSolver/NoCoarsening.H        |  6 +--
 src/RBFMeshMotionSolver/UnitCoarsening.C      | 34 +++++++--------
 src/RBFMeshMotionSolver/UnitCoarsening.H      | 10 ++---
 src/fsi/AdaptiveTimeStepper.C                 |  1 +
 src/fsi/SDC.C                                 |  1 +
 .../testsuite-rbf/test_adaptivecoarsening.C   | 24 +++++------
 .../testsuite-rbf/test_elrbfinterpolation.C   | 38 ++++++++---------
 src/tests/testsuite-rbf/test_nocoarsener.C    | 16 +++----
 src/tests/testsuite-rbf/test_unitcoarsening.C | 28 ++++++-------
 16 files changed, 133 insertions(+), 129 deletions(-)

diff --git a/src/RBFMeshMotionSolver/AdaptiveCoarsening.C b/src/RBFMeshMotionSolver/AdaptiveCoarsening.C
index b030bebe..10639797 100644
--- a/src/RBFMeshMotionSolver/AdaptiveCoarsening.C
+++ b/src/RBFMeshMotionSolver/AdaptiveCoarsening.C
@@ -30,8 +30,8 @@ namespace rbf
 
     void AdaptiveCoarsening::compute(
         std::shared_ptr<RBFFunctionInterface> function,
-        std::unique_ptr<El::DistMatrix<double> > pos,
-        std::unique_ptr<El::DistMatrix<double> > posInterpolation
+        std::unique_ptr<ElDistVector> pos,
+        std::unique_ptr<ElDistVector> posInterpolation
         )
     {
         // Selection is performed as necessary. It is only performed when interpolate()
@@ -45,24 +45,24 @@ namespace rbf
         this->positionsInterpolation = std::move( posInterpolation );
     }
 
-    std::pair<int, double> AdaptiveCoarsening::computeError( const std::unique_ptr<El::DistMatrix<double> > & values )
+    std::pair<int, double> AdaptiveCoarsening::computeError( const std::unique_ptr<ElDistVector> & values )
     {
         // Select a subset of values based on the selected points
 
-        std::unique_ptr<El::DistMatrix<double> > valuesCoarse( new El::DistMatrix<double>() );
+        std::unique_ptr<ElDistVector> valuesCoarse( new ElDistVector() );
         El::Zeros( *valuesCoarse, selectedPositions.size(), values->Width() );
         selectData( values, valuesCoarse );
 
         // Perform the interpolation
 
-        std::unique_ptr<El::DistMatrix<double> > result = rbfCoarse->interpolate( valuesCoarse );
+        std::unique_ptr<ElDistVector> result = rbfCoarse->interpolate( valuesCoarse );
 
         assert( values->Height() == result->Height() );
 
         // Compute error
-        El::DistMatrix<double> diff = *values;
+        ElDistVector diff = *values;
         El::Axpy( -1, *result, diff );
-        El::DistMatrix<double, El::MC, El::STAR> errors;
+        ElDistVector errors;
         El::RowTwoNorms( diff, errors );
 
         // Get location of max error
@@ -78,7 +78,7 @@ namespace rbf
             return std::pair<int, double>( locMax.i, locMax.value );
     }
 
-    void AdaptiveCoarsening::greedySelection( const std::unique_ptr<El::DistMatrix<double> > & values )
+    void AdaptiveCoarsening::greedySelection( const std::unique_ptr<ElDistVector> & values )
     {
         selectedPositions.clear();
 
@@ -90,19 +90,19 @@ namespace rbf
 
         {
             // Find first point: largest value
-            El::DistMatrix<double, El::MC, El::STAR> norms;
+            ElDistVector norms;
             El::RowTwoNorms( *values, norms );
             El::Entry<double> locMax = El::MaxAbsLoc( norms );
             selectedPositions.push_back( locMax.i );
 
             // Find second point: largest distance from the first point
-            El::DistMatrix<double> distance = *positions;
-            El::DistMatrix<double> tmp;
+            ElDistVector distance = *positions;
+            ElDistVector tmp;
             El::Ones( tmp, distance.Height(), distance.Width() );
 
             for ( int iColumn = 0; iColumn < tmp.Width(); iColumn++ )
             {
-                El::DistMatrix<double> view;
+                ElDistVector view;
                 El::View( view, tmp, 0, iColumn, tmp.Height(), 1 );
                 El::Scale( positions->Get( locMax.i, iColumn ), view );
             }
@@ -121,12 +121,12 @@ namespace rbf
         for ( int i = 0; i < maxPoints; i++ )
         {
             // Build the matrices for the RBF interpolation
-            std::unique_ptr<El::DistMatrix<double> > positionsCoarse( new El::DistMatrix<double>() );
+            std::unique_ptr<ElDistVector> positionsCoarse( new ElDistVector() );
             El::Zeros( *positionsCoarse, selectedPositions.size(), positions->Width() );
             selectData( positions, positionsCoarse );
 
             // Perform the RBF interpolation
-            std::unique_ptr<El::DistMatrix<double> > positionsInterpolationCoarse( new El::DistMatrix<double>() );
+            std::unique_ptr<ElDistVector> positionsInterpolationCoarse( new ElDistVector() );
             *positionsInterpolationCoarse = *positions;
             rbfCoarse = std::unique_ptr<ElRBFInterpolation>( new ElRBFInterpolation( rbfFunction, std::move( positionsCoarse ), std::move( positionsInterpolationCoarse ) ) );
 
@@ -154,14 +154,14 @@ namespace rbf
 
         // Initialize interpolator
 
-        std::unique_ptr<El::DistMatrix<double> > positionsCoarse( new El::DistMatrix<double>() );
+        std::unique_ptr<ElDistVector> positionsCoarse( new ElDistVector() );
         El::Zeros( *positionsCoarse, selectedPositions.size(), positions->Width() );
 
         selectData( positions, positionsCoarse );
 
         rbf = std::unique_ptr<ElRBFInterpolation>( new ElRBFInterpolation() );
 
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolationTmp( new El::DistMatrix<double>( *positionsInterpolation ) );
+        std::unique_ptr<ElDistVector> positionsInterpolationTmp( new ElDistVector( *positionsInterpolation ) );
 
         rbf->compute( rbfFunction, std::move( positionsCoarse ), std::move( positionsInterpolationTmp ) );
     }
@@ -171,7 +171,7 @@ namespace rbf
         return rbf->initialized();
     }
 
-    std::unique_ptr<El::DistMatrix<double> > AdaptiveCoarsening::interpolate( const std::unique_ptr<El::DistMatrix<double> > & values )
+    std::unique_ptr<ElDistVector> AdaptiveCoarsening::interpolate( const std::unique_ptr<ElDistVector> & values )
     {
         values->ProcessQueues();
 
@@ -190,7 +190,7 @@ namespace rbf
             }
             else
             {
-                std::unique_ptr<El::DistMatrix<double> > result( new El::DistMatrix<double>() );
+                std::unique_ptr<ElDistVector> result( new ElDistVector() );
                 El::Zeros( *result, positionsInterpolation->Height(), positionsInterpolation->Width() );
                 return result;
             }
@@ -219,7 +219,7 @@ namespace rbf
             greedySelection( values );
         }
 
-        std::unique_ptr<El::DistMatrix<double> > selectedValues( new El::DistMatrix<double>() );
+        std::unique_ptr<ElDistVector> selectedValues( new ElDistVector() );
         El::Zeros( *selectedValues, selectedPositions.size(), values->Width() );
 
         selectData( values, selectedValues );
@@ -228,8 +228,8 @@ namespace rbf
     }
 
     void AdaptiveCoarsening::selectData(
-        const std::unique_ptr<El::DistMatrix<double> > & data,
-        std::unique_ptr<El::DistMatrix<double> > & selection
+        const std::unique_ptr<ElDistVector> & data,
+        std::unique_ptr<ElDistVector> & selection
         )
     {
         assert( selection->Height() == int( selectedPositions.size() ) );
diff --git a/src/RBFMeshMotionSolver/AdaptiveCoarsening.H b/src/RBFMeshMotionSolver/AdaptiveCoarsening.H
index 81f09b11..58cfb903 100644
--- a/src/RBFMeshMotionSolver/AdaptiveCoarsening.H
+++ b/src/RBFMeshMotionSolver/AdaptiveCoarsening.H
@@ -24,22 +24,22 @@ namespace rbf
 
             void compute(
                 std::shared_ptr<RBFFunctionInterface> rbfFunction,
-                std::unique_ptr<El::DistMatrix<double> > positions,
-                std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+                std::unique_ptr<ElDistVector> positions,
+                std::unique_ptr<ElDistVector> positionsInterpolation
                 );
 
-            void greedySelection( const std::unique_ptr<El::DistMatrix<double> > & values );
+            void greedySelection( const std::unique_ptr<ElDistVector> & values );
 
             bool initialized();
 
-            std::unique_ptr<El::DistMatrix<double> > interpolate( const std::unique_ptr<El::DistMatrix<double> > & values );
+            std::unique_ptr<ElDistVector> interpolate( const std::unique_ptr<ElDistVector> & values );
 
         private:
-            std::pair<int, double> computeError( const std::unique_ptr<El::DistMatrix<double> > & values );
+            std::pair<int, double> computeError( const std::unique_ptr<ElDistVector> & values );
 
             void selectData(
-                const std::unique_ptr<El::DistMatrix<double> > & data,
-                std::unique_ptr<El::DistMatrix<double> > & selection
+                const std::unique_ptr<ElDistVector> & data,
+                std::unique_ptr<ElDistVector> & selection
                 );
 
             const double tol;
@@ -51,7 +51,7 @@ namespace rbf
             std::vector<size_t> selectedPositions;
 
             std::shared_ptr<RBFFunctionInterface> rbfFunction;
-            std::unique_ptr<El::DistMatrix<double> > positions;
-            std::unique_ptr<El::DistMatrix<double> > positionsInterpolation;
+            std::unique_ptr<ElDistVector> positions;
+            std::unique_ptr<ElDistVector> positionsInterpolation;
     };
 }
diff --git a/src/RBFMeshMotionSolver/Coarsener.H b/src/RBFMeshMotionSolver/Coarsener.H
index b0057cda..27aa7c0b 100644
--- a/src/RBFMeshMotionSolver/Coarsener.H
+++ b/src/RBFMeshMotionSolver/Coarsener.H
@@ -12,12 +12,12 @@ namespace rbf
 
             virtual void compute(
                 std::shared_ptr<RBFFunctionInterface> rbfFunction,
-                std::unique_ptr<El::DistMatrix<double> > positions,
-                std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+                std::unique_ptr<ElDistVector> positions,
+                std::unique_ptr<ElDistVector> positionsInterpolation
                 ) = 0;
 
             virtual bool initialized() = 0;
 
-            virtual std::unique_ptr<El::DistMatrix<double> > interpolate( const std::unique_ptr<El::DistMatrix<double> > & values ) = 0;
+            virtual std::unique_ptr<ElDistVector> interpolate( const std::unique_ptr<ElDistVector> & values ) = 0;
     };
 }
diff --git a/src/RBFMeshMotionSolver/ElRBFInterpolation.C b/src/RBFMeshMotionSolver/ElRBFInterpolation.C
index 30a60e48..23b058f3 100644
--- a/src/RBFMeshMotionSolver/ElRBFInterpolation.C
+++ b/src/RBFMeshMotionSolver/ElRBFInterpolation.C
@@ -17,8 +17,8 @@ namespace rbf
 
     ElRBFInterpolation::ElRBFInterpolation(
         std::shared_ptr<RBFFunctionInterface> rbfFunction,
-        std::unique_ptr<El::DistMatrix<double> > positions,
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+        std::unique_ptr<ElDistVector> positions,
+        std::unique_ptr<ElDistVector> positionsInterpolation
         )
         :
         H( new El::DistMatrix<double>() ),
@@ -32,8 +32,8 @@ namespace rbf
 
     void ElRBFInterpolation::compute(
         std::shared_ptr<RBFFunctionInterface> rbfFunction,
-        std::unique_ptr<El::DistMatrix<double> > positions,
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+        std::unique_ptr<ElDistVector> positions,
+        std::unique_ptr<ElDistVector> positionsInterpolation
         )
     {
         assert( Phi->Height() == 0 );
@@ -139,13 +139,13 @@ namespace rbf
         return Phi->Height() > 0;
     }
 
-    std::unique_ptr<El::DistMatrix<double> > ElRBFInterpolation::interpolate( const std::unique_ptr<El::DistMatrix<double> > & values )
+    std::unique_ptr<ElDistVector> ElRBFInterpolation::interpolate( const std::unique_ptr<ElDistVector> & values )
     {
         assert( Phi->Height() > 0 );
         assert( H->Height() > 0 );
         assert( values->Height() == Phi->Width() );
 
-        std::unique_ptr<El::DistMatrix<double> > result( new El::DistMatrix<double>() );
+        std::unique_ptr<ElDistVector> result( new ElDistVector() );
 
         values->ProcessQueues();
 
diff --git a/src/RBFMeshMotionSolver/ElRBFInterpolation.H b/src/RBFMeshMotionSolver/ElRBFInterpolation.H
index a404e79b..87721976 100644
--- a/src/RBFMeshMotionSolver/ElRBFInterpolation.H
+++ b/src/RBFMeshMotionSolver/ElRBFInterpolation.H
@@ -12,6 +12,8 @@
 
 namespace rbf
 {
+    typedef El::DistMatrix<double, El::VR, El::STAR> ElDistVector;
+
     class ElRBFInterpolation
     {
         public:
@@ -19,21 +21,21 @@ namespace rbf
 
             explicit ElRBFInterpolation(
                 std::shared_ptr<RBFFunctionInterface> rbfFunction,
-                std::unique_ptr<El::DistMatrix<double> > positions,
-                std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+                std::unique_ptr<ElDistVector> positions,
+                std::unique_ptr<ElDistVector> positionsInterpolation
                 );
 
             ~ElRBFInterpolation();
 
             void compute(
                 std::shared_ptr<RBFFunctionInterface> rbfFunction,
-                std::unique_ptr<El::DistMatrix<double> > positions,
-                std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+                std::unique_ptr<ElDistVector> positions,
+                std::unique_ptr<ElDistVector> positionsInterpolation
                 );
 
             bool initialized();
 
-            std::unique_ptr<El::DistMatrix<double> > interpolate( const std::unique_ptr<El::DistMatrix<double> > & values );
+            std::unique_ptr<ElDistVector> interpolate( const std::unique_ptr<ElDistVector> & values );
 
         private:
             std::unique_ptr<El::DistMatrix<double> > H;
diff --git a/src/RBFMeshMotionSolver/ElRBFMeshMotionSolver.C b/src/RBFMeshMotionSolver/ElRBFMeshMotionSolver.C
index 07263231..ca048537 100644
--- a/src/RBFMeshMotionSolver/ElRBFMeshMotionSolver.C
+++ b/src/RBFMeshMotionSolver/ElRBFMeshMotionSolver.C
@@ -317,8 +317,8 @@ void ElRBFMeshMotionSolver::solve()
 
         boundaryData.resize( std::distance( boundaryData.begin(), it ) );
 
-        std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double> () );
+        std::unique_ptr<rbf::ElDistVector> positions( new rbf::ElDistVector() );
+        std::unique_ptr<rbf::ElDistVector> positionsInterpolation( new rbf::ElDistVector() );
 
         std::vector<size_t> allBoundaryPoints = mxx::allgather( boundaryData.size() );
         size_t totalNbBoundaryPoints = 0;
@@ -442,7 +442,7 @@ void ElRBFMeshMotionSolver::solve()
     for ( size_t n : allBoundaryPoints )
         totalNbBoundaryPoints += n;
 
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<rbf::ElDistVector> data( new rbf::ElDistVector() );
     El::Zeros( *data, totalNbBoundaryPoints, mesh().nGeometricD() );
 
     data->Reserve( boundaryPoints.size() * mesh().nGeometricD() );
@@ -467,7 +467,7 @@ void ElRBFMeshMotionSolver::solve()
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf->interpolate( data );
+    std::unique_ptr<rbf::ElDistVector> result = rbf->interpolate( data );
 
     vectorField valuesInterpolationField( mesh().points().size(), Foam::vector::zero );
 
diff --git a/src/RBFMeshMotionSolver/NoCoarsening.C b/src/RBFMeshMotionSolver/NoCoarsening.C
index 44674892..50589a57 100644
--- a/src/RBFMeshMotionSolver/NoCoarsening.C
+++ b/src/RBFMeshMotionSolver/NoCoarsening.C
@@ -17,8 +17,8 @@ namespace rbf
 
     void NoCoarsening::compute(
         std::shared_ptr<RBFFunctionInterface> rbfFunction,
-        std::unique_ptr<El::DistMatrix<double> > positions,
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+        std::unique_ptr<ElDistVector> positions,
+        std::unique_ptr<ElDistVector> positionsInterpolation
         )
     {
         rbf->compute( rbfFunction, std::move( positions ), std::move( positionsInterpolation ) );
@@ -29,7 +29,7 @@ namespace rbf
         return rbf->initialized();
     }
 
-    std::unique_ptr<El::DistMatrix<double> > NoCoarsening::interpolate( const std::unique_ptr<El::DistMatrix<double> > & values )
+    std::unique_ptr<ElDistVector> NoCoarsening::interpolate( const std::unique_ptr<ElDistVector> & values )
     {
         return rbf->interpolate( values );
     }
diff --git a/src/RBFMeshMotionSolver/NoCoarsening.H b/src/RBFMeshMotionSolver/NoCoarsening.H
index f087561f..8a9ada9f 100644
--- a/src/RBFMeshMotionSolver/NoCoarsening.H
+++ b/src/RBFMeshMotionSolver/NoCoarsening.H
@@ -19,13 +19,13 @@ namespace rbf
 
             void compute(
                 std::shared_ptr<RBFFunctionInterface> rbfFunction,
-                std::unique_ptr<El::DistMatrix<double> > positions,
-                std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+                std::unique_ptr<ElDistVector> positions,
+                std::unique_ptr<ElDistVector> positionsInterpolation
                 );
 
             bool initialized();
 
-            std::unique_ptr<El::DistMatrix<double> > interpolate( const std::unique_ptr<El::DistMatrix<double> > & values );
+            std::unique_ptr<ElDistVector> interpolate( const std::unique_ptr<ElDistVector> & values );
 
         private:
             std::unique_ptr<ElRBFInterpolation> rbf;
diff --git a/src/RBFMeshMotionSolver/UnitCoarsening.C b/src/RBFMeshMotionSolver/UnitCoarsening.C
index 3aa4cfda..9040415a 100644
--- a/src/RBFMeshMotionSolver/UnitCoarsening.C
+++ b/src/RBFMeshMotionSolver/UnitCoarsening.C
@@ -28,8 +28,8 @@ namespace rbf
 
     void UnitCoarsening::compute(
         std::shared_ptr<RBFFunctionInterface> rbfFunction,
-        std::unique_ptr<El::DistMatrix<double> > positions,
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+        std::unique_ptr<ElDistVector> positions,
+        std::unique_ptr<ElDistVector> positionsInterpolation
         )
     {
         // selection is only performed once. Thereafter, the selected points are used for the
@@ -45,19 +45,19 @@ namespace rbf
 
         {
             // Find first point: largest radius from origin
-            El::DistMatrix<double, El::MC, El::STAR> norms;
+            ElDistVector norms;
             El::RowTwoNorms( *positions, norms );
             El::Entry<double> locMax = El::MaxAbsLoc( norms );
             selectedPositions.push_back( locMax.i );
 
             // Find second point: largest distance from the first point
-            El::DistMatrix<double> distance = *positions;
-            El::DistMatrix<double> tmp;
+            ElDistVector distance = *positions;
+            ElDistVector tmp;
             El::Ones( tmp, distance.Height(), distance.Width() );
 
             for ( int iColumn = 0; iColumn < tmp.Width(); iColumn++ )
             {
-                El::DistMatrix<double> view;
+                ElDistVector view;
                 El::View( view, tmp, 0, iColumn, tmp.Height(), 1 );
                 El::Scale( positions->Get( locMax.i, iColumn ), view );
             }
@@ -77,30 +77,30 @@ namespace rbf
         for ( int i = 0; i < maxPoints; i++ )
         {
             // Build the matrices for the RBF interpolation
-            std::unique_ptr<El::DistMatrix<double> > positionsCoarse( new El::DistMatrix<double>() );
-            std::unique_ptr<El::DistMatrix<double> > valuesCoarse( new El::DistMatrix<double>() );
+            std::unique_ptr<ElDistVector> positionsCoarse( new ElDistVector() );
+            std::unique_ptr<ElDistVector> valuesCoarse( new ElDistVector() );
             El::Ones( *valuesCoarse, selectedPositions.size(), positions->Width() );
             El::Zeros( *positionsCoarse, selectedPositions.size(), positions->Width() );
 
             selectData( positions, positionsCoarse );
 
             // Perform the RBF interpolation
-            std::unique_ptr<El::DistMatrix<double> > positionsInterpolationCoarse( new El::DistMatrix<double>() );
+            std::unique_ptr<ElDistVector> positionsInterpolationCoarse( new ElDistVector() );
 
             *positionsInterpolationCoarse = *positions;
 
             ElRBFInterpolation rbf( rbfFunction, std::move( positionsCoarse ), std::move( positionsInterpolationCoarse ) );
 
-            std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( valuesCoarse ) );
+            std::unique_ptr<ElDistVector> result = rbf.interpolate( std::move( valuesCoarse ) );
 
             assert( result->Height() == positions->Height() );
             assert( result->Width() == positions->Width() );
 
             // Compute the error
-            El::DistMatrix<double> diff;
+            ElDistVector diff;
             El::Ones( diff, result->Height(), result->Width() );
             El::Axpy( -1, *result, diff );
-            El::DistMatrix<double, El::MC, El::STAR> errors;
+            ElDistVector errors;
             El::RowTwoNorms( diff, errors );
 
             // Get location of max error
@@ -126,7 +126,7 @@ namespace rbf
             std::cout << " error = " << largestError << ", tol = " << tol << std::endl;
         }
 
-        std::unique_ptr<El::DistMatrix<double> > positionsCoarse( new El::DistMatrix<double>() );
+        std::unique_ptr<ElDistVector> positionsCoarse( new ElDistVector() );
         El::Zeros( *positionsCoarse, selectedPositions.size(), positions->Width() );
 
         selectData( positions, positionsCoarse );
@@ -139,9 +139,9 @@ namespace rbf
         return rbf.initialized();
     }
 
-    std::unique_ptr<El::DistMatrix<double> > UnitCoarsening::interpolate( const std::unique_ptr<El::DistMatrix<double> > & values )
+    std::unique_ptr<ElDistVector> UnitCoarsening::interpolate( const std::unique_ptr<ElDistVector> & values )
     {
-        std::unique_ptr<El::DistMatrix<double> > selectedValues( new El::DistMatrix<double>() );
+        std::unique_ptr<ElDistVector> selectedValues( new ElDistVector() );
         El::Zeros( *selectedValues, selectedPositions.size(), values->Width() );
 
         selectData( values, selectedValues );
@@ -150,8 +150,8 @@ namespace rbf
     }
 
     void UnitCoarsening::selectData(
-        const std::unique_ptr<El::DistMatrix<double> > & data,
-        std::unique_ptr<El::DistMatrix<double> > & selection
+        const std::unique_ptr<ElDistVector> & data,
+        std::unique_ptr<ElDistVector> & selection
         )
     {
         assert( selection->Height() == int( selectedPositions.size() ) );
diff --git a/src/RBFMeshMotionSolver/UnitCoarsening.H b/src/RBFMeshMotionSolver/UnitCoarsening.H
index ff159dd7..6da0a0bb 100644
--- a/src/RBFMeshMotionSolver/UnitCoarsening.H
+++ b/src/RBFMeshMotionSolver/UnitCoarsening.H
@@ -19,18 +19,18 @@ namespace rbf
 
             void compute(
                 std::shared_ptr<RBFFunctionInterface> rbfFunction,
-                std::unique_ptr<El::DistMatrix<double> > positions,
-                std::unique_ptr<El::DistMatrix<double> > positionsInterpolation
+                std::unique_ptr<ElDistVector> positions,
+                std::unique_ptr<ElDistVector> positionsInterpolation
                 );
 
             bool initialized();
 
-            std::unique_ptr<El::DistMatrix<double> > interpolate( const std::unique_ptr<El::DistMatrix<double> > & values );
+            std::unique_ptr<ElDistVector> interpolate( const std::unique_ptr<ElDistVector> & values );
 
         private:
             void selectData(
-                const std::unique_ptr<El::DistMatrix<double> > & data,
-                std::unique_ptr<El::DistMatrix<double> > & selection
+                const std::unique_ptr<ElDistVector> & data,
+                std::unique_ptr<ElDistVector> & selection
                 );
 
             const double tol;
diff --git a/src/fsi/AdaptiveTimeStepper.C b/src/fsi/AdaptiveTimeStepper.C
index 90c01984..30deb4fb 100644
--- a/src/fsi/AdaptiveTimeStepper.C
+++ b/src/fsi/AdaptiveTimeStepper.C
@@ -5,6 +5,7 @@
  */
 
 #include "AdaptiveTimeStepper.H"
+#include "PstreamReduceOps.H"
 
 using namespace sdc;
 
diff --git a/src/fsi/SDC.C b/src/fsi/SDC.C
index 9b668253..7a768359 100644
--- a/src/fsi/SDC.C
+++ b/src/fsi/SDC.C
@@ -4,6 +4,7 @@
  *   David Blom, TU Delft. All rights reserved.
  */
 
+#include "PstreamReduceOps.H"
 #include "SDC.H"
 #include "GaussRadau.H"
 #include "GaussRadau.H"
diff --git a/src/tests/testsuite-rbf/test_adaptivecoarsening.C b/src/tests/testsuite-rbf/test_adaptivecoarsening.C
index 4c6ffdbc..301ffc38 100644
--- a/src/tests/testsuite-rbf/test_adaptivecoarsening.C
+++ b/src/tests/testsuite-rbf/test_adaptivecoarsening.C
@@ -16,9 +16,9 @@ using namespace rbf;
 TEST( AdaptiveCoarsener, interpolate )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 8, 1 );
     El::Zeros( *positions, 4, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -73,15 +73,15 @@ TEST( AdaptiveCoarsener, interpolate )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
 }
 
 TEST( AdaptiveCoarsener, verify_Eigen_coarsening )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 20, 1 );
     El::Zeros( *positions, 10, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -228,7 +228,7 @@ TEST( AdaptiveCoarsener, verify_Eigen_coarsening )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
 
     std::vector<double> buffer;
     result->ReservePulls( result->Height() * result->Width() );
@@ -254,9 +254,9 @@ TEST( AdaptiveCoarsener, verify_Eigen_coarsening )
 TEST( AdaptiveCoarsener, 2d_grid )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positions, 20, 2 );
     El::Zeros( *positionsInterpolation, 50, 2 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -368,7 +368,7 @@ TEST( AdaptiveCoarsener, 2d_grid )
 
     rbf.compute( std::move( rbfFunction ), std::move( positions ), std::move( positionsInterpolation ) );
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( data );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( data );
 
     bufferPos.clear();
     data->ReservePulls( data->Height() * data->Width() );
diff --git a/src/tests/testsuite-rbf/test_elrbfinterpolation.C b/src/tests/testsuite-rbf/test_elrbfinterpolation.C
index f25e08dd..c1dd8898 100644
--- a/src/tests/testsuite-rbf/test_elrbfinterpolation.C
+++ b/src/tests/testsuite-rbf/test_elrbfinterpolation.C
@@ -14,8 +14,8 @@ using namespace rbf;
 TEST( ElRBFInterpolation, constructor )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>( 4, 1 ) );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>( 8, 1 ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>( 4, 1 ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>( 8, 1 ) );
 
     const double dx = 1.0 / (positions->Height() - 1);
     const double dy = 1.0 / (positionsInterpolation->Height() - 1);
@@ -52,9 +52,9 @@ TEST( ElRBFInterpolation, constructor )
 TEST( ElRBFInterpolation, interpolate )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 8, 1 );
     El::Zeros( *positions, 4, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -107,7 +107,7 @@ TEST( ElRBFInterpolation, interpolate )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
 
     for ( int i = 0; i < result->Height(); i++ )
     {
@@ -120,9 +120,9 @@ TEST( ElRBFInterpolation, interpolate )
 TEST( ElRBFInterpolation, verify_Eigen )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 8, 1 );
     El::Zeros( *positions, 4, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -266,7 +266,7 @@ TEST( ElRBFInterpolation, verify_Eigen )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
 
     std::vector<double> buffer;
     result->ReservePulls( result->Height() * result->Width() );
@@ -292,9 +292,9 @@ TEST( ElRBFInterpolation, verify_Eigen )
 TEST( ElRBFInterpolation, interpolate_twice )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 8, 1 );
     El::Zeros( *positions, 4, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -438,7 +438,7 @@ TEST( ElRBFInterpolation, interpolate_twice )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( data );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( data );
 
     std::vector<double> buffer;
     result->ReservePulls( result->Height() * result->Width() );
@@ -460,7 +460,7 @@ TEST( ElRBFInterpolation, interpolate_twice )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result2 = rbf.interpolate( data );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result2 = rbf.interpolate( data );
 
     buffer.clear();
     result->ReservePulls( result2->Height() * result2->Width() );
@@ -486,9 +486,9 @@ TEST( ElRBFInterpolation, interpolate_twice )
 TEST( ElRBFInterpolation, 2d_grid )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positions, 5, 2 );
     El::Zeros( *positionsInterpolation, 5, 2 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -598,7 +598,7 @@ TEST( ElRBFInterpolation, 2d_grid )
 
     ElRBFInterpolation rbf( std::move( rbfFunction ), std::move( positions ), std::move( positionsInterpolation ) );
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( data );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( data );
 
     bufferPos.clear();
     data->ReservePulls( data->Height() * data->Width() );
diff --git a/src/tests/testsuite-rbf/test_nocoarsener.C b/src/tests/testsuite-rbf/test_nocoarsener.C
index 5869b49c..df243dc5 100644
--- a/src/tests/testsuite-rbf/test_nocoarsener.C
+++ b/src/tests/testsuite-rbf/test_nocoarsener.C
@@ -15,9 +15,9 @@ using namespace rbf;
 TEST( NoCoarsener, interpolate )
 {
     std::shared_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 8, 1 );
     El::Zeros( *positions, 4, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -70,13 +70,13 @@ TEST( NoCoarsener, interpolate )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
-    std::unique_ptr<El::DistMatrix<double> > result2;
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result2;
 
     {
-        std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-        std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-        std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+        std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+        std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+        std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
         El::Zeros( *positionsInterpolation, 8, 1 );
         El::Zeros( *positions, 4, 1 );
         El::Zeros( *data, positions->Height(), 1 );
diff --git a/src/tests/testsuite-rbf/test_unitcoarsening.C b/src/tests/testsuite-rbf/test_unitcoarsening.C
index 0733195b..8d0d8b81 100644
--- a/src/tests/testsuite-rbf/test_unitcoarsening.C
+++ b/src/tests/testsuite-rbf/test_unitcoarsening.C
@@ -16,9 +16,9 @@ using namespace rbf;
 TEST( UnitCoarsener, interpolate )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 8, 1 );
     El::Zeros( *positions, 4, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -73,15 +73,15 @@ TEST( UnitCoarsener, interpolate )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
 }
 
 TEST( UnitCoarsener, verify_Eigen_coarsening )
 {
     std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positionsInterpolation, 20, 1 );
     El::Zeros( *positions, 20, 1 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -228,7 +228,7 @@ TEST( UnitCoarsener, verify_Eigen_coarsening )
         }
     }
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( std::move( data ) );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( std::move( data ) );
 
     std::vector<double> buffer;
     result->ReservePulls( result->Height() * result->Width() );
@@ -253,10 +253,10 @@ TEST( UnitCoarsener, verify_Eigen_coarsening )
 
 TEST( UnitCoarsener, 2d_grid )
 {
-    std::unique_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
-    std::unique_ptr<El::DistMatrix<double> > positions( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > positionsInterpolation( new El::DistMatrix<double>() );
-    std::unique_ptr<El::DistMatrix<double> > data( new El::DistMatrix<double>() );
+    std::shared_ptr<RBFFunctionInterface> rbfFunction( new TPSFunction() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positions( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > positionsInterpolation( new El::DistMatrix<double, El::VR, El::STAR>() );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > data( new El::DistMatrix<double, El::VR, El::STAR>() );
     El::Zeros( *positions, 20, 2 );
     El::Zeros( *positionsInterpolation, 50, 2 );
     El::Zeros( *data, positions->Height(), 1 );
@@ -366,9 +366,9 @@ TEST( UnitCoarsener, 2d_grid )
 
     UnitCoarsening rbf( 1e-3, 2, 15 );
 
-    rbf.compute( std::move( rbfFunction ), std::move( positions ), std::move( positionsInterpolation ) );
+    rbf.compute( rbfFunction, std::move( positions ), std::move( positionsInterpolation ) );
 
-    std::unique_ptr<El::DistMatrix<double> > result = rbf.interpolate( data );
+    std::unique_ptr<El::DistMatrix<double, El::VR, El::STAR> > result = rbf.interpolate( data );
 
     bufferPos.clear();
     data->ReservePulls( data->Height() * data->Width() );