diff --git a/Makefile b/Makefile
index 3673340..4746013 100644
--- a/Makefile
+++ b/Makefile
@@ -28,7 +28,7 @@ FLIBS = -lstdc++
 
 ## Includes and dependencies
 INCLUDE = -I $(INC_DIR)/
-DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o
+DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o $(OBJ_DIR)/SM_Helpers.o
 DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
 
 ## Directory structure
@@ -66,6 +66,9 @@ $(BIN_DIR) $(OBJ_DIR):
 $(OBJ_DIR)/%.o: $(TST_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
 	$(CXX) $(CXXFLAGS) $(INCLUDE) -c -o $@ $<
 
+$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
+	$(CXX) $(CXXFLAGS) -fPIE $(INCLUDE) -c -o $@ $<
+
 ## HDF5/C++ objects
 $(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR)
 	$(H5CXX) $(H5CXXFLAGS) $(INCLUDE) -c -o $@ $<
diff --git a/Makefile.verificarlo b/Makefile.verificarlo
index 93ccf0e..e9f135c 100644
--- a/Makefile.verificarlo
+++ b/Makefile.verificarlo
@@ -9,7 +9,7 @@ CXXFLAGS = -O0 -g $(H5FLAGS)
 FFLAGS = -O0 -g $(H5FLAGS)
 
 INCLUDE = -I $(INC_DIR)/
-DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o   $(OBJ_DIR)/SM_Standard.o  
+DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o $(OBJ_DIR)/SM_Helpers.o
 DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
 FLIBS = -lstdc++
 
@@ -49,6 +49,9 @@ $(BIN_DIR) $(OBJ_DIR):
 $(OBJ_DIR)/%.o: $(TST_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
 	$(CXX) $(CXXFLAGS) $(ARCH) $(INCLUDE) -c -o $@ $<
 
+$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp $(INC_DIR)/* | $(OBJ_DIR)
+	$(CXX) $(CXXFLAGS) -fPIE $(INCLUDE) -c -o $@ $<
+
 ## HDF5/C++ objects
 $(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR)
 	$(H5CXX) $(H5CXXFLAGS) $(INCLUDE) -c -o $@ $<
diff --git a/include/Helpers.hpp b/include/SM_Helpers.hpp
similarity index 72%
rename from include/Helpers.hpp
rename to include/SM_Helpers.hpp
index 06ea441..049df9c 100644
--- a/include/Helpers.hpp
+++ b/include/SM_Helpers.hpp
@@ -1,43 +1,48 @@
-// Helpers.hpp
+// SM_Helpers.hpp
 // Some usefull helper functions to support the Maponi algorithm.
 #include <iostream>
-#include <cmath>
 #include <string>
 #include <cstring>
-using namespace std;
+#include <cmath>
+
+void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
+
+void selectLargestDenominator(unsigned int l, unsigned int N_updates, 
+                      unsigned int *Updates_index, unsigned int *p, 
+                      double ***ylk);
 
 template<typename T>
-void showScalar(T scalar, string name) {
-    cout << name << " = " << scalar << endl << endl;
+void showScalar(T scalar, std::string name) {
+    std::cout << name << " = " << scalar << std::endl << std::endl;
 }
 
 template<typename T>
-void showVector(T *vector, unsigned int size, string name) {
-    cout << name << " = " << endl;
+void showVector(T *vector, unsigned int size, std::string name) {
+    std::cout << name << " = " << std::endl;
     for (unsigned int i = 0; i < size; i++) {
-        cout << "[ " << vector[i] << " ]" << endl;
+        std::cout << "[ " << vector[i] << " ]" << std::endl;
     }
-    cout << endl;
+    std::cout << std::endl;
 }
 
 template<typename T>
-void showMatrix(T *matrix, unsigned int M, string name) {
-    cout.precision(17);
-    cout << name << " = [" << endl;
+void showMatrix(T *matrix, unsigned int M, std::string name) {
+    std::cout.precision(17);
+    std::cout << name << " = [" << std::endl;
     for (unsigned int i = 0; i < M; i++) {
-        cout << "[";
+        std::cout << "[";
         for (unsigned int j = 0; j < M; j++) {
             if (matrix[i*M + j] >= 0) {
-                cout << "  " << matrix[i*M + j] << ",";
+                std::cout << "  " << matrix[i*M + j] << ",";
             }
             else {
-                cout << " " << matrix[i*M + j] << "," ;
+                std::cout << " " << matrix[i*M + j] << "," ;
             }
         }
-        cout << " ]," << endl;
+        std::cout << " ]," << std::endl;
     }
-    cout << "]" << endl;
-    cout << endl;
+    std::cout << "]" << std::endl;
+    std::cout << std::endl;
 }
 
 template<typename T>
diff --git a/smvars.sh b/smvars.sh
index b9695aa..5807a32 100644
--- a/smvars.sh
+++ b/smvars.sh
@@ -43,7 +43,7 @@ case $ENV in
 	;;
   *)
     echo "Unknown environment descriptor given."
-	echo "Usage: source smvars.sh {intel | llvm | gnu}"
+	echo "Usage: source smvars.sh {intel | llvm | vfc | gnu}"
 	return 1
 	;;
 esac
diff --git a/src/SM_Helpers.cpp b/src/SM_Helpers.cpp
new file mode 100644
index 0000000..185c367
--- /dev/null
+++ b/src/SM_Helpers.cpp
@@ -0,0 +1,31 @@
+#include "SM_Helpers.hpp"
+
+void Switch(unsigned int *p, unsigned int l, unsigned int lbar) {
+  unsigned int tmp = p[l+1];
+  p[l+1] = p[lbar];
+  p[lbar] = tmp;
+}
+
+void selectLargestDenominator(unsigned int l, unsigned int N_updates, 
+                      unsigned int *Updates_index, unsigned int *p, 
+                      double ***ylk) {
+  unsigned int lbar = l+1, max =0;
+  unsigned int index = 0, component = 0;
+  unsigned int tmp = 0;
+  double breakdown = 0;
+  for (unsigned int j = lbar; j < N_updates + 1; j++) {
+    index = p[j];
+    component = Updates_index[index - 1];
+    breakdown = abs(1 + ylk[l][index][component]);
+    #ifdef DEBUG
+    std::cout << "Inside selectLargestDenominator()" << std::endl;
+    std::cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << std::endl;
+    std::cout << std::endl;
+    #endif
+    if (breakdown > max) {
+      max = breakdown;
+      lbar = j;
+    }
+  }
+  Switch(p, l, lbar);
+}
diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp
index df4540f..1196ed5 100644
--- a/src/SM_MaponiA3.cpp
+++ b/src/SM_MaponiA3.cpp
@@ -2,33 +2,9 @@
 // Algorithm 3 from P. Maponi,
 // p. 283, doi:10.1016/j.laa.2006.07.007
 #include "SM_MaponiA3.hpp"
-#include "Helpers.hpp"
+#include "SM_Helpers.hpp"
 
-void selectBestUpdate(unsigned int l, unsigned int N_updates, 
-                      unsigned int *Updates_index, unsigned int *p, 
-                      double ***ylk) {
-  unsigned int lbar = l+1, max =0;
-  unsigned int index = 0, component = 0;
-  unsigned int tmp = 0;
-  double breakdown = 0;
-  for (unsigned int j = lbar; j < N_updates + 1; j++) {
-    index = p[j];
-    component = Updates_index[index - 1];
-    breakdown = abs(1 + ylk[l][index][component]);
-    #ifdef DEBUG
-    cout << "Inside selectBestUpdate()" << endl;
-    cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << endl;
-    cout << endl;
-    #endif
-    if (breakdown > max) {
-      max = breakdown;
-      lbar = j;
-    }
-  }
-  tmp = p[l+1];
-  p[l+1] = p[lbar];
-  p[lbar] = tmp; 
-}
+// #define DEBUG
 
 void MaponiA3(double *Slater_inv, unsigned int Dim,
               unsigned int N_updates, double *Updates,
@@ -68,9 +44,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
   // Calculate the {y_{0,k}}
   for (k = 1; k < N_updates + 1; k++) {
     #ifdef DEBUG
-    cout << "Compute y0k: " << endl;
-    cout << "ylk[0][" << k << "][:]";
-    cout << endl;
+    std::cout << "Compute y0k: " << std::endl;
+    std::cout << "ylk[0][" << k << "][:]";
+    std::cout << std::endl;
     #endif
     for (i = 1; i < Dim + 1; i++) {
       for (j = 1; j < Dim + 1; j++) {
@@ -86,33 +62,33 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
   // Calculate the {y_{l,k}} from the {y_{0,k}}
   for (l = 0; l < N_updates; l++) {
     #ifdef DEBUG
-    cout << "In outer compute-ylk-loop: l = " << l << endl;
-    cout << endl;
+    std::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
+    std::cout << std::endl;
     #endif
 
     // For given l select intermediate update with largest break-down val
-    selectBestUpdate(l, N_updates, Updates_index, p, ylk);
+    selectLargestDenominator(l, N_updates, Updates_index, p, ylk);
 
     // Select component and comp. bd-condition.
     component = Updates_index[p[l+1] - 1];
     beta = 1 + ylk[l][p[l+1]][component];
     #ifdef DEBUG
-    cout << "p[l+1] = " << p[l+1] << endl;
-    cout << "component = " << component << endl;
-    cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl;
-    cout << endl;
+    std::cout << "p[l+1] = " << p[l+1] << std::endl;
+    std::cout << "component = " << component << std::endl;
+    std::cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << std::endl;
+    std::cout << std::endl;
     #endif
     if (fabs(beta) < 1e-3) {
-      cerr << "Break-down occured." << endl;
+      std::cerr << "Break-down occured." << std::endl;
     }
 
     // Compute intermediate update to Slater_inv
     #ifdef DEBUG
-    cout << "Compute intermediate update to Slater_inv" << endl;
-    cout << "component = " << component << endl;
-    cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl;
-    cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << endl;
-    cout << endl;
+    std::cout << "Compute intermediate update to Slater_inv" << std::endl;
+    std::cout << "component = " << component << std::endl;
+    std::cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << std::endl;
+    std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << std::endl;
+    std::cout << std::endl;
     #endif
     for (i = 0; i < Dim; i++) {
       for (j = 0; j < Dim; j++) {
@@ -132,9 +108,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     for (k = l+2; k < N_updates+1; k++) {
       alpha = ylk[l][p[k]][component] / beta;
       #ifdef DEBUG
-      cout << "Inside k-loop: k = " << k << endl;
-      cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << endl;
-      cout << endl;
+      std::cout << "Inside k-loop: k = " << k << std::endl;
+      std::cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << std::endl;
+      std::cout << std::endl;
       #endif
       for (i = 1; i < Dim + 1; i++) {
         ylk[l+1][p[k]][i] = ylk[l][p[k]][i]
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index 8ac8b72..2354450 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -1,8 +1,7 @@
 // SM-Standard.cpp
 // Standard Sherman Morrison with multiple updates
 #include "SM_Standard.hpp"
-#include "Helpers.hpp"
-#include <iostream>
+#include "SM_Helpers.hpp"
 
 // Set common break-down threshold
 double threshold() {
diff --git a/tests/cMaponiA3_test_3x3_3.cpp b/tests/cMaponiA3_test_3x3_3.cpp
index 525d15d..20f7d45 100644
--- a/tests/cMaponiA3_test_3x3_3.cpp
+++ b/tests/cMaponiA3_test_3x3_3.cpp
@@ -1,6 +1,6 @@
 // main.cpp
 #include "SM_MaponiA3.hpp"
-#include "Helpers.hpp"
+#include "SM_Helpers.hpp"
 
 int main() {
 
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index eb719e5..12d9b38 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -1,11 +1,9 @@
-#include <iostream>
-#include <string>
 #include "hdf5/serial/hdf5.h"
 #include "hdf5/serial/H5Cpp.h"
 
 #include "SM_MaponiA3.hpp"
 #include "SM_Standard.hpp"
-#include "Helpers.hpp"
+#include "SM_Helpers.hpp"
 
 using namespace H5;
 // #define DEBUG