diff --git a/mex/sources/block_kalman_filter/block_kalman_filter.cc b/mex/sources/block_kalman_filter/block_kalman_filter.cc
index c179a5485..e6a82d673 100644
--- a/mex/sources/block_kalman_filter/block_kalman_filter.cc
+++ b/mex/sources/block_kalman_filter/block_kalman_filter.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2017 Dynare Team
+ * Copyright © 2007-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -16,16 +16,14 @@
  * You should have received a copy of the GNU General Public License
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
-#include <cstring>
+
 #include <cmath>
-#include <iostream>
-#include <fstream>
-#include <vector>
+#include <algorithm>
 #ifdef USE_OMP
 # include <omp.h>
 #endif
 #include "block_kalman_filter.hh"
-using namespace std;
+
 #define BLAS
 //#define CUBLAS
 
@@ -33,17 +31,18 @@ using namespace std;
 # include <cuda_runtime.h>
 # include <cublas_v2.h>
 #endif
+
 void
-mexDisp(mxArray *P)
+mexDisp(const mxArray *P)
 {
-  unsigned int n = mxGetN(P);
-  unsigned int m = mxGetM(P);
-  double *M = mxGetPr(P);
+  size_t n = mxGetN(P);
+  size_t m = mxGetM(P);
+  const double *M = mxGetPr(P);
   mexPrintf("%d x %d\n", m, n);
   mexEvalString("drawnow;");
-  for (unsigned int i = 0; i < m; i++)
+  for (size_t i = 0; i < m; i++)
     {
-      for (unsigned int j = 0; j < n; j++)
+      for (size_t j = 0; j < n; j++)
         mexPrintf(" %9.4f", M[i+ j * m]);
       mexPrintf("\n");
     }
@@ -51,7 +50,7 @@ mexDisp(mxArray *P)
 }
 
 void
-mexDisp(double *M, int m, int n)
+mexDisp(const double *M, int m, int n)
 {
   mexPrintf("%d x %d\n", m, n);
   mexEvalString("drawnow;");
@@ -63,6 +62,7 @@ mexDisp(double *M, int m, int n)
     }
   mexEvalString("drawnow;");
 }
+
 /*if block
   %nz_state_var = M_.nz_state_var;
   while notsteady && t<smpl
@@ -118,13 +118,12 @@ mexDisp(double *M, int m, int n)
 */
 
 bool
-not_all_abs_F_bellow_crit(double *F, int size, double crit)
+not_all_abs_F_bellow_crit(const double *F, int size, double crit)
 {
   int i = 0;
   while (i < size && abs(F[i]) < crit)
-    {
-      i++;
-    }
+    i++;
+
   if (i < size)
     return false;
   else
@@ -132,7 +131,7 @@ not_all_abs_F_bellow_crit(double *F, int size, double crit)
 }
 
 double
-det(double *F, int dim, lapack_int *ipiv)
+det(const double *F, int dim, const lapack_int *ipiv)
 {
   double det = 1.0;
   for (int i = 0; i < dim; i++)
@@ -143,7 +142,7 @@ det(double *F, int dim, lapack_int *ipiv)
   return det;
 }
 
-BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[], double *P_mf[], double *v_pp[], double *K[], double *v_n[], double *a[], double *K_P[], double *P_t_t1[], double *tmp[], double *P[])
+BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
   if (nlhs > 3)
     DYN_MEX_FUNC_ERR_MSG_TXT("block_kalman_filter provides at most 3 output argument.");
@@ -203,7 +202,7 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
   R = mxGetPr(pR);
   Q = mxGetPr(pQ);
   H = mxGetPr(pH);
-  *P = mxGetPr(pP);
+  P = mxGetPr(pP);
   Y = mxGetPr(pY);
 
   n   = mxGetN(pT);           // Number of state variables.
@@ -222,13 +221,13 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
     if (mxGetNumberOfElements(pdata_index) != static_cast<unsigned int>(smpl))
       DYN_MEX_FUNC_ERR_MSG_TXT("the number of element in the cell array passed to block_missing_observation_kalman_filter as first argument has to be equal to the smpl size");
 
-  i_nz_state_var = static_cast<int *>(mxMalloc(n*sizeof(int)));
+  i_nz_state_var = std::make_unique<int[]>(n);
   for (int i = 0; i < n; i++)
     i_nz_state_var[i] = nz_state_var[i];
 
   pa = mxCreateDoubleMatrix(n, 1, mxREAL);         // State vector.
-  *a = mxGetPr(pa);
-  tmp_a = static_cast<double *>(mxMalloc(n * sizeof(double)));
+  a = mxGetPr(pa);
+  tmp_a = std::make_unique<double[]>(n);
   dF = 0.0;                                            // det(F).
 
   p_tmp1 = mxCreateDoubleMatrix(n, n_shocks, mxREAL);
@@ -240,12 +239,11 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
   LIK  = 0.0;                                          // Default value of the log likelihood.
   notsteady   = true;                                 // Steady state flag.
   F_singular  = true;
-  *v_pp = static_cast<double *>(mxMalloc(pp * sizeof(double)));
-  *v_n = static_cast<double *>(mxMalloc(n * sizeof(double)));
-  mf = static_cast<int *>(mxMalloc(pp * sizeof(int)));
+  v_pp = std::make_unique<double[]>(pp);
+  v_n = std::make_unique<double[]>(n);
+  mf = std::make_unique<int[]>(pp);
   for (int i = 0; i < pp; i++)
     mf[i] = mfd[i] - 1;
-  pi = atan2(0.0, -1.0);
 
   /*compute QQ = R*Q*transpose(R)*/                        // Variance of R times the vector of structural innovations.;
   // tmp = R * Q;
@@ -278,109 +276,104 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
   piF =  mxCreateDoubleMatrix(pp, pp, mxREAL);
   iF = mxGetPr(piF);
   lw = pp * 4;
-  w = static_cast<double *>(mxMalloc(lw * sizeof(double)));
-  iw = static_cast<lapack_int *>(mxMalloc(pp * sizeof(lapack_int)));
-  ipiv = static_cast<lapack_int *>(mxMalloc(pp * sizeof(lapack_int)));
+  w = std::make_unique<double[]>(lw);
+  iw = std::make_unique<lapack_int[]>(pp);
+  ipiv = std::make_unique<lapack_int[]>(pp);
   info = 0;
 #if defined(BLAS) || defined(CUBLAS)
   p_tmp = mxCreateDoubleMatrix(n, n, mxREAL);
-  *tmp = mxGetPr(p_tmp);
+  tmp = mxGetPr(p_tmp);
   p_P_t_t1 = mxCreateDoubleMatrix(n, n, mxREAL);
-  *P_t_t1 = mxGetPr(p_P_t_t1);
+  P_t_t1 = mxGetPr(p_P_t_t1);
   pK = mxCreateDoubleMatrix(n, n, mxREAL);
-  *K = mxGetPr(pK);
+  K = mxGetPr(pK);
   p_K_P = mxCreateDoubleMatrix(n, n, mxREAL);
-  *K_P = mxGetPr(p_K_P);
-  oldK  = static_cast<double *>(mxMalloc(n * n * sizeof(double)));
-  *P_mf = static_cast<double *>(mxMalloc(n * n * sizeof(double)));
+  K_P = mxGetPr(p_K_P);
+  oldK = std::make_unique<double[]>(n * n);
+  P_mf = std::make_unique<double[]>(n * n);
   for (int i = 0; i < n  * n; i++)
     oldK[i] = Inf;
 #else
   p_tmp = mxCreateDoubleMatrix(n, n_state, mxREAL);
-  *tmp = mxGetPr(p_tmp);
+  tmp = mxGetPr(p_tmp);
   p_P_t_t1 = mxCreateDoubleMatrix(n_state, n_state, mxREAL);
-  *P_t_t1 = mxGetPr(p_P_t_t1);
+  P_t_t1 = mxGetPr(p_P_t_t1);
   pK = mxCreateDoubleMatrix(n, pp, mxREAL);
-  *K = mxGetPr(pK);
+  K = mxGetPr(pK);
   p_K_P = mxCreateDoubleMatrix(n_state, n_state, mxREAL);
-  *K_P = mxGetPr(p_K_P);
-  oldK  = static_cast<double *>(mxMalloc(n * pp * sizeof(double)));
-  *P_mf = static_cast<double *>(mxMalloc(n * pp * sizeof(double)));
+  K_P = mxGetPr(p_K_P);
+  oldK = std::make_unique<double[]>(n * pp);
+  P_mf = std::make_unique<double[]>(n * pp);
   for (int i = 0; i < n  * pp; i++)
     oldK[i] = Inf;
 #endif
 }
 
 void
-BlockKalmanFilter::block_kalman_filter_ss(double *P_mf, double *v_pp, double *K, double *v_n, double *a, double *K_P, double *P_t_t1, double *tmp, double *P)
+BlockKalmanFilter::block_kalman_filter_ss()
 {
   if (t+1 < smpl)
-    {
-      while (t < smpl)
-        {
-          //v = Y(:,t)-a(mf);
-          for (int i = 0; i < pp; i++)
-            v[i]  = Y[i + t * pp] - a[mf[i]];
+    while (t < smpl)
+      {
+        //v = Y(:,t)-a(mf);
+        for (int i = 0; i < pp; i++)
+          v[i]  = Y[i + t * pp] - a[mf[i]];
 
-          //a = T*(a+K*v);
-          for (int i = pure_obs; i < n; i++)
-            {
-              double res = 0.0;
-              for (int j = 0; j < pp; j++)
-                res += K[j  * n + i] * v[j];
-              v_n[i] = res + a[i];
-            }
-          for (int i = 0; i < n; i++)
-            {
-              double res = 0.0;
-              for (int j = pure_obs; j < n; j++)
-                res += T[j  * n + i] * v_n[j];
-              a[i] = res;
-            }
+        //a = T*(a+K*v);
+        for (int i = pure_obs; i < n; i++)
+          {
+            double res = 0.0;
+            for (int j = 0; j < pp; j++)
+              res += K[j  * n + i] * v[j];
+            v_n[i] = res + a[i];
+          }
+        for (int i = 0; i < n; i++)
+          {
+            double res = 0.0;
+            for (int j = pure_obs; j < n; j++)
+              res += T[j  * n + i] * v_n[j];
+            a[i] = res;
+          }
 
-          //lik(t) = transpose(v)*iF*v;
-          for (int i = 0; i < pp; i++)
-            {
-              double res = 0.0;
-              for (int j = 0; j < pp; j++)
-                res += v[j] * iF[j * pp + i];
-              v_pp[i] = res;
-            }
-          double res = 0.0;
-          for (int i = 0; i < pp; i++)
-            res += v_pp[i] * v[i];
+        //lik(t) = transpose(v)*iF*v;
+        for (int i = 0; i < pp; i++)
+          {
+            double res = 0.0;
+            for (int j = 0; j < pp; j++)
+              res += v[j] * iF[j * pp + i];
+            v_pp[i] = res;
+          }
+        double res = 0.0;
+        for (int i = 0; i < pp; i++)
+          res += v_pp[i] * v[i];
 
-          lik[t] = (log(dF) + res + pp * log(2.0*pi))/2;
-          if (t + 1 > start)
-            LIK += lik[t];
+        lik[t] = (log(dF) + res + pp * log(2.0*M_PI))/2;
+        if (t + 1 > start)
+          LIK += lik[t];
 
-          t++;
-        }
-    }
+        t++;
+      }
 }
 
 bool
-BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf, double *v_pp, double *K, double *v_n, double *a, double *K_P, double *P_t_t1, double *tmp, double *P)
+BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[])
 {
   while (notsteady && t < smpl)
     {
       if (missing_observations)
         {
-
           // retrieve the d_index
           pd_index = mxGetCell(pdata_index, t);
           dd_index = static_cast<double *>(mxGetData(pd_index));
           size_d_index = mxGetM(pd_index);
           d_index.resize(size_d_index);
           for (int i = 0; i < size_d_index; i++)
-            {
-              d_index[i] = ceil(dd_index[i]) - 1;
-            }
+            d_index[i] = ceil(dd_index[i]) - 1;
 
           //v = Y(:,t) - a(mf)
           int i_i = 0;
           //#pragma omp parallel for shared(v, i_i, d_index) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
-          for (vector<int>::const_iterator i = d_index.begin(); i != d_index.end(); i++)
+          for (auto i = d_index.begin(); i != d_index.end(); i++)
             {
               //mexPrintf("i_i=%d, omp_get_max_threads()=%d\n",i_i,omp_get_max_threads());
               v[i_i]  = Y[*i + t * pp] - a[mf[*i]];
@@ -390,26 +383,21 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
           //F  = P(mf,mf) + H;
           i_i = 0;
           if (H_size == 1)
-            {
-              //#pragma omp parallel for shared(iF, F, i_i) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
-              for (vector<int>::const_iterator i = d_index.begin(); i != d_index.end(); i++, i_i++)
-                {
-                  int j_j = 0;
-                  for (vector<int>::const_iterator j = d_index.begin(); j != d_index.end(); j++, j_j++)
-                    iF[i_i + j_j * size_d_index] = F[i_i + j_j * size_d_index] = P[mf[*i] + mf[*j] * n] + H[0];
-                }
-            }
+            //#pragma omp parallel for shared(iF, F, i_i) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
+            for (auto i = d_index.begin(); i != d_index.end(); i++, i_i++)
+              {
+                int j_j = 0;
+                for (auto j = d_index.begin(); j != d_index.end(); j++, j_j++)
+                  iF[i_i + j_j * size_d_index] = F[i_i + j_j * size_d_index] = P[mf[*i] + mf[*j] * n] + H[0];
+              }
           else
-            {
-              //#pragma omp parallel for shared(iF, F, P, H, mf, i_i) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
-              for (vector<int>::const_iterator i = d_index.begin(); i != d_index.end(); i++, i_i++)
-                {
-                  int j_j = 0;
-                  for (vector<int>::const_iterator j = d_index.begin(); j != d_index.end(); j++, j_j++)
-                    iF[i_i + j_j * size_d_index] = F[i_i + j_j * size_d_index] = P[mf[*i] + mf[*j] * n] + H[*i + *j * pp];
-                }
-            }
-
+            //#pragma omp parallel for shared(iF, F, P, H, mf, i_i) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
+            for (auto i = d_index.begin(); i != d_index.end(); i++, i_i++)
+              {
+                int j_j = 0;
+                for (auto j = d_index.begin(); j != d_index.end(); j++, j_j++)
+                  iF[i_i + j_j * size_d_index] = F[i_i + j_j * size_d_index] = P[mf[*i] + mf[*j] * n] + H[*i + *j * pp];
+              }
         }
       else
         {
@@ -421,30 +409,26 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
 
           //F  = P(mf,mf) + H;
           if (H_size == 1)
-            {
-              for (int i = 0; i < pp; i++)
-                for (int j = 0; j < pp; j++)
-                  iF[i + j * pp] = F[i + j * pp] = P[mf[i] + mf[j] * n] + H[0];
-            }
+            for (int i = 0; i < pp; i++)
+              for (int j = 0; j < pp; j++)
+                iF[i + j * pp] = F[i + j * pp] = P[mf[i] + mf[j] * n] + H[0];
           else
-            {
-              for (int i = 0; i < pp; i++)
-                for (int j = 0; j < pp; j++)
-                  iF[i + j * pp] = F[i + j * pp] = P[mf[i] + mf[j] * n] + H[i + j * pp];
-            }
+            for (int i = 0; i < pp; i++)
+              for (int j = 0; j < pp; j++)
+                iF[i + j * pp] = F[i + j * pp] = P[mf[i] + mf[j] * n] + H[i + j * pp];
         }
 
       /* Computes the norm of iF */
-      double anorm = dlange("1", &size_d_index, &size_d_index, iF, &size_d_index, w);
+      double anorm = dlange("1", &size_d_index, &size_d_index, iF, &size_d_index, w.get());
       //mexPrintf("anorm = %f\n",anorm);
 
       /* Modifies F in place with a LU decomposition */
-      dgetrf(&size_d_index, &size_d_index, iF, &size_d_index, ipiv, &info);
+      dgetrf(&size_d_index, &size_d_index, iF, &size_d_index, ipiv.get(), &info);
       if (info != 0)
         mexPrintf("dgetrf failure with error %d\n", static_cast<int>(info));
 
       /* Computes the reciprocal norm */
-      dgecon("1", &size_d_index, iF, &size_d_index, &anorm, &rcond, w, iw, &info);
+      dgecon("1", &size_d_index, iF, &size_d_index, &anorm, &rcond, w.get(), iw.get(), &info);
       if (info != 0)
         mexPrintf("dgecon failure with error %d\n", static_cast<int>(info));
 
@@ -454,12 +438,10 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
             mexPrintf("error: F singular\n");
             LIK  = Inf;
             if (nlhs == 3)
-              {
-                for (int i = t; i < smpl; i++)
-                  lik[i] = Inf;
-              }
+              for (int i = t; i < smpl; i++)
+                lik[i] = Inf;
             // info = 0
-            return_results_and_clean(nlhs, plhs, &P_mf, &v_pp, &K, &K_P, &a, &K_P, &P_t_t1, &tmp, &P);
+            return_results_and_clean(nlhs, plhs);
             return false;
           }
         else
@@ -474,10 +456,10 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
                   res += T[i + j *n] * a[j];
                 tmp_a[i] = res;
               }
-            memcpy(a, tmp_a, n * sizeof(double));
+            std::copy_n(tmp_a.get(), n, a);
 
             //P = T*P*transpose(T)+QQ;
-            memset(tmp, 0, n * n_state * sizeof(double));
+            std::fill_n(tmp, 0, n * n_state);
 
             for (int i = 0; i < n; i++)
               for (int j = pure_obs; j < n; j++)
@@ -488,17 +470,15 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
                     tmp[i + j1 * n ] += T[i + k * n] * P[k + j1_n_state];
                 }
 
-            memset(P, 0, n * n * sizeof(double));
+            std::fill_n(P, 0, n * n);
             int n_n_obs = n * pure_obs;
             for (int i = 0; i < n; i++)
               for (int j = i; j < n; j++)
-                {
-                  for (int k = pure_obs; k < i_nz_state_var[j]; k++)
-                    {
-                      int k_n = k * n;
-                      P[i * n + j] += tmp[i + k_n - n_n_obs] * T[j + k_n];
-                    }
-                }
+                for (int k = pure_obs; k < i_nz_state_var[j]; k++)
+                  {
+                    int k_n = k * n;
+                    P[i * n + j] += tmp[i + k_n - n_n_obs] * T[j + k_n];
+                  }
 
             for (int i = 0; i < n; i++)
               {
@@ -513,11 +493,11 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
           F_singular = false;
 
           //dF     = det(F);
-          dF     = det(iF, size_d_index, ipiv);
+          dF     = det(iF, size_d_index, ipiv.get());
 
           //iF     = inv(F);
           //int lwork = 4/*2*/* pp;
-          dgetri(&size_d_index, iF, &size_d_index, ipiv, w, &lw, &info);
+          dgetri(&size_d_index, iF, &size_d_index, ipiv.get(), w.get(), &lw, &info);
           if (info != 0)
             mexPrintf("dgetri failure with error %d\n", static_cast<int>(info));
 
@@ -536,31 +516,27 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
           for (int i = 0; i < size_d_index; i++)
             res += v_pp[i] * v[i];
 
-          lik[t] = (log(dF) + res + size_d_index * log(2.0*pi))/2;
+          lik[t] = (log(dF) + res + size_d_index * log(2.0*M_PI))/2;
           if (t + 1 >= start)
             LIK += lik[t];
 
           if (missing_observations)
-            {
               //K      = P(:,mf)*iF;
 #ifdef USE_OMP
 # pragma omp parallel for shared(P_mf) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
 #endif
-              for (int i = 0; i < n; i++)
-                {
-                  int j_j = 0;
-                  //for (int j = 0; j < pp; j++)
-                  for (vector<int>::const_iterator j = d_index.begin(); j != d_index.end(); j++, j_j++)
-                    P_mf[i + j_j * n] = P[i + mf[*j] * n];
-                }
-            }
+            for (int i = 0; i < n; i++)
+              {
+                int j_j = 0;
+                //for (int j = 0; j < pp; j++)
+                for (auto j = d_index.begin(); j != d_index.end(); j++, j_j++)
+                  P_mf[i + j_j * n] = P[i + mf[*j] * n];
+              }
           else
-            {
-              //K      = P(:,mf)*iF;
-              for (int i = 0; i < n; i++)
-                for (int j = 0; j < pp; j++)
-                  P_mf[i + j * n] = P[i + mf[j] * n];
-            }
+            //K      = P(:,mf)*iF;
+            for (int i = 0; i < n; i++)
+              for (int j = 0; j < pp; j++)
+                P_mf[i + j * n] = P[i + mf[j] * n];
 
 #ifdef USE_OMP
 # pragma omp parallel for shared(K) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
@@ -603,20 +579,18 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
               //P      = T*(P-K*P(mf,:))*transpose(T)+QQ;
               int i_i = 0;
               //#pragma omp parallel for shared(P_mf) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
-              for (vector<int>::const_iterator i = d_index.begin(); i != d_index.end(); i++, i_i++)
+              for (auto i = d_index.begin(); i != d_index.end(); i++, i_i++)
                 for (int j = pure_obs; j < n; j++)
                   P_mf[i_i + j * size_d_index] = P[mf[*i] + j * n];
             }
           else
-            {
-              //P      = T*(P-K*P(mf,:))*transpose(T)+QQ;
+            //P      = T*(P-K*P(mf,:))*transpose(T)+QQ;
 #ifdef USE_OMP
 # pragma omp parallel for shared(P_mf) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
 #endif
-              for (int i = 0; i < pp; i++)
-                for (int j = pure_obs; j < n; j++)
-                  P_mf[i + j * pp] = P[mf[i] + j * n];
-            }
+            for (int i = 0; i < pp; i++)
+              for (int j = pure_obs; j < n; j++)
+                P_mf[i + j * pp] = P[mf[i] + j * n];
 
 #ifdef BLAS
 # ifdef USE_OMP
@@ -640,7 +614,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
               }
           double one  = 1.0;
           double zero = 0.0;
-          memcpy(P, QQ, n * n *sizeof(double));
+          std::copy_n(QQ, n * n, P);
           blas_int n_b = n;
           /*mexPrintf("sizeof(n_b)=%d, n_b=%d, sizeof(n)=%d, n=%d\n",sizeof(n_b),n_b,sizeof(n),n);
             mexEvalString("drawnow;");*/
@@ -685,27 +659,24 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
           /*int device;
             cudaGetDevice(&device);*/
           int n2 = n * n;
-          double *d_A = 0;
-          double *d_B = 0;
-          double *d_C = 0;
-          double *d_D = 0;
+          double *d_A = nullptr, *d_B = nullptr, *d_C = nullptr, *d_D = nullptr;
           // Allocate device memory for the matrices
-          if (cudaMalloc((void **) &d_A, n2 * sizeof(double)) != cudaSuccess)
+          if (cudaMalloc(static_cast<void **>(&d_A), n2 * sizeof(double)) != cudaSuccess)
             {
               mexPrintf("!!!! device memory allocation error (allocate A)\n");
               return false;
             }
-          if (cudaMalloc((void **) &d_B, n2 * sizeof(d_B[0])) != cudaSuccess)
+          if (cudaMalloc(static_cast<void **>(&d_B), n2 * sizeof(d_B[0])) != cudaSuccess)
             {
               mexPrintf("!!!! device memory allocation error (allocate B)\n");
               return false;
             }
-          if (cudaMalloc((void **) &d_C, n2 * sizeof(d_C[0])) != cudaSuccess)
+          if (cudaMalloc(static_cast<void **>(&d_C), n2 * sizeof(d_C[0])) != cudaSuccess)
             {
               mexPrintf("!!!! device memory allocation error (allocate C)\n");
               return false;
             }
-          if (cudaMalloc((void **) &d_D, n2 * sizeof(d_D[0])) != cudaSuccess)
+          if (cudaMalloc(static_cast<void **>(&d_D), n2 * sizeof(d_D[0])) != cudaSuccess)
             {
               mexPrintf("!!!! device memory allocation error (allocate D)\n");
               return false;
@@ -793,7 +764,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
                 }
             }
 
-          memset(tmp, 0, n * n_state * sizeof(double));
+          fill_n(tmp, 0, n * n_state);
 
 #  ifdef USE_OMP
 #   pragma omp parallel for shared(tmp) num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
@@ -811,7 +782,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
                 }
             }
 
-          memset(P, 0, n * n * sizeof(double));
+          fill_n(P, 0, n * n);
 
           int n_n_obs = -n * pure_obs;
 #  ifdef USE_OMP
@@ -849,14 +820,13 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
               for (int i = 0; i < n * size_d_index; i++)
                 {
                   double res = abs(K[i] - oldK[i]);
-                  if (res > max_abs)
-                    max_abs = res;
+                  max_abs = std::max(res, max_abs);
                 }
               notsteady = max_abs > riccati_tol;
 
               //oldK = K(:);
 
-              memcpy(oldK, K, n * pp * sizeof(double));
+              std::copy_n(K, n * pp, oldK.get());
             }
         }
       t++;
@@ -865,12 +835,12 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf,
   if (F_singular)
     mexErrMsgTxt("The variance of the forecast error remains singular until the end of the sample\n");
   if (t < smpl)
-    block_kalman_filter_ss(P_mf, v_pp, K, K_P, a, K_P, P_t_t1, tmp, P);
+    block_kalman_filter_ss();
   return true;
 }
 
 void
-BlockKalmanFilter::return_results_and_clean(int nlhs, mxArray *plhs[], double *P_mf[], double *v_pp[], double *K[], double *v_n[], double *a[], double *K_P[], double *P_t_t1[], double *tmp[], double *P[])
+BlockKalmanFilter::return_results_and_clean(int nlhs, mxArray *plhs[])
 {
   plhs[0] = mxCreateDoubleScalar(0);
 
@@ -886,17 +856,6 @@ BlockKalmanFilter::return_results_and_clean(int nlhs, mxArray *plhs[], double *P
   else
     mxDestroyArray(plik);
 
-  mxFree(w);
-  mxFree(i_nz_state_var);
-  mxFree(tmp_a);
-  //mxFree(v_pp);
-  //mxFree(v_n);
-  mxFree(mf);
-  mxFree(iw);
-  mxFree(ipiv);
-  mxFree(oldK);
-  //mxFree(P_mf);
-
   mxDestroyArray(pa);
   mxDestroyArray(p_tmp);
   mxDestroyArray(pQQ);
@@ -911,8 +870,7 @@ BlockKalmanFilter::return_results_and_clean(int nlhs, mxArray *plhs[], double *P
 void
 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
-  double *P_mf, *v_pp, *v_n, *a, *K, *K_P, *P_t_t1, *tmp, *P;
-  BlockKalmanFilter block_kalman_filter(nlhs, plhs, nrhs, prhs, &P_mf, &v_pp, &K, &v_n, &a, &K_P, &P_t_t1, &tmp, &P);
-  if (block_kalman_filter.block_kalman_filter(nlhs, plhs, P_mf, v_pp, K, K_P, a, K_P, P_t_t1, tmp, P))
-    block_kalman_filter.return_results_and_clean(nlhs, plhs, &P_mf, &v_pp, &K, &K_P, &a, &K_P, &P_t_t1, &tmp, &P);
+  BlockKalmanFilter block_kalman_filter(nlhs, plhs, nrhs, prhs);
+  if (block_kalman_filter.block_kalman_filter(nlhs, plhs))
+    block_kalman_filter.return_results_and_clean(nlhs, plhs);
 }
diff --git a/mex/sources/block_kalman_filter/block_kalman_filter.hh b/mex/sources/block_kalman_filter/block_kalman_filter.hh
index 52d8296ce..18c79c975 100644
--- a/mex/sources/block_kalman_filter/block_kalman_filter.hh
+++ b/mex/sources/block_kalman_filter/block_kalman_filter.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2017 Dynare Team
+ * Copyright © 2007-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -26,9 +26,11 @@
 # include "mex_interface.hh"
 #endif
 
+#include <vector>
+#include <memory>
+
 #include <dynblas.h>
 #include <dynlapack.h>
-using namespace std;
 
 class BlockKalmanFilter
 {
@@ -40,25 +42,27 @@ public:
   lapack_int pp, lw, info;
 
   double *nz_state_var;
-  int *i_nz_state_var, *mf;
+  std::unique_ptr<int[]> i_nz_state_var, mf;
   int n_diag, t;
   mxArray *M_;
   mxArray *pa, *p_tmp, *p_tmp1, *plik;
-  double *tmp_a, *tmp1, *lik, *v_n, *w, *oldK;
+  std::unique_ptr<double[]> tmp_a;
+  double *tmp1, *lik;
+  std::unique_ptr<double[]> v_n, v_pp, w, oldK, P_mf;
   bool notsteady, F_singular, missing_observations;
-  lapack_int *iw, *ipiv;
+  std::unique_ptr<lapack_int[]> iw, ipiv;
   double anorm, rcond;
   lapack_int size_d_index;
   int no_more_missing_observations, number_of_observations;
   const mxArray *pdata_index;
-  vector<int> d_index;
+  std::vector<int> d_index;
   const mxArray *pd_index;
   double *dd_index;
-
+  double *K, *a, *K_P, *P_t_t1, *tmp, *P;
 public:
-  BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[], double *P_mf[], double *v_pp[], double *K[], double *v_n[], double *a[], double *K_P[], double *P_t_t1[], double *tmp[], double *P[]);
-  bool block_kalman_filter(int nlhs, mxArray *plhs[], double *P_mf, double *v_pp, double *K, double *v_n, double *a, double *K_P, double *P_t_t1, double *tmp, double *P);
-  void block_kalman_filter_ss(double *P_mf, double *v_pp, double *K, double *v_n, double *a, double *K_P, double *P_t_t1, double *tmp, double *P);
-  void return_results_and_clean(int nlhs, mxArray *plhs[], double *P_mf[], double *v_pp[], double *K[], double *v_n[], double *a[], double *K_P[], double *P_t_t1[], double *tmp[], double *P[]);
+  BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);
+  bool block_kalman_filter(int nlhs, mxArray *plhs[]);
+  void block_kalman_filter_ss();
+  void return_results_and_clean(int nlhs, mxArray *plhs[]);
 };
 #endif