diff --git a/dynare++/integ/testing/tests.cc b/dynare++/integ/testing/tests.cc
index 3d2ecc8a3ba5ead902fa491c9b9ed005975b5b10..b99db42ddd72c572232d4631e7af24536c76dfc9 100644
--- a/dynare++/integ/testing/tests.cc
+++ b/dynare++/integ/testing/tests.cc
@@ -256,7 +256,7 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level
 
   // check against theoretical moments
   UNormalMoments moments(imom, msq);
-  smol_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
+  smol_out.add(-1.0, moments.get(Symmetry{imom}).getData());
   std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
   return smol_out.getMax() < 1.e-7;
 }
@@ -285,7 +285,7 @@ TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level
 
   // check against theoretical moments
   UNormalMoments moments(imom, msq);
-  prod_out.add(-1.0, (moments.get(Symmetry{imom}))->getData());
+  prod_out.add(-1.0, moments.get(Symmetry{imom}).getData());
   std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
   return prod_out.getMax() < 1.e-7;
 }
diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index fbda760250d572265f6699ba95809262b4238c2c..7de6b69f1554c65e5ff42a7285f22263d2935f75 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -82,7 +82,7 @@ Approximation::approxAtSteady()
 {
   model.calcDerivativesAtSteady();
   FirstOrder fo(model.nstat(), model.npred(), model.nboth(), model.nforw(),
-                model.nexog(), *(model.getModelDerivatives().get(Symmetry{1})),
+                model.nexog(), model.getModelDerivatives().get(Symmetry{1}),
                 journal, qz_criterium);
   KORD_RAISE_IF_X(!fo.isStable(),
                   "The model is not Blanchard-Kahn stable",
@@ -302,14 +302,14 @@ Approximation::calcStochShift(Vector &out, double at_sigma) const
           ten->zeros();
           for (int l = 1; l <= d; l++)
             {
-              const FSSparseTensor *f = model.getModelDerivatives().get(Symmetry{l});
-              zaux.multAndAdd(*f, *ten);
+              const FSSparseTensor &f = model.getModelDerivatives().get(Symmetry{l});
+              zaux.multAndAdd(f, *ten);
             }
 
           // multiply with shocks and add to result
           FGSTensor *tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry{0, 0, 0, 0}, nvs));
           tmp->zeros();
-          ten->contractAndAdd(1, *tmp, *(mom.get(Symmetry{d})));
+          ten->contractAndAdd(1, *tmp, mom.get(Symmetry{d}));
 
           out.add(pow(at_sigma, d)/dfac, tmp->getData());
           delete ten;
@@ -362,8 +362,8 @@ Approximation::check(double at_sigma) const
 TwoDMatrix *
 Approximation::calcYCov() const
 {
-  const TwoDMatrix &gy = *(rule_ders->get(Symmetry{1, 0, 0, 0}));
-  const TwoDMatrix &gu = *(rule_ders->get(Symmetry{0, 1, 0, 0}));
+  const TwoDMatrix &gy = rule_ders->get(Symmetry{1, 0, 0, 0});
+  const TwoDMatrix &gu = rule_ders->get(Symmetry{0, 1, 0, 0});
   TwoDMatrix G(model.numeq(), model.numeq());
   G.zeros();
   G.place(gy, 0, model.nstat());
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index 841154ebc90df836901774aedcfe1c91b28ee82a..efb9564f386431ec7f8c77ffadf71dce76c1fe42 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -224,7 +224,7 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
               if (g.check(sym))
                 {
                   double mult = pow(sigma, k)/dfact/kfact;
-                  tmp.add(mult, *(g.get(sym)));
+                  tmp.add(mult, g.get(sym));
                 }
             }
           g_yud->addSubTensor(tmp);
@@ -538,9 +538,9 @@ DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
   fillTensors(g, sigma);
   _Tparent yspol(ypart.nstat, ypart.nys(), *this);
   bigf = new _Tparent((const _Tparent &) yspol);
-  _Ttensym *frst = bigf->get(Symmetry{1});
+  _Ttensym &frst = bigf->get(Symmetry{1});
   for (int i = 0; i < ypart.nys(); i++)
-    frst->get(i, i) = frst->get(i, i) - 1;
+    frst.get(i, i) = frst.get(i, i) - 1;
   bigfder = new _Tparent(*bigf, 0);
 }
 
@@ -573,9 +573,9 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
         {
           if (g.check(Symmetry{d, 0, 0, k}))
             {
-              const _Ttensor *ten = g.get(Symmetry{d, 0, 0, k});
+              const _Ttensor &ten = g.get(Symmetry{d, 0, 0, k});
               double mult = pow(sigma, k)/dfact/kfact;
-              g_yd->add(mult, *ten);
+              g_yd->add(mult, ten);
             }
         }
       this->insert(std::move(g_yd));
diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc
index 1d3b01b7eaa63d86597228e8767fc50c55132049..dc6849a4b8c7a61a2dace024ccec16be9c6b2725 100644
--- a/dynare++/kord/global_check.cc
+++ b/dynare++/kord/global_check.cc
@@ -97,7 +97,7 @@ ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
                          model->nboth()+model->nforw());
   if (hss->check(Symmetry{0}))
     {
-      hss->get(Symmetry{0})->getData().add(1.0, ysteady_ss);
+      hss->get(Symmetry{0}).getData().add(1.0, ysteady_ss);
     }
   else
     {
diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc
index dcf3db35e79ad8ada47d1680c6e9f533bb1bc5ce..fec3b0ac8efae3a014ec5731339cecb4f34733af 100644
--- a/dynare++/kord/korder.cc
+++ b/dynare++/kord/korder.cc
@@ -239,9 +239,9 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     _um(maxk, v), _fm(_um), f(fcont),
-    matA(*(f.get(Symmetry{1})), _uZstack.getStackSizes(), gy, ypart),
-    matS(*(f.get(Symmetry{1})), _uZstack.getStackSizes(), gy, ypart),
-    matB(*(f.get(Symmetry{1})), _uZstack.getStackSizes()),
+    matA(f.get(Symmetry{1}), _uZstack.getStackSizes(), gy, ypart),
+    matS(f.get(Symmetry{1}), _uZstack.getStackSizes(), gy, ypart),
+    matB(f.get(Symmetry{1}), _uZstack.getStackSizes()),
     journal(jr)
 {
   KORD_RAISE_IF(gy.ncols() != ypart.nys(),
@@ -303,7 +303,7 @@ KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
     {
       KORD_RAISE_IF(!der.isFinite(),
                     "RHS of Sylverster is not finite");
-      TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry{1, 0, 0, 0})));
+      TwoDMatrix gs_y(gs<unfold>().get(Symmetry{1, 0, 0, 0}));
       GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(),
                             ypart.nstat+ypart.npred,
                             matA.getData(), matB.getData(),
@@ -343,7 +343,7 @@ KOrder::switchToFolded()
       {
         if (si[2] == 0 && g<unfold>().check(si))
           {
-            auto ft = std::make_unique<FGSTensor>(*(g<unfold>().get(si)));
+            auto ft = std::make_unique<FGSTensor>(g<unfold>().get(si));
             insertDerivative<fold>(std::move(ft));
             if (dim > 1)
               {
@@ -354,7 +354,7 @@ KOrder::switchToFolded()
           }
         if (G<unfold>().check(si))
           {
-            auto ft = std::make_unique<FGSTensor>(*(G<unfold>().get(si)));
+            auto ft = std::make_unique<FGSTensor>(G<unfold>().get(si));
             G<fold>().insert(std::move(ft));
             if (dim > 1)
               {
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index d4b51dcfaf17100c01a8f8dd67bafca7e2ed63e4..11cb6ebe9706eb496a4639d89fd929b781a82394 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -478,10 +478,10 @@ KOrder::recover_y(int i)
 
   insertDerivative<t>(std::move(g_yi));
 
-  _Ttensor *gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
-  gs<t>().multAndAdd(*gss_y, *G_yi_ptr);
-  _Ttensor *gss_yi = gss<t>().get(sym);
-  gs<t>().multAndAdd(*gss_yi, *G_yi_ptr);
+  _Ttensor &gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
+  gs<t>().multAndAdd(gss_y, *G_yi_ptr);
+  _Ttensor &gss_yi = gss<t>().get(sym);
+  gs<t>().multAndAdd(gss_yi, *G_yi_ptr);
 }
 
 /* Here we solve $\left[F_{y^iu^j}\right]=0$ to obtain $g_{y^iu^j}$ for
@@ -512,7 +512,7 @@ KOrder::recover_yu(int i, int j)
   matA.multInv(*g_yiuj);
   insertDerivative<t>(std::move(g_yiuj));
 
-  gs<t>().multAndAdd(*(gss<t>().get(Symmetry{1, 0, 0, 0})), *G_yiuj_ptr);
+  gs<t>().multAndAdd(gss<t>().get(Symmetry{1, 0, 0, 0}), *G_yiuj_ptr);
 }
 
 /* Here we solve
@@ -739,7 +739,7 @@ KOrder::calcD_ijk(int i, int j, int k) const
   if (is_even(k))
     {
       auto tmp = faaDiBrunoZ<t>(Symmetry{i, j, k, 0});
-      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{k})));
+      tmp->contractAndAdd(2, *res, m<t>().get(Symmetry{k}));
     }
   return res;
 }
@@ -761,7 +761,7 @@ KOrder::calcE_ijk(int i, int j, int k) const
     {
       auto tmp = faaDiBrunoZ<t>(Symmetry{i, j, n, k-n});
       tmp->mult((double) (PascalTriangle::noverk(k, n)));
-      tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry{n})));
+      tmp->contractAndAdd(2, *res, m<t>().get(Symmetry{n}));
     }
   return res;
 }
diff --git a/dynare++/kord/korder_stoch.cc b/dynare++/kord/korder_stoch.cc
index 020e389b15225f18f422807bcdeb2759733c1ff6..924f797a7cb811793f08d045750fc30566a7a0f4 100644
--- a/dynare++/kord/korder_stoch.cc
+++ b/dynare++/kord/korder_stoch.cc
@@ -35,7 +35,7 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     f(fcont),
-    matA(*(fcont.get(Symmetry{1})), _uZstack.getStackSizes(), *(hh.get(Symmetry{1, 0, 0, 0})),
+    matA(fcont.get(Symmetry{1}), _uZstack.getStackSizes(), hh.get(Symmetry{1, 0, 0, 0}),
          ypart)
 {
 }
@@ -52,7 +52,7 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
     _uGstack(&_ugs, ypart.nys(), nu),
     _fGstack(&_fgs, ypart.nys(), nu),
     f(fcont),
-    matA(*(fcont.get(Symmetry{1})), _uZstack.getStackSizes(), *(hh.get(Symmetry{1, 0, 0, 0})),
+    matA(fcont.get(Symmetry{1}), _uZstack.getStackSizes(), hh.get(Symmetry{1, 0, 0, 0}),
          ypart)
 {
 }
diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh
index 684ee36e8491488baa08a4c7e856afa5f69ab38c..7f37bc7e84491e286ee94bd4f0e6e982d175dd1e 100644
--- a/dynare++/kord/korder_stoch.hh
+++ b/dynare++/kord/korder_stoch.hh
@@ -110,12 +110,12 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
                   double mult = (pow(at_sigma, m)*povern)/mfac;
                   Symmetry sym_mn{i, m+n, 0, k};
                   if (m+n == 0 && g.check(sym_mn))
-                    ten->add(mult, *(g.get(sym_mn)));
+                    ten->add(mult, g.get(sym_mn));
                   if (m+n > 0 && KOrder::is_even(m+n) && g.check(sym_mn))
                     {
-                      _Ttensor gtmp(*(g.get(sym_mn)));
+                      _Ttensor gtmp(g.get(sym_mn));
                       gtmp.mult(mult);
-                      gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry{m+n})));
+                      gtmp.contractAndAdd(1, *ten, mom.get(Symmetry{m+n}));
                     }
                 }
             }
@@ -190,7 +190,7 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
         {
           int k = d-i;
           if (g_int.check(Symmetry{i, 0, 0, k}))
-            ten->addSubTensor(*(g_int.get(Symmetry{i, 0, 0, k})));
+            ten->addSubTensor(g_int.get(Symmetry{i, 0, 0, k}));
         }
       g_int_sym.insert(std::move(ten));
     }
@@ -231,7 +231,7 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
             {
               Symmetry sym{i, 0, 0, d-i};
               IntSequence coor(sym, pp);
-              auto ten = std::make_unique<_Ttensor>(*(g_int_cent.get(Symmetry{d})), ss, coor,
+              auto ten = std::make_unique<_Ttensor>(g_int_cent.get(Symmetry{d}), ss, coor,
                                                     TensorDimens(sym, true_nvs));
               this->insert(std::move(ten));
             }
diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc
index 10e3483e11c974996ddc7b3d44820490fca983db..baec0bb1b14cf096c54a0a7c7bf63b0a9f3314d8 100644
--- a/dynare++/kord/tests.cc
+++ b/dynare++/kord/tests.cc
@@ -289,7 +289,7 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
   for (int d = 1; d <= maxdim; d++)
     std::cout << "\ttensor fill for dim=" << d << " is:   "
               << std::setprecision(2) << std::setw(6) << std::fixed
-              << c.get(Symmetry{d})->getFillFactor()*100.0 << " %\n"
+              << c.get(Symmetry{d}).getFillFactor()*100.0 << " %\n"
               << std::defaultfloat;
   Journal jr("out.txt");
   KOrder kord(nstat, npred, nboth, nforw, c, gy, gu, v, jr);
diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc
index 2b5f85942184f2b7d7feb47d9e821d534176cef6..973e36f0e7a0a1194675364c6032d246ccc19ecd 100644
--- a/dynare++/src/dynare3.cc
+++ b/dynare++/src/dynare3.cc
@@ -353,11 +353,11 @@ DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms &a,
 void
 DynareDerEvalLoader::load(int i, int iord, const int *vars, double res)
 {
-  FSSparseTensor *t = md.get(Symmetry{iord});
+  FSSparseTensor &t = md.get(Symmetry{iord});
   IntSequence s(iord, 0);
   for (int j = 0; j < iord; j++)
     s[j] = atoms.get_pos_of_all(vars[j]);
-  t->insert(s, i, res);
+  t.insert(s, i, res);
 }
 
 DynareJacobian::DynareJacobian(Dynare &dyn)
diff --git a/dynare++/tl/cc/stack_container.cc b/dynare++/tl/cc/stack_container.cc
index 3bd2ef64f170ca74dca0e55af7b108485f266481..dc410aeb228bfd27f4db04bca2809eed228f139d 100644
--- a/dynare++/tl/cc/stack_container.cc
+++ b/dynare++/tl/cc/stack_container.cc
@@ -58,8 +58,8 @@ WorkerFoldMAADense::operator()(std::mutex &mut)
 {
   Permutation iden(dense_cont.num());
   IntSequence coor(sym, iden.getMap());
-  const FGSTensor *g = dense_cont.get(sym);
-  cont.multAndAddStacks(coor, *g, out, mut);
+  const FGSTensor &g = dense_cont.get(sym);
+  cont.multAndAddStacks(coor, g, out, mut);
 }
 
 WorkerFoldMAADense::WorkerFoldMAADense(const FoldedStackContainer &container,
@@ -407,8 +407,8 @@ WorkerUnfoldMAADense::operator()(std::mutex &mut)
 {
   Permutation iden(dense_cont.num());
   IntSequence coor(sym, iden.getMap());
-  const UGSTensor *g = dense_cont.get(sym);
-  cont.multAndAddStacks(coor, *g, out, mut);
+  const UGSTensor &g = dense_cont.get(sym);
+  cont.multAndAddStacks(coor, g, out, mut);
 }
 
 WorkerUnfoldMAADense::WorkerUnfoldMAADense(const UnfoldedStackContainer &container,
diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh
index a960e3d1f0aa9fb0381e2bcc1f2d57327e61d39e..670ac29d5c4720f1e479fa1d1222ad8271629995 100644
--- a/dynare++/tl/cc/stack_container.hh
+++ b/dynare++/tl/cc/stack_container.hh
@@ -112,7 +112,7 @@ public:
   virtual itype getType(int i, const Symmetry &s) const = 0;
   virtual int numStacks() const = 0;
   virtual bool isZero(int i, const Symmetry &s) const = 0;
-  virtual const _Ttype *getMatrix(int i, const Symmetry &s) const = 0;
+  virtual const _Ttype &getMatrix(int i, const Symmetry &s) const = 0;
   virtual int getLengthOfMatrixStacks(const Symmetry &s) const = 0;
   virtual int getUnitPos(const Symmetry &s) const = 0;
   virtual Vector *createPackedColumn(const Symmetry &s,
@@ -197,7 +197,7 @@ public:
             || (getType(i, s) == _Stype::matrix && !conts[i]->check(s)));
   }
 
-  const _Ttype *
+  const _Ttype &
   getMatrix(int i, const Symmetry &s) const override
   {
     TL_RAISE_IF(isZero(i, s) || getType(i, s) == _Stype::unit,
@@ -246,10 +246,10 @@ public:
     i = 0;
     while (i < numStacks() && getType(i, s) == _Stype::matrix)
       {
-        const _Ttype *t = getMatrix(i, s);
-        Tensor::index ind(*t, coor);
+        const _Ttype &t = getMatrix(i, s);
+        Tensor::index ind(t, coor);
         Vector subres(*res, stack_offsets[i], stack_sizes[i]);
-        subres = ConstGeneralMatrix(*t).getCol(*ind);
+        subres = ConstGeneralMatrix(t).getCol(*ind);
         i++;
       }
     if (iu != -1)
@@ -285,7 +285,7 @@ public:
              FGSTensor &out) const
   {
     if (c.check(Symmetry{dim}))
-      multAndAdd(*(c.get(Symmetry{dim})), out);
+      multAndAdd(c.get(Symmetry{dim}), out);
   }
   void multAndAdd(const FSSparseTensor &t, FGSTensor &out) const;
   void multAndAdd(int dim, const FGSContainer &c, FGSTensor &out) const;
@@ -315,7 +315,7 @@ public:
              UGSTensor &out) const
   {
     if (c.check(Symmetry{dim}))
-      multAndAdd(*(c.get(Symmetry{dim})), out);
+      multAndAdd(c.get(Symmetry{dim}), out);
   }
   void multAndAdd(const FSSparseTensor &t, UGSTensor &out) const;
   void multAndAdd(int dim, const UGSContainer &c, UGSTensor &out) const;
@@ -557,7 +557,7 @@ public:
     return stack_cont.getType(is, syms[ip]);
   }
 
-  const _Ttype *
+  const _Ttype &
   getMatrix(int is, int ip) const
   {
     return stack_cont.getMatrix(is, syms[ip]);
@@ -636,10 +636,10 @@ public:
           setUnit(i, sp.getSize(istack[i]));
         if (sp.getType(istack[i], i) == _Stype::matrix)
           {
-            const TwoDMatrix *m = sp.getMatrix(istack[i], i);
-            TL_RAISE_IF(m->nrows() != sp.getSize(istack[i]),
+            const TwoDMatrix &m = sp.getMatrix(istack[i], i);
+            TL_RAISE_IF(m.nrows() != sp.getSize(istack[i]),
                         "Wrong size of returned matrix in KronProdStack constructor");
-            setMat(i, *m);
+            setMat(i, m);
           }
       }
   }
diff --git a/dynare++/tl/cc/t_container.hh b/dynare++/tl/cc/t_container.hh
index 7f39749fbc4a1da67f978e950aba118a8b1fd4ea..c8238863ad736ba75785cd7023eed615d30c1136 100644
--- a/dynare++/tl/cc/t_container.hh
+++ b/dynare++/tl/cc/t_container.hh
@@ -62,6 +62,8 @@
 #include <memory>
 #include <utility>
 
+#include <cstring>
+
 #include <matio.h>
 
 // |ltsym| predicate
@@ -96,8 +98,6 @@ template<class _Ttype>
 class TensorContainer
 {
 protected:
-  using _const_ptr = const _Ttype *;
-  using _ptr = _Ttype *;
   using _Map = std::map<Symmetry, std::unique_ptr<_Ttype>, ltsym>;
 public:
   using iterator = typename _Map::iterator;
@@ -131,34 +131,24 @@ public:
       insert(std::make_unique<_Ttype>(first_row, num, *(it.second)));
   }
 
-  _const_ptr
+  const _Ttype &
   get(const Symmetry &s) const
   {
     TL_RAISE_IF(s.num() != num(),
                 "Incompatible symmetry lookup in TensorContainer::get");
     auto it = m.find(s);
-    if (it == m.end())
-      {
-        TL_RAISE("Symmetry not found in TensorContainer::get");
-        return nullptr;
-      }
-    else
-      return it->second.get();
+    TL_RAISE_IF(it == m.end(), "Symmetry not found in TensorContainer::get");
+    return *(it->second);
   }
 
-  _ptr
+  _Ttype &
   get(const Symmetry &s)
   {
     TL_RAISE_IF(s.num() != num(),
                 "Incompatible symmetry lookup in TensorContainer::get");
     auto it = m.find(s);
-    if (it == m.end())
-      {
-        TL_RAISE("Symmetry not found in TensorContainer::get");
-        return nullptr;
-      }
-    else
-      return it->second.get();
+    TL_RAISE_IF(it == m.end(), "Symmetry not found in TensorContainer::get");
+    return *(it->second);
   }
 
   bool
@@ -260,16 +250,16 @@ public:
      through all equivalence classes, calculate implied symmetry, and
      fetch its tensor storing it in the same order to the vector. */
 
-  std::vector<_const_ptr>
+  std::vector<const _Ttype *>
   fetchTensors(const Symmetry &rsym, const Equivalence &e) const
   {
-    std::vector<_const_ptr> res(e.numClasses());
+    std::vector<const _Ttype *> res(e.numClasses());
     int i = 0;
     for (auto it = e.begin();
          it != e.end(); ++it, i++)
       {
         Symmetry s(rsym, *it);
-        res[i] = get(s);
+        res[i] = &get(s);
       }
     return res;
   }
diff --git a/dynare++/tl/cc/t_polynomial.hh b/dynare++/tl/cc/t_polynomial.hh
index 0745dedfeb0405c9865127002d21eedcdeb50d10..3ec0b224e0650765b3202573d5473c986002c0ea 100644
--- a/dynare++/tl/cc/t_polynomial.hh
+++ b/dynare++/tl/cc/t_polynomial.hh
@@ -99,7 +99,6 @@ class TensorPolynomial : public TensorContainer<_Ttype>
   int nv;
   int maxdim;
   using _Tparent = TensorContainer<_Ttype>;
-  using _ptr = typename _Tparent::_ptr;
 public:
   TensorPolynomial(int rows, int vars)
     : TensorContainer<_Ttype>(1),
@@ -175,7 +174,7 @@ public:
                 /* The pointer |ten| is either a new tensor or got from |this| container. */
                 _Ttype *ten;
                 if (_Tparent::check(Symmetry{j}))
-                  ten = _Tparent::get(Symmetry{j});
+                  ten = &_Tparent::get(Symmetry{j});
                 else
                   {
                     auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
@@ -186,7 +185,7 @@ public:
 
                 Symmetry sym{i, j};
                 IntSequence coor(sym, pp);
-                _TGStype slice(*(tp.get(Symmetry{i+j})), ss, coor, TensorDimens(sym, ss));
+                _TGStype slice(tp.get(Symmetry{i+j}), ss, coor, TensorDimens(sym, ss));
                 slice.mult(PascalTriangle::noverk(i+j, j));
                 _TGStype tmp(*ten);
                 slice.contractAndAdd(0, tmp, xpow);
@@ -207,7 +206,7 @@ public:
             /* Same code as above */
             _Ttype *ten;
             if (_Tparent::check(Symmetry{j}))
-              ten = _Tparent::get(Symmetry{j});
+              ten = &_Tparent::get(Symmetry{j});
             else
               {
                 auto ten_smart = std::make_unique<_Ttype>(nrows(), nvars(), j);
@@ -218,7 +217,7 @@ public:
 
             Symmetry sym{0, j};
             IntSequence coor(sym, pp);
-            _TGStype slice(*(tp.get(Symmetry{j})), ss, coor, TensorDimens(sym, ss));
+            _TGStype slice(tp.get(Symmetry{j}), ss, coor, TensorDimens(sym, ss));
             ten->add(1.0, slice);
           }
       }
@@ -247,7 +246,7 @@ public:
   evalTrad(Vector &out, const ConstVector &v) const
   {
     if (_Tparent::check(Symmetry{0}))
-      out = _Tparent::get(Symmetry{0})->getData();
+      out = _Tparent::get(Symmetry{0}).getData();
     else
       out.zeros();
 
@@ -258,8 +257,8 @@ public:
         Symmetry cs{d};
         if (_Tparent::check(cs))
           {
-            const _Ttype *t = _Tparent::get(cs);
-            t->multaVec(out, p.getData());
+            const _Ttype &t = _Tparent::get(cs);
+            t.multaVec(out, p.getData());
           }
       }
   }
@@ -271,7 +270,7 @@ public:
   evalHorner(Vector &out, const ConstVector &v) const
   {
     if (_Tparent::check(Symmetry{0}))
-      out = _Tparent::get(Symmetry{0})->getData();
+      out = _Tparent::get(Symmetry{0}).getData();
     else
       out.zeros();
 
@@ -280,16 +279,16 @@ public:
 
     _Ttype *last;
     if (maxdim == 1)
-      last = new _Ttype(*(_Tparent::get(Symmetry{1})));
+      last = new _Ttype(_Tparent::get(Symmetry{1}));
     else
-      last = new _Ttype(*(_Tparent::get(Symmetry{maxdim})), v);
+      last = new _Ttype(_Tparent::get(Symmetry{maxdim}), v);
     for (int d = maxdim-1; d >= 1; d--)
       {
         Symmetry cs{d};
         if (_Tparent::check(cs))
           {
-            const _Ttype *nt = _Tparent::get(cs);
-            last->add(1.0, ConstTwoDMatrix(*nt));
+            const _Ttype &nt = _Tparent::get(cs);
+            last->add(1.0, ConstTwoDMatrix(nt));
           }
         if (d > 1)
           {
@@ -338,8 +337,8 @@ public:
       {
         if (_Tparent::check(Symmetry{d}))
           {
-            _Ttype *ten = _Tparent::get(Symmetry{d});
-            ten->mult((double) std::max((d-k), 0));
+            _Ttype &ten = _Tparent::get(Symmetry{d});
+            ten.mult((double) std::max((d-k), 0));
           }
       }
   }
@@ -367,13 +366,13 @@ public:
     res->zeros();
 
     if (_Tparent::check(Symmetry{s}))
-      res->add(1.0, *(_Tparent::get(Symmetry{s})));
+      res->add(1.0, _Tparent::get(Symmetry{s}));
 
     for (int d = s+1; d <= maxdim; d++)
       {
         if (_Tparent::check(Symmetry{d}))
           {
-            const _Ttype &ltmp = *(_Tparent::get(Symmetry{d}));
+            const _Ttype &ltmp = _Tparent::get(Symmetry{d});
             auto last = std::make_unique<_Ttype>(ltmp);
             for (int j = 0; j < d - s; j++)
               {
@@ -474,7 +473,7 @@ public:
         if (pol.check(Symmetry{d}))
           {
             TwoDMatrix subt(*this, offset, dumgs.ncols());
-            subt.add(1.0, *(pol.get(Symmetry{d})));
+            subt.add(1.0, pol.get(Symmetry{d}));
           }
         offset += dumgs.ncols();
       }
diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc
index 6bbb97f8527e5e775db62ef053d00a4684263542..e88a3f3066c7c800f9b840bf5812d5d38c99ec7e 100644
--- a/dynare++/tl/testing/tests.cc
+++ b/dynare++/tl/testing/tests.cc
@@ -269,8 +269,8 @@ TestRunnable::folded_monomial(int ng, int nx, int ny, int nu, int dim)
       stime = clock() - stime;
       std::cout << "\t\ttime for symmetry: "
                 << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
-      const FGSTensor *mres = gen.rcont->get(s);
-      res.add(-1.0, *mres);
+      const FGSTensor &mres = gen.rcont->get(s);
+      res.add(-1.0, mres);
       double normtmp = res.getData().getMax();
       std::cout << "\t\terror normMax:     " << normtmp << '\n';
       if (normtmp > maxnorm)
@@ -307,9 +307,9 @@ TestRunnable::unfolded_monomial(int ng, int nx, int ny, int nu, int dim)
       stime = clock() - stime;
       std::cout << "\t\ttime for symmetry: "
                 << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
-      const FGSTensor *mres = gen.rcont->get(s);
+      const FGSTensor &mres = gen.rcont->get(s);
       FGSTensor foldres(res);
-      foldres.add(-1.0, *mres);
+      foldres.add(-1.0, mres);
       double normtmp = foldres.getData().getMax();
       std::cout << "\t\terror normMax:     " << normtmp << '\n';
       if (normtmp > maxnorm)
@@ -353,8 +353,8 @@ TestRunnable::fold_zcont(int nf, int ny, int nu, int nup, int nbigg,
         stime = clock() - stime;
         std::cout << "\t\ttime for symmetry: "
                   << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
-        const FGSTensor *mres = dg.rcont.get(si);
-        res.add(-1.0, *mres);
+        const FGSTensor &mres = dg.rcont.get(si);
+        res.add(-1.0, mres);
         double normtmp = res.getData().getMax();
         std::cout << "\t\terror normMax:     " << normtmp << '\n';
         if (normtmp > maxnorm)
@@ -407,8 +407,8 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg,
         std::cout << "\t\ttime for symmetry: "
                   << static_cast<double>(stime)/CLOCKS_PER_SEC << '\n';
         FGSTensor fold_res(res);
-        const FGSTensor *mres = dg.rcont.get(si);
-        fold_res.add(-1.0, *mres);
+        const FGSTensor &mres = dg.rcont.get(si);
+        fold_res.add(-1.0, mres);
         double normtmp = fold_res.getData().getMax();
         std::cout << "\t\terror normMax:     " << normtmp << '\n';
         if (normtmp > maxnorm)
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index f2581c79b670fd4cf7299f4c2247eee656c3753e..cb81d608b3975356a0569b5a3c4feee7084497d1 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -65,11 +65,11 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, std
 void
 copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer *derivs, const std::string &fieldname)
 {
-  const TwoDMatrix *x = derivs->get(sym);
-  int n = x->numRows();
-  int m = x->numCols();
+  const TwoDMatrix &x = derivs->get(sym);
+  int n = x.numRows();
+  int m = x.numCols();
   mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
-  memcpy(mxGetPr(tmp), x->getData().base(), n*m*sizeof(double));
+  memcpy(mxGetPr(tmp), x.getData().base(), n*m*sizeof(double));
   mxSetField(destin, 0, fieldname.c_str(), tmp);
 }