diff --git a/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.cpp b/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.cpp
index 7f539692b09ce0fa3ce8b263255503037de6e99b..959f89d6871fd65c8d0a8c15238e6c6d62cbe7ce 100644
--- a/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.cpp
+++ b/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.cpp
@@ -35,17 +35,6 @@ namespace MoReFEM::TestNS::PetscNS::MatrixOperationsNS
         decltype(auto) numbering_subset = god_of_dof.GetNumberingSubset(sole);
         decltype(auto) mpi = god_of_dof.GetMpi();
 
-        vector_init_1_ = std::make_unique<GlobalVector>(numbering_subset);
-        vector_init_2_ = std::make_unique<GlobalVector>(numbering_subset);
-        vector_init_3_ = std::make_unique<GlobalVector>(numbering_subset);
-        
-        auto& vector_init_1 = *vector_init_1_;
-        auto& vector_init_2 = *vector_init_2_;
-        auto& vector_init_3 = *vector_init_3_;
-        
-        AllocateGlobalVector(god_of_dof, vector_init_1);
-        AllocateGlobalVector(god_of_dof, vector_init_2);
-        AllocateGlobalVector(god_of_dof, vector_init_3);
 
         for (auto i = 0ul; i < Utilities::ArraySize<decltype(initialized_matrix_list_)>::GetValue(); ++i)
         {
@@ -54,36 +43,18 @@ namespace MoReFEM::TestNS::PetscNS::MatrixOperationsNS
         }
 
 
-        matrix_init_2_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_init_3_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        
-        matrix_4_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_5_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_6_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_7_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_8_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_9_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        matrix_10_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
-        
-        auto& matrix_init_2 = *matrix_init_2_;
-        auto& matrix_init_3 = *matrix_init_3_;
-        
-        // no initialization for the following matrices:
-        auto& matrix_4 = *matrix_4_;
-        auto& matrix_5 = *matrix_5_;
-        auto& matrix_6 = *matrix_6_;
-        auto& matrix_7 = *matrix_7_;
-        auto& matrix_8 = *matrix_8_;
-        auto& matrix_9 = *matrix_9_;
-        auto& matrix_10 = *matrix_10_;
-        
-        // only three matrices allocated, all the others are set to PETSC_NULL
-        AllocateGlobalMatrix(god_of_dof, matrix_init_2);
-        AllocateGlobalMatrix(god_of_dof, matrix_init_3);
+        for (auto i = 0ul; i < Utilities::ArraySize<decltype(initialized_vector_list_)>::GetValue(); ++i)
+        {
+            initialized_vector_list_[i] = std::make_unique<GlobalVector>(numbering_subset);
+            AllocateGlobalVector(god_of_dof, *initialized_vector_list_[i]);
+        }
 
-        // Init vector_init_1
+
+        // Put values into the first vector.
         {
-            Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_and_write> content(vector_init_1,
+            auto& vector = GetNonCstInitializedVector<0>();
+
+            Wrappers::Petsc::AccessVectorContent<Utilities::Access::read_and_write> content(vector,
                                                                                             __FILE__, __LINE__);
 
             const auto size = content.GetSize(__FILE__, __LINE__);
diff --git a/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hpp b/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hpp
index 0d198384763e13bdd27c3ec2278aa4a40aee8780..329143436f2d94ae6a358d0c83bef535a95b75a5 100644
--- a/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hpp
+++ b/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hpp
@@ -145,6 +145,29 @@ namespace MoReFEM
             GlobalMatrix& GetNonCstInitializedMatrix() noexcept;
 
 
+            /*!
+             * \brief Constant access to a vector that was fully initialized (i.e. with layout also defined, not just the constructor called).
+             *
+             * \tparam I Index of the requested vector (some tests require several so several are stored). There are static check for the consistency of this
+             * index.
+             *
+             * \return Constant accessor to an initialized vector.
+             */
+            template<std::size_t I>
+            const GlobalVector& GetInitializedVector() const noexcept;
+
+            /*!
+             * \brief Non-constant access to a vector that was fully initialized (i.e. with layout also defined, not just the constructor called).
+             *
+             * \tparam I Index of the requested vector (some tests require several so several are stored). There are static check for the consistency of this
+             * index.
+             *
+             * \return Non-constant accessor to an initialized vector.
+             */
+            template<std::size_t I>
+            GlobalVector& GetNonCstInitializedVector() noexcept;
+
+
             //! Get the \a NumberingSubset used in the tests.
             const NumberingSubset& GetNumberingSubset() const noexcept;
 
@@ -169,49 +192,14 @@ namespace MoReFEM
         private:
 
             
-            //! Global vector which is used in the tests.
-            GlobalVector::unique_ptr vector_init_1_ = nullptr;
-            
-            //! Global vector which is used in the tests.
-            GlobalVector::unique_ptr vector_init_2_ = nullptr;
-
-            //! Global vector which is used in the tests.
-            GlobalVector::unique_ptr vector_init_3_ = nullptr;
-
-            /// \name Global matrices used by the PETSc wrappers
-            ///@{
-
-            //! Matrix used in the test.
-            GlobalMatrix::array_unique_ptr<2> initialized_matrix_list_
-                = Utilities::NullptrArray<GlobalMatrix::unique_ptr, 2>();
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_init_2_ = nullptr;
+            //! Matrices used in the test that have been fully initialized (mpi layout is assigned).
+            GlobalMatrix::array_unique_ptr<3> initialized_matrix_list_
+                = Utilities::NullptrArray<GlobalMatrix::unique_ptr, 3>();
 
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_init_3_ = nullptr;
+            //! Vectors used in the test that have been fully initialized (mpi layout is assigned).
+            GlobalVector::array_unique_ptr<3> initialized_vector_list_
+                = Utilities::NullptrArray<GlobalVector::unique_ptr, 3>();
 
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_4_ = nullptr;
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_5_ = nullptr;
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_6_ = nullptr;
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_7_ = nullptr;
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_8_ = nullptr;
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_9_ = nullptr;
-
-            //! Matrix used in the test.
-            GlobalMatrix::unique_ptr matrix_10_ = nullptr;
-            ///@}
             
         };
 
diff --git a/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hxx b/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hxx
index 88909cff418385d67ec67788634f9240d8f54575..a5367434840548482264fc30ea8c560413f70599 100644
--- a/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hxx
+++ b/Sources/Test/ThirdParty/PETSc/MatrixOperations/Model.hxx
@@ -59,6 +59,22 @@ namespace MoReFEM
         }
 
 
+        template<std::size_t I>
+        inline const GlobalVector& Model::GetInitializedVector() const noexcept
+        {
+            static_assert(I < Utilities::ArraySize<decltype(initialized_vector_list_)>::GetValue());
+            assert(!(!initialized_vector_list_[I]));
+            return *initialized_vector_list_[I];
+        }
+
+
+        template<std::size_t I>
+        inline GlobalVector& Model::GetNonCstInitializedVector() noexcept
+        {
+            return const_cast<GlobalVector&>(GetInitializedVector<I>());
+        }
+
+
     } // namespace  TestNS::PetscNS::MatrixOperationsNS
 
 
diff --git a/Sources/Test/ThirdParty/PETSc/MatrixOperations/main.cpp b/Sources/Test/ThirdParty/PETSc/MatrixOperations/main.cpp
index d14ea62190e8b13e37e6fbf2c40d662a646ce295..a153e0a7612adeb04d676f572766cd665bf78c5d 100644
--- a/Sources/Test/ThirdParty/PETSc/MatrixOperations/main.cpp
+++ b/Sources/Test/ThirdParty/PETSc/MatrixOperations/main.cpp
@@ -87,177 +87,202 @@ BOOST_FIXTURE_TEST_CASE(axpy, fixture_type)
 }
 
 
-//BOOST_FIXTURE_TEST_CASE(mat_shift, fixture_type)
-//{
-//    decltype(auto) model = GetModel();
-//
-//       Wrappers::Petsc::MatShift(1.,
-//                                 model.GetInitializedMatrix<0>(),
-//                                 __FILE__, __LINE__);
-//}
+BOOST_FIXTURE_TEST_CASE(mat_shift, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
 
+       Wrappers::Petsc::MatShift(1.,
+                                 model.GetNonCstInitializedMatrix<0>(),
+                                 __FILE__, __LINE__);
+}
 
 
+BOOST_FIXTURE_TEST_CASE(mat_mult, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
 
-/*
+    Wrappers::Petsc::MatMult(model.GetInitializedMatrix<0>(),
+                             model.GetInitializedVector<0>(),
+                             model.GetNonCstInitializedVector<1>(),
+                             __FILE__, __LINE__);
+}
 
-  {
+
+BOOST_FIXTURE_TEST_CASE(mat_mult_add, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    Wrappers::Petsc::MatMultAdd(model.GetInitializedMatrix<0>(),
+                                model.GetInitializedVector<0>(),
+                                model.GetInitializedVector<1>(),
+                                model.GetNonCstInitializedVector<2>(),
+                                __FILE__, __LINE__);
+}
 
 
+BOOST_FIXTURE_TEST_CASE(mat_mult_transpose, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
 
-             Wrappers::Petsc::MatMult(matrix_init_1,
-                                      vector_init_1,
-                                      vector_init_2,
+    Wrappers::Petsc::MatMultTranspose(model.GetInitializedMatrix<0>(),
+                                      model.GetInitializedVector<0>(),
+                                      model.GetNonCstInitializedVector<1>(),
                                       __FILE__, __LINE__);
+}
 
-             Wrappers::Petsc::MatMultAdd(matrix_init_1,
-                                         vector_init_1,
-                                         vector_init_2,
-                                         vector_init_3,
+
+BOOST_FIXTURE_TEST_CASE(mat_mult_transpose_add, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    Wrappers::Petsc::MatMultTransposeAdd(model.GetInitializedMatrix<0>(),
+                                         model.GetInitializedVector<0>(),
+                                         model.GetInitializedVector<1>(),
+                                         model.GetNonCstInitializedVector<2>(),
                                          __FILE__, __LINE__);
+}
+
+
+BOOST_FIXTURE_TEST_CASE(mat_transpose, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    decltype(auto) numbering_subset = model.GetNumberingSubset();
+    GlobalMatrix non_allocated_matrix(numbering_subset, numbering_subset);
+
+    Wrappers::Petsc::MatTranspose(model.GetInitializedMatrix<0>(),
+                                  non_allocated_matrix,
+                                  __FILE__, __LINE__,
+                                  Wrappers::Petsc::DoReuseMatrix::no);
+
+    Wrappers::Petsc::MatTranspose(model.GetInitializedMatrix<0>(),
+                                  non_allocated_matrix,
+                                  __FILE__, __LINE__,
+                                  Wrappers::Petsc::DoReuseMatrix::yes);
+}
+
+
+
+BOOST_FIXTURE_TEST_CASE(mat_transpose_inplace, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    Wrappers::Petsc::MatTranspose(model.GetNonCstInitializedMatrix<0>(),
+                                  model.GetNonCstInitializedMatrix<0>(),
+                                  __FILE__, __LINE__,
+                                  Wrappers::Petsc::DoReuseMatrix::in_place);
+}
 
-             Wrappers::Petsc::MatMultTranspose(matrix_init_1,
-                                               vector_init_1,
-                                               vector_init_2,
-                                               __FILE__, __LINE__);
-
-             Wrappers::Petsc::MatMultTransposeAdd(matrix_init_1,
-                                                  vector_init_1,
-                                                  vector_init_2,
-                                                  vector_init_3,
-                                                  __FILE__, __LINE__);
-
-             // Petsc wrappers that have a DoReuseMatrix policy
-             {
-                 Wrappers::Petsc::MatTranspose(matrix_init_1,
-                                               matrix_5,
-                                               __FILE__, __LINE__,
-                                               Wrappers::Petsc::DoReuseMatrix::no);
-
-                 Wrappers::Petsc::MatTranspose(matrix_init_1,
-                                               matrix_5,
-                                               __FILE__, __LINE__,
-                                               Wrappers::Petsc::DoReuseMatrix::yes);
-
-                 Wrappers::Petsc::MatTranspose(matrix_init_1,
-                                               matrix_init_1,
-                                               __FILE__, __LINE__,
-                                               Wrappers::Petsc::DoReuseMatrix::in_place);
-             }
-             {
-                 Wrappers::Petsc::MatMatMult(matrix_init_1,
-                                             matrix_init_2,
-                                             matrix_6,
+
+BOOST_FIXTURE_TEST_CASE(mat_mat_mult, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    decltype(auto) numbering_subset = model.GetNumberingSubset();
+    GlobalMatrix non_allocated_matrix(numbering_subset, numbering_subset);
+
+    Wrappers::Petsc::MatMatMult(model.GetInitializedMatrix<0>(),
+                                model.GetInitializedMatrix<1>(),
+                                non_allocated_matrix,
+                                __FILE__, __LINE__,
+                                Wrappers::Petsc::DoReuseMatrix::no);
+
+    Wrappers::Petsc::MatMatMult(model.GetInitializedMatrix<0>(),
+                                model.GetInitializedMatrix<1>(),
+                                non_allocated_matrix,
+                                __FILE__, __LINE__,
+                                Wrappers::Petsc::DoReuseMatrix::yes);
+}
+
+
+BOOST_FIXTURE_TEST_CASE(mat_transpose_mat_mult, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    decltype(auto) numbering_subset = model.GetNumberingSubset();
+    GlobalMatrix non_allocated_matrix(numbering_subset, numbering_subset);
+
+    Wrappers::Petsc::MatTransposeMatMult(model.GetInitializedMatrix<0>(),
+                                         model.GetInitializedMatrix<1>(),
+                                         non_allocated_matrix,
+                                         __FILE__, __LINE__,
+                                         Wrappers::Petsc::DoReuseMatrix::no);
+
+    Wrappers::Petsc::MatTransposeMatMult(model.GetInitializedMatrix<0>(),
+                                         model.GetInitializedMatrix<1>(),
+                                         non_allocated_matrix,
+                                         __FILE__, __LINE__,
+                                         Wrappers::Petsc::DoReuseMatrix::yes);
+}
+
+
+BOOST_FIXTURE_TEST_CASE(mat_mat_transpose_mat_mult, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+    decltype(auto) mpi = model.GetMpi();
+
+    decltype(auto) numbering_subset = model.GetNumberingSubset();
+    GlobalMatrix non_allocated_matrix(numbering_subset, numbering_subset);
+
+    if (mpi.Nprocessor<int>() == 1)
+    {
+        Wrappers::Petsc::MatMatTransposeMult(model.GetInitializedMatrix<0>(),
+                                             model.GetInitializedMatrix<1>(),
+                                             non_allocated_matrix,
                                              __FILE__, __LINE__,
                                              Wrappers::Petsc::DoReuseMatrix::no);
 
-                 Wrappers::Petsc::MatMatMult(matrix_init_1,
-                                             matrix_init_2,
-                                             matrix_6,
+        Wrappers::Petsc::MatMatTransposeMult(model.GetInitializedMatrix<0>(),
+                                             model.GetInitializedMatrix<1>(),
+                                             non_allocated_matrix,
                                              __FILE__, __LINE__,
                                              Wrappers::Petsc::DoReuseMatrix::yes);
-             }
-             {
-                 Wrappers::Petsc::MatTransposeMatMult(matrix_init_1,
-                                                      matrix_init_2,
-                                                      matrix_7,
-                                                      __FILE__, __LINE__,
-                                                      Wrappers::Petsc::DoReuseMatrix::no);
-
-                 Wrappers::Petsc::MatTransposeMatMult(matrix_init_1,
-                                                      matrix_init_2,
-                                                      matrix_7,
-                                                      __FILE__, __LINE__,
-                                                      Wrappers::Petsc::DoReuseMatrix::yes);
-             }
-             {
-                 Wrappers::Petsc::MatMatTransposeMult(matrix_init_1,
-                                                      matrix_init_2,
-                                                      matrix_8,
-                                                      __FILE__, __LINE__,
-                                                      Wrappers::Petsc::DoReuseMatrix::no);
-
-                 Wrappers::Petsc::MatMatTransposeMult(matrix_init_1,
-                                                      matrix_init_2,
-                                                      matrix_8,
-                                                      __FILE__, __LINE__,
-                                                      Wrappers::Petsc::DoReuseMatrix::yes);
-             }
-             {
-                 Wrappers::Petsc::PtAP(matrix_init_1,
-                                       matrix_init_2,
-                                       matrix_9,
-                                       __FILE__, __LINE__,
-                                       Wrappers::Petsc::DoReuseMatrix::no);
-
-                 Wrappers::Petsc::PtAP(matrix_init_1,
-                                       matrix_init_2,
-                                       matrix_9,
-                                       __FILE__, __LINE__,
-                                       Wrappers::Petsc::DoReuseMatrix::yes);
-             }
-             {
-                 Wrappers::Petsc::MatMatMatMult(matrix_init_1,
-                                                matrix_init_2,
-                                                matrix_init_3,
-                                                matrix_10,
-                                                __FILE__, __LINE__,
-                                                Wrappers::Petsc::DoReuseMatrix::no);
-
-                 Wrappers::Petsc::MatMatMatMult(matrix_init_1,
-                                                matrix_init_2,
-                                                matrix_init_3,
-                                                matrix_10,
-                                                __FILE__, __LINE__,
-                                                Wrappers::Petsc::DoReuseMatrix::yes);
-             }
-             {
-                 IS rperm;
-                 IS cperm;
-
-                 Wrappers::Petsc::GetOrdering(matrix_init_1,
-                                              MATORDERINGNATURAL,
-                                              &rperm,
-                                              &cperm,
-                                              __FILE__, __LINE__);
-                 ISDestroy(&rperm);
-                 ISDestroy(&cperm);
-             }
-
-             // \todo
- //
- //            {
- //                IS perm;
- //
- //                Wrappers::Petsc::CholeskyFactor(matrix_init_1,
- //                                                perm,
- //                                                NULL,
- //                                                __FILE__, __LINE__);
- //            }
- //
- //            {
- //                IS row;
- //                IS col;
- //
- //                Wrappers::Petsc::LUFactor(matrix_init_1,
- //                                          NULL,
- //                                          NULL,
- //                                          NULL,
- //                                          __FILE__, __LINE__);
- //            }
- //
- //            {
- //
- //                        Wrappers::Petsc::MatMatSolve(matrix_init_1,
- //                                                     matrix_init_2,
- //                                                     matrix_4,
- //                                                     __FILE__, __LINE__);
- //
- //            }
+    }
+}
 
 
+BOOST_FIXTURE_TEST_CASE(pt_a_p, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
 
- */
+    decltype(auto) numbering_subset = model.GetNumberingSubset();
+    GlobalMatrix non_allocated_matrix(numbering_subset, numbering_subset);
+
+    Wrappers::Petsc::PtAP(model.GetInitializedMatrix<0>(),
+                                         model.GetInitializedMatrix<1>(),
+                                         non_allocated_matrix,
+                                         __FILE__, __LINE__,
+                                         Wrappers::Petsc::DoReuseMatrix::no);
+
+    Wrappers::Petsc::PtAP(model.GetInitializedMatrix<0>(),
+                                         model.GetInitializedMatrix<1>(),
+                                         non_allocated_matrix,
+                                         __FILE__, __LINE__,
+                                         Wrappers::Petsc::DoReuseMatrix::yes);
+}
+
+
+BOOST_FIXTURE_TEST_CASE(mat_mat_mat_mult, fixture_type)
+{
+    decltype(auto) model = GetNonCstModel();
+
+    decltype(auto) numbering_subset = model.GetNumberingSubset();
+    GlobalMatrix non_allocated_matrix(numbering_subset, numbering_subset);
+
+    Wrappers::Petsc::MatMatMatMult(model.GetInitializedMatrix<0>(),
+                                   model.GetInitializedMatrix<1>(),
+                                   model.GetInitializedMatrix<2>(),
+                                   non_allocated_matrix,
+                                   __FILE__, __LINE__,
+                                   Wrappers::Petsc::DoReuseMatrix::no);
+
+    Wrappers::Petsc::MatMatMatMult(model.GetInitializedMatrix<0>(),
+                                   model.GetInitializedMatrix<1>(),
+                                   model.GetInitializedMatrix<2>(),
+                                   non_allocated_matrix,
+                                   __FILE__, __LINE__,
+                                   Wrappers::Petsc::DoReuseMatrix::yes);
+}
 
 
 PRAGMA_DIAGNOSTIC(pop)