Merge pull request #34 from fmgjcoppens/fix/slagel_splitting

Various fixes

Makefile

@@ -29,13 +29,12 @@ FLIBS = -lstdc++
 
 ## Includes and dependencies
 INCLUDE = -I $(INC_DIR)/
-DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
+DEPS_CXX = $(OBJ_DIR)/SM_Maponi.o \
-		   $(OBJ_DIR)/SM_MaponiA3S.o \
 		   $(OBJ_DIR)/SM_Standard.o \
 		   $(OBJ_DIR)/SM_Helpers.o
 DEPS_F = $(DEPS_CXX) \
-		 $(OBJ_DIR)/SM_mod.o \
+		 $(OBJ_DIR)/finterface_mod.o \
-		 $(OBJ_DIR)/Helpers_mod.o
+		 $(OBJ_DIR)/helpers_mod.o
 
 ## Directory structure
 SRC_DIR := src

Makefile.verificarlo

@@ -16,13 +16,12 @@ CXXFLAGS = -O0 -g $(H5FLAGS)
 FFLAGS = $(CXXFLAGS)
 
 INCLUDE = -I $(INC_DIR)/
-DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
+DEPS_CXX = $(OBJ_DIR)/SM_Maponi.o \
-		   $(OBJ_DIR)/SM_MaponiA3S.o \
 		   $(OBJ_DIR)/SM_Standard.o \
 		   $(OBJ_DIR)/SM_Helpers.o
 DEPS_F = $(DEPS_CXX) \
-		 $(OBJ_DIR)/SM_mod.o \
+		 $(OBJ_DIR)/finterface_mod.o \
-		 $(OBJ_DIR)/Helpers_mod.o
+		 $(OBJ_DIR)/helpers_mod.o
 FLIBS = -lstdc++
 
 SRC_DIR := src

include/SM_Helpers.hpp

@@ -127,3 +127,53 @@ template <typename T> bool is_identity(T *A, unsigned int M, double tolerance) {
   }
   return true;
 }
+
+template <typename T> T norm_max(T * A, unsigned int Dim) {
+  T res = 0;
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      T delta = A[i * Dim + j];
+      delta = fabs(delta);
+      if (delta > res) {
+        res = delta;
+      }
+    }
+  }
+  return res;
+}
+
+template  <typename T> T norm_frobenius2(T * A, unsigned int Dim) {
+  T res = 0;
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      T delta = A[i * Dim + j];
+      res += delta*delta;
+    }
+  }
+  return res;
+}
+
+template <typename T> T residual_max(T * A, unsigned int Dim) {
+  T res = 0;
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      T  delta = A[i * Dim + j] - (i == j);
+      delta = fabs(delta);
+      if (delta > res) {
+        res = delta;
+      }
+    }
+  }
+  return res;
+}
+
+template <typename T> T residual_frobenius2(T * A, unsigned int Dim) {
+  T res = 0;
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      T delta = A[i * Dim + j] - (i == j);
+      res += delta*delta;
+    }
+  }
+  return res;
+}

include/SM_Maponi.hpp

@@ -0,0 +1,8 @@
+// P. Maponi Algorithm 3
+void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+              double *Updates, unsigned int *Updates_index);
+
+// P. Maponi Algorithm 3 with J. Slagel splitting
+// http://hdl.handle.net/10919/52966
+void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+               double *Updates, unsigned int *Updates_index);

include/SM_MaponiA3.hpp

@@ -1,2 +0,0 @@
-void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-              double *Updates, unsigned int *Updates_index);

include/SM_MaponiA3S.hpp

@@ -1,2 +0,0 @@
-void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-               double *Updates, unsigned int *Updates_index);

smvars.sh

@@ -45,6 +45,7 @@ case $ENV in
   *)
     echo "Unknown environment descriptor given."
 	echo "Usage: source smvars.sh {intel | llvm | vfc | gnu}"
+	echo "Usage: source smvars.sh {intel | llvm | vfc | gnu} [threshold]"
 	return 1
 	;;
 esac

src/SM_Maponi.cpp

@@ -1,17 +1,146 @@
-// SM-MaponiA3_f.cpp
+// SM_Maponi.cpp
 // Algorithm 3 from P. Maponi,
 // p. 283, doi:10.1016/j.laa.2006.07.007
-#include "SM_MaponiA3S.hpp"
+#include "SM_Maponi.hpp"
 #include "SM_Helpers.hpp"
 
 // #define DEBUG
 
-void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-               double *Updates, unsigned int *Updates_index) {
+              double *Updates, unsigned int *Updates_index) {
+  std::cerr << "Called MaponiA3 with " << N_updates << " updates" << std::endl;
+
+  /*
+  DECLARE AND INITIALISE ARRAYS
+  */
+  unsigned int k, l, i, j, component;
+  unsigned int *p = new unsigned int[N_updates + 1]{0};
+  double alpha, beta;
+  double *Al = new double[Dim * Dim];
+  double *next = new double[Dim * Dim]{0};
+  double *last = Slater_inv, *tmp;
+
+  // Populate update-order vector
+  for (i = 0; i < N_updates; i++) {
+    p[i + 1] = i + 1;
+  }
+
+  // Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim]
+  double ***ylk = new double **[N_updates];
+  for (l = 0; l < N_updates; l++) {
+    ylk[l] = new double *[N_updates + 1];
+    for (k = 0; k < N_updates + 1; k++) {
+      ylk[l][k] = new double[Dim + 1]{0};
+    }
+  }
+
+  /*
+  START ALGORITHM
+  */
+
+  // Calculate the {y_{0,k}}
+  for (k = 1; k < N_updates + 1; k++) {
+#ifdef DEBUG
+    std::cerr << "Compute y0k: " << std::endl;
+    std::cerr << "ylk[0][" << k << "][:]";
+    std::cerr << std::endl;
+#endif
+    for (i = 1; i < Dim + 1; i++) {
+      for (j = 1; j < Dim + 1; j++) {
+        ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] *
+                        Updates[(k - 1) * Dim + (j - 1)];
+      }
+    }
+#ifdef DEBUG
+    showVector(ylk[0][k], Dim, "");
+#endif
+  }
+
+  // Calculate the {y_{l,k}} from the {y_{0,k}}
+  for (l = 0; l < N_updates; l++) {
+#ifdef DEBUG
+    std::cerr << "In outer compute-ylk-loop: l = " << l << std::endl;
+    std::cerr << std::endl;
+#endif
+
+    // For given l select intermediate update with largest break-down val
+    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
+    std::cerr << "p[l+1] = " << p[l + 1] << std::endl;
+    std::cerr << "component = " << component << std::endl;
+    std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
+              << "] = " << beta << std::endl;
+    std::cerr << std::endl;
+#endif
+    if (fabs(beta) < threshold()) {
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l]
+                << std::endl;
+    }
+    double ibeta = 1.0 / beta; 
+
+// Compute intermediate update to Slater_inv
+#ifdef DEBUG
+    std::cerr << "Compute intermediate update to Slater_inv" << std::endl;
+    std::cerr << "component = " << component << std::endl;
+    std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
+              << "]" << std::endl;
+    std::cerr << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
+              << std::endl;
+    std::cerr << std::endl;
+#endif
+    for (i = 0; i < Dim; i++) {
+      for (j = 0; j < Dim; j++) {
+        Al[i * Dim + j] =
+            (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
+      }
+    }
+    matMul(Al, last, next, Dim);
+    tmp = next;
+    next = last;
+    last = tmp;
+#ifdef DEBUG
+    showMatrix(last, Dim, "last");
+#endif
+
+    // For given l != 0 compute the next {y_{l,k}}
+    for (k = l + 2; k < N_updates + 1; k++) {
+      alpha = ylk[l][p[k]][component] * ibeta;
+#ifdef DEBUG
+      std::cerr << "Inside k-loop: k = " << k << std::endl;
+      std::cerr << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;
+      std::cerr << std::endl;
+#endif
+      for (i = 1; i < Dim + 1; i++) {
+        ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i];
+      }
+    }
+  }
+  memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
+
+  /*
+  CLEANUP MEMORY
+  */
+
+  for (l = 0; l < N_updates; l++) {
+    for (k = 0; k < N_updates + 1; k++) {
+      delete[] ylk[l][k];
+    }
+    delete[] ylk[l];
+  }
+  delete[] Al, next, p, ylk;
+}
+
+void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+               double *Updates, unsigned int *Updates_index) {
+  std::cerr << "Called MaponiA3S with " << N_updates << " updates" << std::endl;
+
   /*
   DECLARE AND INITIALISE ARRAYS
   */
-  showVector(Updates_index, Dim, "Updates_index");
   unsigned int k, l, i, j, component;
   unsigned int *p = new unsigned int[N_updates + 1]{0};
   double alpha, beta;
@@ -40,7 +169,6 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   /*
   START ALGORITHM
   */
- showMatrix(Slater_inv, Dim, "S^{-1}");
 
   // Calculate the y_{0,k}  = S_0^{-1} u_k
   for (k = 1; k < N_updates + 1; k++) {
@@ -81,27 +209,17 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     std::cerr << std::endl;
 #endif
     if (fabs(beta) < threshold()) {
-#ifdef DEBUG
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l]
-      std::cerr << "Breakdown condition triggered at l = " << l << " and component "<< component
                 << std::endl;
-#endif
       for (unsigned int i = 1; i < Dim + 1; i++) {
-        std::cout << "Let's split..." << std::endl;
         ylk[l][p[l + 1]][i] *= 0.5;
-        // later_updates[later * Dim + i - 1] =  ylk[l][p[l + 1]][i]; // Updates[l * Dim + i] * 0.5;
         later_updates[later * Dim + i - 1] = Updates[l * Dim + i - 1] * 0.5;
       }
       later_index[later] = Updates_index[p[l + 1] - 1];
-      std::cout << "Later = " << later << std::endl;
       later++;
-      std::cout << "Later = " << later << std::endl;
       beta = 1 + ylk[l][p[l + 1]][component];
-      std::cout << "New beta is: " << beta << std::endl;
     }
     double ibeta = 1.0 / beta;
-    showVector(ylk[l][p[l + 1]], Dim+1, "");
-    showVector(later_updates+(Dim*0), Dim, "later_updateds");
-    
 
 // Compute intermediate update to Slater_inv
 #ifdef DEBUG
@@ -135,7 +253,6 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
       for (i = 1; i < Dim + 1; i++) {
         ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i];
       }
-      showVector(ylk[l + 1][p[k]], Dim+1, "");
     }
   }
   memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
@@ -153,12 +270,18 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   delete[] Al, next, p, ylk;
 
   if (later > 0) {
-    std::cerr << "Entering recursive loop with " << later << " updates" << std::endl;
-    std::cout << later_index[0] << std::endl;
     MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
   }
 }
 
+extern "C" {
+  void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
+}
+
 extern "C" {
   void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
                   unsigned int *N_updates, double **linUpdates,

src/SM_MaponiA3.cpp

@@ -1,142 +0,0 @@
-// SM-MaponiA3_f.cpp
-// Algorithm 3 from P. Maponi,
-// p. 283, doi:10.1016/j.laa.2006.07.007
-#include "SM_MaponiA3.hpp"
-#include "SM_Helpers.hpp"
-
-// #define DEBUG
-
-void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-              double *Updates, unsigned int *Updates_index) {
-
-  /*
-  DECLARE AND INITIALISE ARRAYS
-  */
-
-  unsigned int k, l, i, j, component;
-  unsigned int *p = new unsigned int[N_updates + 1]{0};
-  double alpha, beta;
-  double *Al = new double[Dim * Dim];
-  double *next = new double[Dim * Dim]{0};
-  double *last = Slater_inv, *tmp;
-
-  // Populate update-order vector
-  for (i = 0; i < N_updates; i++) {
-    p[i + 1] = i + 1;
-  }
-
-  // Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim]
-  double ***ylk = new double **[N_updates];
-  for (l = 0; l < N_updates; l++) {
-    ylk[l] = new double *[N_updates + 1];
-    for (k = 0; k < N_updates + 1; k++) {
-      ylk[l][k] = new double[Dim + 1]{0};
-    }
-  }
-
-  /*
-  START ALGORITHM
-  */
-
-  // Calculate the {y_{0,k}}
-  for (k = 1; k < N_updates + 1; k++) {
-#ifdef DEBUG
-    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++) {
-        ylk[0][k][i] += Slater_inv[(i - 1) * Dim + (j - 1)] *
-                        Updates[(k - 1) * Dim + (j - 1)];
-      }
-    }
-#ifdef DEBUG
-    showVector(ylk[0][k], Dim, "");
-#endif
-  }
-
-  // Calculate the {y_{l,k}} from the {y_{0,k}}
-  for (l = 0; l < N_updates; l++) {
-#ifdef DEBUG
-    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
-    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
-    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) < threshold()) {
-      std::cerr << "Break-down occured." << std::endl;
-    }
-    double ibeta = 1.0 / beta; 
-
-// Compute intermediate update to Slater_inv
-#ifdef DEBUG
-    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++) {
-        Al[i * Dim + j] =
-            (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
-      }
-    }
-    matMul(Al, last, next, Dim);
-    tmp = next;
-    next = last;
-    last = tmp;
-#ifdef DEBUG
-    showMatrix(last, Dim, "last");
-#endif
-
-    // For given l != 0 compute the next {y_{l,k}}
-    for (k = l + 2; k < N_updates + 1; k++) {
-      alpha = ylk[l][p[k]][component] * ibeta;
-#ifdef DEBUG
-      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] - alpha * ylk[l][p[l + 1]][i];
-      }
-    }
-  }
-  memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
-
-  /*
-  CLEANUP MEMORY
-  */
-
-  for (l = 0; l < N_updates; l++) {
-    for (k = 0; k < N_updates + 1; k++) {
-      delete[] ylk[l][k];
-    }
-    delete[] ylk[l];
-  }
-  delete[] Al, next, p, ylk;
-}
-
-extern "C" {
-  void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
-                  unsigned int *N_updates, double **linUpdates,
-                  unsigned int **Updates_index) {
-    MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
-  }
-}

src/SM_Standard.cpp

@@ -50,7 +50,7 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // http://hdl.handle.net/10919/52966
 void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
-  std::cerr << "Called SM2 with updates " << N_updates << std::endl;
+  std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl;
   double C[Dim];
   double D[Dim];
 
@@ -110,7 +110,7 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // Sherman Morrison, leaving zero denominators for later
 void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
-  std::cerr << "Called SM3 with updates " << N_updates << std::endl;
+  std::cerr << "Called SM3 with " << N_updates << " updates" << std::endl;
   double C[Dim];
   double D[Dim];
 
@@ -176,7 +176,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // Leave zero denominators for later (SM3), and when none are left then split (SM2)
 void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
-  std::cerr << "Called SM4 with updates " << N_updates << std::endl;
+  std::cerr << "Called SM4 with " << N_updates << " updates" << std::endl;
   double C[Dim];
   double D[Dim];
 
@@ -254,6 +254,7 @@ extern "C" {
                   unsigned int **Updates_index) {
     SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
   }
+
   void SM4_f(double **linSlater_inv, unsigned int *Dim,
                   unsigned int *N_updates, double **linUpdates,
                   unsigned int **Updates_index) {

tests/cMaponiA3_test_3x3_3.cpp

@@ -1,5 +1,5 @@
 // main.cpp
-#include "SM_MaponiA3.hpp"
+#include "SM_Maponi.hpp"
 #include "SM_Helpers.hpp"
 
 int main() {

tests/run_test.sh

@@ -12,7 +12,8 @@ ALGO=$1
 START=$2
 STOP=$3
 BLOCKSIZE=329
-TOLERANCE=1e-3
+#TOLERANCE=1e-3
+TOLERANCE=$4
 
 BLOCKS=$((STOP / BLOCKSIZE))
 LEFT=$((STOP % BLOCKSIZE))

tests/test_h5.cpp

@@ -1,39 +1,15 @@
 #include "hdf5/serial/hdf5.h"
 #include "hdf5/serial/H5Cpp.h"
 
-#include "SM_MaponiA3.hpp"
+#include "SM_Maponi.hpp"
-#include "SM_MaponiA3S.hpp"
 #include "SM_Standard.hpp"
 #include "SM_Helpers.hpp"
 
 using namespace H5;
-#define DEBUG
+// #define DEBUG
 
 const H5std_string FILE_NAME( "dataset.hdf5" );
 
-double residual_max(double * A, unsigned int Dim) {
-  double max = 0.0;
-  for (unsigned int i = 0; i < Dim; i++) {
-    for (unsigned int j = 0; j < Dim; j++) {
-      double delta = (A[i * Dim + j] - (i == j));
-      delta = abs(delta);
-      if (delta > max) max = delta;
-    }
-  }
-  return max;
-}
-
-double residual2(double * A, unsigned int Dim) {
-  double res = 0.0;
-  for (unsigned int i = 0; i < Dim; i++) {
-    for (unsigned int j = 0; j < Dim; j++) {
-      double delta = (A[i * Dim + j] - (i == j));
-      res += delta*delta;
-    }
-  }
-  return res;
-}
-
 void read_int(H5File file, std::string key, unsigned int * data) {
   DataSet ds = file.openDataSet(key);
   ds.read(data, PredType::STD_U32LE);
@@ -92,8 +68,6 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   showMatrix(u, dim, "Updates");
 #endif
 
-
-
   if (version == "maponia3") {
     MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "maponia3s") {
@@ -124,7 +98,7 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   bool ok = is_identity(res, dim, tolerance);
 
   double res_max = residual_max(res, dim);
-  double res2 = residual2(res, dim);
+  double res2 = residual_frobenius2(res, dim);
   std::cout << "Residual = " << version << " " << cycle << " " << res_max << " " << res2 << std::endl;
 
 #ifdef DEBUG

tests/vfc_test_h5.cpp

@@ -11,7 +11,7 @@
 #include <sstream>
 
 
-#include "SM_MaponiA3.hpp"
+#include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
 #include "SM_Helpers.hpp"
 #include "vfc_probe.h"
@@ -21,29 +21,6 @@ using namespace H5;
 
 const H5std_string FILE_NAME( "datasets/ci_dataset.hdf5" );
 
-double residual_max(double * A, unsigned int Dim) {
-  double max = 0.0;
-  for (unsigned int i = 0; i < Dim; i++) {
-    for (unsigned int j = 0; j < Dim; j++) {
-      double delta = (A[i * Dim + j] - (i == j));
-      delta = abs(delta);
-      if (delta > max) max = delta;
-    }
-  }
-  return max;
-}
-
-double residual2(double * A, unsigned int Dim) {
-  double res = 0.0;
-  for (unsigned int i = 0; i < Dim; i++) {
-    for (unsigned int j = 0; j < Dim; j++) {
-      double delta = (A[i * Dim + j] - (i == j));
-      res += delta*delta;
-    }
-  }
-  return res;
-}
-
 void read_int(H5File file, std::string key, unsigned int * data) {
   DataSet ds = file.openDataSet(key);
   ds.read(data, PredType::STD_U32LE);
@@ -56,7 +33,6 @@ void read_double(H5File file, std::string key, double * data) {
   ds.close();
 }
 
-
 /* Return a vector containing all cycles to execute by reading a data file */
 std::vector<int> get_cycles_list(std::string path) {
   std::ifstream file_stream(path);
@@ -157,7 +133,7 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
   bool ok = is_identity(res, dim, 1e-3);
 
   double res_max = residual_max(res, dim);
-  double res2 = residual2(res, dim);
+  double res2 = residual_frobenius2(res, dim);
 
 #ifdef DEBUG
   showMatrix(res, dim, "Result");