diff --git a/src/heat/README.md b/src/heat/README.md
new file mode 100644
index 0000000..0d642a0
--- /dev/null
+++ b/src/heat/README.md
@@ -0,0 +1,9 @@
+The heat benchmark solves the steady heat equation in a regular grid of NxM
+elements using an iterative solver.
+
+The solver is either the Gauss-Seidel or Successive-over-relaxation with a wiven
+relaxation parameter (--relax).
+
+In every iteration the relative error of the solution is computed by using the
+infinite norm until the tolerance limit is reached.
+
diff --git a/src/heat/results b/src/heat/results
new file mode 100644
index 0000000..019865e
--- /dev/null
+++ b/src/heat/results
@@ -0,0 +1,8 @@
+In MareNostrum 4:
+24cores  -s_16384_-t_500_-b_1024   3906.56
+24cores  -s_16384_-t_500_-b_2048   3890.20
+24cores  -s_16384_-t_500_-b_128    3752.77
+
+48cores  -s_16384_-t_500_-b_1024   4843.91
+48cores  -s_16384_-t_500_-b_128    4367.97
+48cores  -s_16384_-t_500_-b_2048   4275.91
diff --git a/src/heat/smp/solver_ompss2_residual.c b/src/heat/smp/solver_ompss2_residual.c
index 03baa7f..17d59b3 100644
--- a/src/heat/smp/solver_ompss2_residual.c
+++ b/src/heat/smp/solver_ompss2_residual.c
@@ -1,4 +1,5 @@
 #include <string.h>
+#include <stdio.h>
 
 #include "common/heat.h"
 
@@ -9,41 +10,74 @@ summary(void)
 		"residual";
 }
 
-static inline void gaussSeidelSolver(int64_t rows, int64_t cols, int rbs, int cbs, int nrb, int ncb, double M[rows][cols], char reps[nrb][ncb], double *residual)
+static inline void gaussSeidelSolver(int64_t rows, int64_t cols, int rbs, int cbs, int nrb, int ncb, double M[rows][cols], char reps[nrb][ncb], double *residual, double *max_elem, double relax)
 {
 	for (int R = 1; R < nrb-1; ++R) {
 		for (int C = 1; C < ncb-1; ++C) {
 			#pragma oss task label("block computation") \
 					in(reps[R-1][C]) in(reps[R+1][C]) \
 					in(reps[R][C-1]) in(reps[R][C+1]) \
-					inout(reps[R][C]) reduction(+: [1]residual)
-			*residual += computeBlockResidual(rows, cols, (R-1)*rbs+1, R*rbs, (C-1)*cbs+1, C*cbs, M);
+					inout(reps[R][C]) \
+					reduction(max: [1]residual) \
+					reduction(max: [1]max_elem)
+			{
+				double lresidual = 0.0;
+				double lmax_elem = 0.0;
+
+				computeBlockResidual(rows, cols, (R-1)*rbs+1,
+						R*rbs, (C-1)*cbs+1, C*cbs, M,
+						relax, &lresidual, &lmax_elem);
+
+				*residual = fmax(*residual, lresidual);
+				*max_elem = fmax(*max_elem, lmax_elem);
+			}
 		}
 	}
 }
 
 double solve(HeatConfiguration *conf, int64_t rows, int64_t cols, int timesteps, void *extraData)
 {
+	FILE *f = fopen("convergence.csv", "w");
+	fprintf(f, "iter error time\n");
+
 	(void) extraData;
 	double (*matrix)[cols] = (double (*)[cols]) conf->matrix;
 	const double delta = conf->delta;
 	const int rbs = conf->rbs;
 	const int cbs = conf->cbs;
 
-	const int N = 10;
+	const int N = 4;
 	double results[N];
-	for (int i = 0; i < N; ++i)
-		results[i] = delta;
+	double max_elem[N];
+	double residual[N];
+
+	for (int i = 0; i < N; ++i) {
+		results[i] = 666;
+		max_elem[i] = 666;
+		residual[i] = 666;
+	}
 
 	const int nrb = (rows-2)/rbs+2;
 	const int ncb = (cols-2)/cbs+2;
 	char representatives[nrb][ncb];
 
+	double t0 = getTime();
+
 	int t = 0;
 	while (t < timesteps) {
 		results[t%N] = 0.0f;
+		max_elem[t%N] = 0.0f;
+		residual[t%N] = 0.0f;
 
-		gaussSeidelSolver(rows, cols, rbs, cbs, nrb, ncb, matrix, representatives, &results[t%N]);
+		gaussSeidelSolver(rows, cols, rbs, cbs, nrb, ncb, matrix,
+				representatives, &residual[t%N], &max_elem[t%N], conf->relax);
+
+		#pragma oss task in(residual[t%N], max_elem[t%N]) out(results[t%N])
+		{
+			results[t%N] = residual[t%N] / max_elem[t%N];
+			fprintf(f, "%d %e %e\n", t, results[t%N], getTime() - t0);
+			//fprintf(stderr, "t=%d error=%e\n", t, results[t%N]);
+		}
 
 		// Advance to the next timestep
 		++t;
@@ -60,5 +94,7 @@ double solve(HeatConfiguration *conf, int64_t rows, int64_t cols, int timesteps,
 	// Save the number of performed timesteps
 	conf->convergenceTimesteps = t;
 
+	fclose(f);
+
 	return results[(t-1)%N];
 }