From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 6 Dec 2017 20:00:15 +0000 (-0700)
Subject: Provide more instantiations for std::complex data types.
X-Git-Tag: v9.0.0-rc1~674^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5590%2Fhead;p=dealii.git

Provide more instantiations for std::complex data types.
---

diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h
index a58fe540fc..7c12c0d0c5 100644
--- a/include/deal.II/lac/constraint_matrix.h
+++ b/include/deal.II/lac/constraint_matrix.h
@@ -1464,13 +1464,13 @@ private:
   /**
    * Internal helper function for distribute_local_to_global function.
    */
-  template <typename LocalType>
-  LocalType
+  template <typename MatrixScalar, typename VectorScalar>
+  typename ProductType<VectorScalar,MatrixScalar>::type
   resolve_vector_entry (const size_type                       i,
                         const internals::GlobalRowsFromLocal &global_rows,
-                        const Vector<LocalType>              &local_vector,
+                        const Vector<VectorScalar>           &local_vector,
                         const std::vector<size_type>         &local_dof_indices,
-                        const FullMatrix<LocalType>          &local_matrix) const;
+                        const FullMatrix<MatrixScalar>       &local_matrix) const;
 };
 
 
diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h
index 6c397f2e24..4c732d2e2f 100644
--- a/include/deal.II/lac/constraint_matrix.templates.h
+++ b/include/deal.II/lac/constraint_matrix.templates.h
@@ -1399,7 +1399,7 @@ namespace internals
    * Access to the scratch data is only through the accessor class which
    * handles the access as well as marking the data as used.
    */
-  template <typename Number>
+  template <typename MatrixScalar, typename VectorScalar = MatrixScalar>
   class ConstraintMatrixData
   {
   public:
@@ -1434,7 +1434,7 @@ namespace internals
       /**
        * Temporary array for column values
        */
-      std::vector<Number>    values;
+      std::vector<MatrixScalar>    values;
 
       /**
        * Temporary array for block start indices
@@ -1449,7 +1449,7 @@ namespace internals
       /**
        * Temporary array for vector values
        */
-      std::vector<Number> vector_values;
+      std::vector<VectorScalar> vector_values;
 
       /**
        * Data array for reorder row/column indices.
@@ -2187,19 +2187,19 @@ make_sorted_row_list (const std::vector<size_type> &local_dof_indices,
 
 
 // Resolve the constraints from the vector and apply inhomogeneities.
-template < typename LocalType>
+template <typename MatrixScalar, typename VectorScalar>
 inline
-LocalType
+typename ProductType<VectorScalar,MatrixScalar>::type
 ConstraintMatrix::
 resolve_vector_entry (const size_type                       i,
                       const internals::GlobalRowsFromLocal &global_rows,
-                      const Vector<LocalType>              &local_vector,
+                      const Vector<VectorScalar>           &local_vector,
                       const std::vector<size_type>         &local_dof_indices,
-                      const FullMatrix<LocalType>          &local_matrix) const
+                      const FullMatrix<MatrixScalar>       &local_matrix) const
 {
   const size_type loc_row = global_rows.local_row(i);
   const size_type n_inhomogeneous_rows = global_rows.n_inhomogeneities();
-  LocalType val = 0;
+  typename ProductType<VectorScalar,MatrixScalar>::type val = 0;
   // has a direct contribution from some local entry. If we have inhomogeneous
   // constraints, compute the contribution of the inhomogeneity in the current
   // row.
@@ -2217,7 +2217,7 @@ resolve_vector_entry (const size_type                       i,
   for (size_type q=0; q<global_rows.size(i); ++q)
     {
       const size_type loc_row_q = global_rows.local_row(i,q);
-      LocalType add_this = local_vector (loc_row_q);
+      typename ProductType<VectorScalar,MatrixScalar>::type add_this = local_vector (loc_row_q);
       for (size_type k=0; k<n_inhomogeneous_rows; ++k)
         add_this -= (local_matrix(loc_row_q,global_rows.constraint_origin(k)) *
                      lines[lines_cache[calculate_line_index
@@ -2263,7 +2263,7 @@ ConstraintMatrix::distribute_local_to_global (
 
   const size_type n_local_dofs = local_dof_indices.size();
 
-  typename internals::ConstraintMatrixData<number>::ScratchDataAccessor
+  typename internals::ConstraintMatrixData<typename MatrixType::value_type,typename VectorType::value_type>::ScratchDataAccessor
   scratch_data;
 
   internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
@@ -2282,7 +2282,7 @@ ConstraintMatrix::distribute_local_to_global (
   std::vector<number>    &vals           = scratch_data->values;
   // create arrays for writing into the vector as well
   std::vector<size_type> &vector_indices = scratch_data->vector_indices;
-  std::vector<number>    &vector_values  = scratch_data->vector_values;
+  std::vector<typename VectorType::value_type> &vector_values  = scratch_data->vector_values;
   vector_indices.resize(n_actual_dofs);
   vector_values.resize(n_actual_dofs);
   SparseMatrix<number> *sparse_matrix
@@ -2328,13 +2328,14 @@ ConstraintMatrix::distribute_local_to_global (
       // hand side.
       if (use_vectors == true)
         {
-          const number val = resolve_vector_entry (i, global_rows,
-                                                   local_vector,
-                                                   local_dof_indices,
-                                                   local_matrix);
+          const typename VectorType::value_type
+          val = resolve_vector_entry (i, global_rows,
+                                      local_vector,
+                                      local_dof_indices,
+                                      local_matrix);
           AssertIsFinite(val);
 
-          if (val != number ())
+          if (val != typename VectorType::value_type ())
             {
               vector_indices[local_row_n] = row;
               vector_values[local_row_n] = val;
@@ -2406,7 +2407,7 @@ distribute_local_to_global (
     }
   Assert (sorted == true, ExcMatrixNotClosed());
 
-  typename internals::ConstraintMatrixData<number>::ScratchDataAccessor
+  typename internals::ConstraintMatrixData<typename MatrixType::value_type,typename VectorType::value_type>::ScratchDataAccessor
   scratch_data;
 
   const size_type n_local_dofs = local_dof_indices.size();
@@ -2528,7 +2529,7 @@ ConstraintMatrix::distribute_local_to_global (
   const size_type n_local_row_dofs = row_indices.size();
   const size_type n_local_col_dofs = col_indices.size();
 
-  typename internals::ConstraintMatrixData<number>::ScratchDataAccessor
+  typename internals::ConstraintMatrixData<typename MatrixType::value_type>::ScratchDataAccessor
   scratch_data;
   internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
   global_rows.reinit(n_local_row_dofs);
diff --git a/include/deal.II/numerics/matrix_creator.templates.h b/include/deal.II/numerics/matrix_creator.templates.h
index 578a8c054b..39b48200ba 100644
--- a/include/deal.II/numerics/matrix_creator.templates.h
+++ b/include/deal.II/numerics/matrix_creator.templates.h
@@ -142,7 +142,7 @@ namespace MatrixCreator
       {
         Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
                  const UpdateFlags         update_flags,
-                 const Function<spacedim,number> *coefficient,
+                 const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *coefficient,
                  const Function<spacedim,number> *rhs_function,
                  const ::dealii::hp::QCollection<dim> &quadrature,
                  const ::dealii::hp::MappingCollection<dim,spacedim> &mapping)
@@ -156,7 +156,7 @@ namespace MatrixCreator
                        update_flags),
           coefficient_values(quadrature_collection.max_n_quadrature_points()),
           coefficient_vector_values (quadrature_collection.max_n_quadrature_points(),
-                                     dealii::Vector<number> (fe_collection.n_components())),
+                                     dealii::Vector<typename numbers::NumberTraits<number>::real_type> (fe_collection.n_components())),
           rhs_values(quadrature_collection.max_n_quadrature_points()),
           rhs_vector_values(quadrature_collection.max_n_quadrature_points(),
                             dealii::Vector<number> (fe_collection.n_components())),
@@ -196,13 +196,13 @@ namespace MatrixCreator
 
         ::dealii::hp::FEValues<dim,spacedim> x_fe_values;
 
-        std::vector<number>                  coefficient_values;
-        std::vector<dealii::Vector<number> > coefficient_vector_values;
-        std::vector<number>                  rhs_values;
-        std::vector<dealii::Vector<number> > rhs_vector_values;
+        std::vector<typename numbers::NumberTraits<number>::real_type>                  coefficient_values;
+        std::vector<dealii::Vector<typename numbers::NumberTraits<number>::real_type> > coefficient_vector_values;
+        std::vector<number>                                                             rhs_values;
+        std::vector<dealii::Vector<number> >                                            rhs_vector_values;
 
-        const Function<spacedim,number>   *coefficient;
-        const Function<spacedim,number>   *rhs_function;
+        const Function<spacedim,typename numbers::NumberTraits<number>::real_type>     *coefficient;
+        const Function<spacedim,number>                                                *rhs_function;
 
         const UpdateFlags update_flags;
       };
@@ -211,10 +211,10 @@ namespace MatrixCreator
       template <typename number>
       struct CopyData
       {
-        std::vector<types::global_dof_index> dof_indices;
-        FullMatrix<number>        cell_matrix;
-        dealii::Vector<number>    cell_rhs;
-        const ConstraintMatrix   *constraints;
+        std::vector<types::global_dof_index>                          dof_indices;
+        FullMatrix<typename numbers::NumberTraits<number>::real_type> cell_matrix;
+        dealii::Vector<number>                                        cell_rhs;
+        const ConstraintMatrix                                       *constraints;
       };
     }
 
@@ -225,7 +225,7 @@ namespace MatrixCreator
               typename number>
     void mass_assembler (const CellIterator &cell,
                          MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim,number> &data,
-                         MatrixCreator::internal::AssemblerData::CopyData<number>      &copy_data)
+                         MatrixCreator::internal::AssemblerData::CopyData<number>             &copy_data)
     {
       data.x_fe_values.reinit (cell);
       const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
@@ -280,14 +280,13 @@ namespace MatrixCreator
           else
             {
               data.coefficient_vector_values.resize (n_q_points,
-                                                     dealii::Vector<number>(n_components));
+                                                     dealii::Vector<typename numbers::NumberTraits<number>::real_type>(n_components));
               data.coefficient->vector_value_list (fe_values.get_quadrature_points(),
                                                    data.coefficient_vector_values);
             }
         }
 
 
-      number add_data;
       const std::vector<double> &JxW = fe_values.get_JxW_values();
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         if (fe.is_primitive ())
@@ -295,7 +294,6 @@ namespace MatrixCreator
             const unsigned int component_i =
               fe.system_to_component_index(i).first;
             const double *phi_i = &fe_values.shape_value(i,0);
-            add_data = 0;
 
             // use symmetry in the mass matrix here:
             // just need to calculate the diagonal
@@ -307,7 +305,7 @@ namespace MatrixCreator
                    component_i))
                 {
                   const double *phi_j = &fe_values.shape_value(j,0);
-                  add_data = 0;
+                  typename numbers::NumberTraits<number>::real_type add_data = 0;
                   if (use_coefficient)
                     {
                       if (data.coefficient->n_components==1)
@@ -331,7 +329,7 @@ namespace MatrixCreator
 
             if (use_rhs_function)
               {
-                add_data = 0;
+                number add_data = 0;
                 if (data.rhs_function->n_components==1)
                   for (unsigned int point=0; point<n_q_points; ++point)
                     add_data += data.rhs_values[point] *
@@ -349,7 +347,7 @@ namespace MatrixCreator
             // symmetry again
             for (unsigned int j=i; j<dofs_per_cell; ++j)
               {
-                add_data = 0;
+                typename numbers::NumberTraits<number>::real_type add_data = 0;
                 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
                   if (fe.get_nonzero_components(i)[comp_i] &&
                       fe.get_nonzero_components(j)[comp_i])
@@ -381,7 +379,7 @@ namespace MatrixCreator
 
             if (use_rhs_function)
               {
-                add_data = 0;
+                number add_data = 0;
                 for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
                   if (fe.get_nonzero_components(i)[comp_i])
                     {
@@ -408,7 +406,7 @@ namespace MatrixCreator
               typename CellIterator>
     void laplace_assembler (const CellIterator &cell,
                             MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim,double> &data,
-                            MatrixCreator::internal::AssemblerData::CopyData<double>      &copy_data)
+                            MatrixCreator::internal::AssemblerData::CopyData<double>             &copy_data)
     {
       data.x_fe_values.reinit (cell);
       const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
@@ -630,12 +628,12 @@ namespace MatrixCreator
 
         CopyData(CopyData const &data);
 
-        unsigned int dofs_per_cell;
-        std::vector<types::global_dof_index> dofs;
-        std::vector<std::vector<bool> > dof_is_on_face;
-        typename DoFHandlerType::active_cell_iterator cell;
-        std::vector<FullMatrix<number> > cell_matrix;
-        std::vector<Vector<number> > cell_vector;
+        unsigned int                                                                dofs_per_cell;
+        std::vector<types::global_dof_index>                                        dofs;
+        std::vector<std::vector<bool> >                                             dof_is_on_face;
+        typename DoFHandlerType::active_cell_iterator                               cell;
+        std::vector<FullMatrix<typename numbers::NumberTraits<number>::real_type> > cell_matrix;
+        std::vector<Vector<number> >                                                cell_vector;
       };
 
 
@@ -723,10 +721,10 @@ namespace MatrixCreator
   void create_mass_matrix (const Mapping<dim,spacedim>       &mapping,
                            const DoFHandler<dim,spacedim>    &dof,
                            const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number>      &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const coefficient,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const coefficient,
                            const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
@@ -754,7 +752,7 @@ namespace MatrixCreator
                      static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
                      &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>,
                      std::bind(&MatrixCreator::internal::
-                               copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
+                               copy_local_to_global<number,SparseMatrix<typename numbers::NumberTraits<number>::real_type>, Vector<number> >,
                                std::placeholders::_1, &matrix, &rhs_vector),
                      assembler_data,
                      copy_data);
@@ -765,10 +763,10 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
                            const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number>      &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const coefficient,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const coefficient,
                            const ConstraintMatrix   &constraints)
   {
     create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping,
@@ -833,10 +831,10 @@ namespace MatrixCreator
   void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
                            const hp::DoFHandler<dim,spacedim> &dof,
                            const hp::QCollection<dim>    &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number>      &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const coefficient,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const coefficient,
                            const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
@@ -861,7 +859,7 @@ namespace MatrixCreator
                      static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
                      &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>,
                      std::bind (&MatrixCreator::internal::
-                                copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
+                                copy_local_to_global<number,SparseMatrix<typename numbers::NumberTraits<number>::real_type>, Vector<number> >,
                                 std::placeholders::_1, &matrix, &rhs_vector),
                      assembler_data,
                      copy_data);
@@ -872,10 +870,10 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
                            const hp::QCollection<dim> &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number> &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const coefficient,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const coefficient,
                            const ConstraintMatrix   &constraints)
   {
     create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
@@ -892,12 +890,12 @@ namespace MatrixCreator
     create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell,
                                    MatrixCreator::internal::AssemblerBoundary::Scratch const &,
                                    MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
-                                   spacedim>,number > &copy_data,
+                                   spacedim>,number> &copy_data,
                                    Mapping<dim, spacedim> const &mapping,
                                    FiniteElement<dim,spacedim> const &fe,
                                    Quadrature<dim-1> const &q,
                                    std::map<types::boundary_id, const Function<spacedim,number>*> const &boundary_functions,
-                                   Function<spacedim,number> const *const coefficient,
+                                   Function<spacedim,typename numbers::NumberTraits<number>::real_type> const *const coefficient,
                                    std::vector<unsigned int> const &component_mapping)
 
     {
@@ -921,9 +919,9 @@ namespace MatrixCreator
 
       // two variables for the coefficient, one for the two cases
       // indicated in the name
-      std::vector<number>          coefficient_values (fe_values.n_quadrature_points, 1.);
-      std::vector<Vector<number> > coefficient_vector_values (fe_values.n_quadrature_points,
-                                                              Vector<number>(n_components));
+      std::vector<typename numbers::NumberTraits<number>::real_type>          coefficient_values (fe_values.n_quadrature_points, 1.);
+      std::vector<Vector<typename numbers::NumberTraits<number>::real_type> > coefficient_vector_values (fe_values.n_quadrature_points,
+          Vector<typename numbers::NumberTraits<number>::real_type>(n_components));
       const bool coefficient_is_vector = (coefficient != nullptr && coefficient->n_components != 1);
 
       std::vector<number>          rhs_values_scalar (fe_values.n_quadrature_points);
@@ -1096,7 +1094,7 @@ namespace MatrixCreator
                                      spacedim>,number > const &copy_data,
                                      std::map<types::boundary_id, const Function<spacedim,number>*> const &boundary_functions,
                                      std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
-                                     SparseMatrix<number> &matrix,
+                                     SparseMatrix<typename numbers::NumberTraits<number>::real_type> &matrix,
                                      Vector<number> &rhs_vector)
     {
       // now transfer cell matrix and vector to the whole boundary matrix
@@ -1204,11 +1202,11 @@ namespace MatrixCreator
   create_boundary_mass_matrix (const Mapping<dim, spacedim>  &mapping,
                                const DoFHandler<dim,spacedim> &dof,
                                const Quadrature<dim-1>  &q,
-                               SparseMatrix<number>  &matrix,
+                               SparseMatrix<typename numbers::NumberTraits<number>::real_type>  &matrix,
                                const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
                                Vector<number>            &rhs_vector,
                                std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                               const Function<spacedim,number> *const coefficient,
+                               const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const coefficient,
                                std::vector<unsigned int> component_mapping)
   {
     // what would that be in 1d? the
@@ -1283,7 +1281,7 @@ namespace MatrixCreator
                                       hp::FECollection<dim,spacedim> const &fe_collection,
                                       hp::QCollection<dim-1> const &q,
                                       const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
-                                      Function<spacedim,number> const *const coefficient,
+                                      Function<spacedim,typename numbers::NumberTraits<number>::real_type> const *const coefficient,
                                       std::vector<unsigned int> const &component_mapping)
     {
       const unsigned int n_components  = fe_collection.n_components();
@@ -1308,8 +1306,8 @@ namespace MatrixCreator
       // two variables for the coefficient,
       // one for the two cases indicated in
       // the name
-      std::vector<number>          coefficient_values;
-      std::vector<Vector<number> > coefficient_vector_values;
+      std::vector<typename numbers::NumberTraits<number>::real_type>          coefficient_values;
+      std::vector<Vector<typename numbers::NumberTraits<number>::real_type> > coefficient_vector_values;
 
       const bool coefficient_is_vector = (coefficient != nullptr && coefficient->n_components != 1);
 
@@ -1347,7 +1345,7 @@ namespace MatrixCreator
               // FE has several components
               {
                 coefficient_vector_values.resize (fe_values.n_quadrature_points,
-                                                  Vector<number>(n_components));
+                                                  Vector<typename numbers::NumberTraits<number>::real_type>(n_components));
                 coefficient_values.resize (fe_values.n_quadrature_points, 1.);
 
                 rhs_values_system.resize (fe_values.n_quadrature_points,
@@ -1498,7 +1496,7 @@ namespace MatrixCreator
                                         ::CopyData<hp::DoFHandler<dim,spacedim>,number > const &copy_data,
                                         std::map<types::boundary_id, const Function<spacedim,number>*> const &boundary_functions,
                                         std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
-                                        SparseMatrix<number> &matrix,
+                                        SparseMatrix<typename numbers::NumberTraits<number>::real_type> &matrix,
                                         Vector<number> &rhs_vector)
     {
       // now transfer cell matrix and vector to the whole boundary matrix
@@ -1591,11 +1589,11 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_boundary_mass_matrix (const DoFHandler<dim,spacedim>     &dof,
                                     const Quadrature<dim-1>   &q,
-                                    SparseMatrix<number>      &matrix,
+                                    SparseMatrix<typename numbers::NumberTraits<number>::real_type>      &matrix,
                                     const std::map<types::boundary_id, const Function<spacedim,number>*> &rhs,
                                     Vector<number>            &rhs_vector,
                                     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim,number> *const a,
+                                    const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a,
                                     std::vector<unsigned int> component_mapping)
   {
     create_boundary_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
@@ -1609,11 +1607,11 @@ namespace MatrixCreator
   create_boundary_mass_matrix (const hp::MappingCollection<dim,spacedim>        &mapping,
                                const hp::DoFHandler<dim,spacedim>     &dof,
                                const hp::QCollection<dim-1>   &q,
-                               SparseMatrix<number>      &matrix,
+                               SparseMatrix<typename numbers::NumberTraits<number>::real_type>      &matrix,
                                const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
                                Vector<number>            &rhs_vector,
                                std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                               const Function<spacedim,number> *const coefficient,
+                               const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const coefficient,
                                std::vector<unsigned int> component_mapping)
   {
     // what would that be in 1d? the
@@ -1678,11 +1676,11 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_boundary_mass_matrix (const hp::DoFHandler<dim,spacedim>     &dof,
                                     const hp::QCollection<dim-1>   &q,
-                                    SparseMatrix<number>      &matrix,
+                                    SparseMatrix<typename numbers::NumberTraits<number>::real_type>      &matrix,
                                     const std::map<types::boundary_id, const Function<spacedim,number>*> &rhs,
                                     Vector<number>            &rhs_vector,
                                     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim,number> *const a,
+                                    const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a,
                                     std::vector<unsigned int> component_mapping)
   {
     create_boundary_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
diff --git a/include/deal.II/numerics/matrix_tools.h b/include/deal.II/numerics/matrix_tools.h
index 06d6ff3bc0..b3fe2c22e9 100644
--- a/include/deal.II/numerics/matrix_tools.h
+++ b/include/deal.II/numerics/matrix_tools.h
@@ -277,10 +277,10 @@ namespace MatrixCreator
   void create_mass_matrix (const Mapping<dim, spacedim>   &mapping,
                            const DoFHandler<dim,spacedim> &dof,
                            const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number> &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const a = nullptr,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                            const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
@@ -290,10 +290,10 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_mass_matrix (const DoFHandler<dim,spacedim> &dof,
                            const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number> &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const a = nullptr,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                            const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
@@ -324,10 +324,10 @@ namespace MatrixCreator
   void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
                            const hp::DoFHandler<dim,spacedim> &dof,
                            const hp::QCollection<dim> &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number> &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const a = nullptr,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                            const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
@@ -336,10 +336,10 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
                            const hp::QCollection<dim> &q,
-                           SparseMatrix<number>     &matrix,
+                           SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                            const Function<spacedim,number> &rhs,
                            Vector<number>           &rhs_vector,
-                           const Function<spacedim,number> *const a = nullptr,
+                           const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                            const ConstraintMatrix   &constraints = ConstraintMatrix());
 
 
@@ -372,11 +372,11 @@ namespace MatrixCreator
   void create_boundary_mass_matrix (const Mapping<dim, spacedim>       &mapping,
                                     const DoFHandler<dim,spacedim>    &dof,
                                     const Quadrature<dim-1>  &q,
-                                    SparseMatrix<number>     &matrix,
+                                    SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                                     const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
                                     Vector<number>           &rhs_vector,
                                     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim,number> *const weight = 0,
+                                    const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const weight = 0,
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
 
@@ -387,11 +387,11 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_boundary_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
                                     const Quadrature<dim-1>  &q,
-                                    SparseMatrix<number>     &matrix,
+                                    SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                                     const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
                                     Vector<number>           &rhs_vector,
                                     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim,number> *const a = nullptr,
+                                    const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
   /**
@@ -401,11 +401,11 @@ namespace MatrixCreator
   void create_boundary_mass_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
                                     const hp::DoFHandler<dim,spacedim>    &dof,
                                     const hp::QCollection<dim-1>  &q,
-                                    SparseMatrix<number>     &matrix,
+                                    SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                                     const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
                                     Vector<number>           &rhs_vector,
                                     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim,number> *const a = nullptr,
+                                    const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
   /**
@@ -414,11 +414,11 @@ namespace MatrixCreator
   template <int dim, int spacedim, typename number>
   void create_boundary_mass_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
                                     const hp::QCollection<dim-1>  &q,
-                                    SparseMatrix<number>     &matrix,
+                                    SparseMatrix<typename numbers::NumberTraits<number>::real_type>     &matrix,
                                     const std::map<types::boundary_id, const Function<spacedim,number>*> &boundary_functions,
                                     Vector<number>           &rhs_vector,
                                     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim,number> *const a = nullptr,
+                                    const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *const a = nullptr,
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
   /**
diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 65b671e924..a134b95429 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -22,6 +22,7 @@
 #include <deal.II/base/polynomials_piecewise.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/quadrature.h>
+#include <deal.II/base/numbers.h>
 #include <deal.II/distributed/tria_base.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
@@ -37,6 +38,7 @@
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
 #include <deal.II/lac/vector_memory.h>
 #include <deal.II/lac/filtered_matrix.h>
 #include <deal.II/lac/constraint_matrix.h>
@@ -922,7 +924,7 @@ namespace VectorTools
     }
 
     template <typename number>
-    void invert_mass_matrix(const SparseMatrix<std::complex<number> > &/*mass_matrix*/,
+    void invert_mass_matrix(const SparseMatrix<number> &/*mass_matrix*/,
                             const Vector<std::complex<number> > &/*rhs*/,
                             Vector<std::complex<number> > &/*solution*/)
     {
@@ -2854,9 +2856,25 @@ namespace VectorTools
       return a > b;
     }
 
+    template <typename number1, typename number2>
+    bool real_part_bigger_than(const number1               a,
+                               const std::complex<number2> b)
+    {
+      Assert(std::abs(b.imag()) <= 1e-15*std::abs(b), ExcInternalError());
+      return a > b.real();
+    }
+
+    template <typename number1, typename number2>
+    bool real_part_bigger_than(const std::complex<number1> a,
+                               const number2               b)
+    {
+      Assert(std::abs(a.imag()) <= 1e-15*std::abs(a), ExcInternalError());
+      return a.real() > b;
+    }
+
     template <typename number1, typename number2>
     bool real_part_bigger_than(const std::complex<number1> a,
-                               const std::complex<number2> &b)
+                               const std::complex<number2> b)
     {
       Assert(std::abs(a.imag()) <= 1e-15*std::abs(a), ExcInternalError());
       Assert(std::abs(b.imag()) <= 1e-15*std::abs(b), ExcInternalError());
@@ -2906,13 +2924,13 @@ namespace VectorTools
 
     template <typename number>
     void invert_mass_matrix
-    (const SparseMatrix<std::complex<number> > &/*mass_matrix*/,
-     const FilteredMatrix<Vector<std::complex<number> > > &/*filtered_mass_matrix*/,
-     FilteredMatrix<Vector<std::complex<number> > > &/*filtered_preconditioner*/,
-     const Vector<std::complex<number> > &/*rhs*/,
-     Vector<std::complex<number> > &/*boundary_projection*/)
+    (const SparseMatrix<number> &mass_matrix,
+     const FilteredMatrix<Vector<std::complex<number> > > &filtered_mass_matrix,
+     FilteredMatrix<Vector<std::complex<number> > > &filtered_preconditioner,
+     const Vector<std::complex<number> > &rhs,
+     Vector<std::complex<number> > &boundary_projection)
     {
-      Assert(false, ExcNotImplemented());
+      Assert (false, ExcInternalError());
     }
 
 
@@ -3026,14 +3044,14 @@ namespace VectorTools
 
 
       // make mass matrix and right hand side
-      SparseMatrix<number> mass_matrix(sparsity);
+      SparseMatrix<typename numbers::NumberTraits<number>::real_type> mass_matrix(sparsity);
       Vector<number>       rhs(sparsity.n_rows());
 
 
       MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
                                                   mass_matrix, boundary_functions,
                                                   rhs, dof_to_boundary_mapping,
-                                                  (const Function<spacedim,number> *) nullptr,
+                                                  (const Function<spacedim,typename numbers::NumberTraits<number>::real_type> *) nullptr,
                                                   component_mapping);
 
       // For certain weird elements,
diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc
index 6b159a292f..5939293efa 100644
--- a/source/lac/constraint_matrix.cc
+++ b/source/lac/constraint_matrix.cc
@@ -1324,13 +1324,13 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
                                                       VectorType                      &, \
                                                       bool                             , \
                                                       std::integral_constant<bool, false>) const
-#define MATRIX_FUNCTIONS(MatrixType) \
+#define MATRIX_FUNCTIONS(MatrixType,VectorScalar) \
   template void ConstraintMatrix:: \
-  distribute_local_to_global<MatrixType,Vector<MatrixType::value_type> > (const FullMatrix<MatrixType::value_type>        &, \
-      const Vector<MatrixType::value_type>            &, \
+  distribute_local_to_global<MatrixType,Vector<VectorScalar> > (const FullMatrix<MatrixType::value_type>        &, \
+      const Vector<VectorScalar>            &, \
       const std::vector<ConstraintMatrix::size_type> &, \
       MatrixType                      &, \
-      Vector<MatrixType::value_type>                  &, \
+      Vector<VectorScalar>                  &, \
       bool                             , \
       std::integral_constant<bool, false>) const
 #define BLOCK_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType)   \
@@ -1352,37 +1352,42 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
       bool                             , \
       std::integral_constant<bool, true>) const
 
-MATRIX_FUNCTIONS(SparseMatrix<double>);
-MATRIX_FUNCTIONS(SparseMatrix<float>);
-MATRIX_FUNCTIONS(FullMatrix<double>);
-MATRIX_FUNCTIONS(FullMatrix<float>);
-MATRIX_FUNCTIONS(FullMatrix<std::complex<double> >);
-MATRIX_FUNCTIONS(SparseMatrix<std::complex<double> >);
-MATRIX_FUNCTIONS(SparseMatrix<std::complex<float> >);
+MATRIX_FUNCTIONS(FullMatrix<double>,double);
+MATRIX_FUNCTIONS(FullMatrix<float>,float);
+MATRIX_FUNCTIONS(FullMatrix<double>,std::complex<double>);
+MATRIX_FUNCTIONS(FullMatrix<std::complex<double> >,std::complex<double>);
+
+MATRIX_FUNCTIONS(SparseMatrix<double>,double);
+MATRIX_FUNCTIONS(SparseMatrix<float>,float);
+MATRIX_FUNCTIONS(SparseMatrix<double>,std::complex<double>);
+MATRIX_FUNCTIONS(SparseMatrix<float>,std::complex<float>);
+MATRIX_FUNCTIONS(SparseMatrix<std::complex<double> >,std::complex<double>);
+MATRIX_FUNCTIONS(SparseMatrix<std::complex<float> >,std::complex<float>);
+
+MATRIX_FUNCTIONS(SparseMatrixEZ<double>,double);
+MATRIX_FUNCTIONS(SparseMatrixEZ<float>,float);
+MATRIX_FUNCTIONS(ChunkSparseMatrix<double>,double);
+MATRIX_FUNCTIONS(ChunkSparseMatrix<float>,float);
+
 
 BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix<double>);
 BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix<float>);
 BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<double>, BlockVector<double>);
 BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<float>,  BlockVector<float>);
 
-MATRIX_FUNCTIONS(SparseMatrixEZ<double>);
-MATRIX_FUNCTIONS(SparseMatrixEZ<float>);
-MATRIX_FUNCTIONS(ChunkSparseMatrix<double>);
-MATRIX_FUNCTIONS(ChunkSparseMatrix<float>);
-
 // BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrixEZ<double>);
 // BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ<float>,  Vector<float>);
 
 #ifdef DEAL_II_WITH_PETSC
-MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix);
-MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix);
+MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix,PetscScalar);
+MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix,PetscScalar);
 BLOCK_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix);
 MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector);
 BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix,PETScWrappers::MPI::BlockVector);
 #endif
 
 #ifdef DEAL_II_WITH_TRILINOS
-MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix);
+MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix,double);
 BLOCK_MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix);
 MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::SparseMatrix, TrilinosWrappers::MPI::Vector);
 BLOCK_MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix, TrilinosWrappers::MPI::BlockVector);
@@ -1468,15 +1473,17 @@ ONLY_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix);
 // constructor of scratch_data (it won't allow one to be constructed in place).
 namespace internals
 {
-#define SCRATCH_INITIALIZER(Number,Name)                                \
-  ConstraintMatrixData<Number>::ScratchData scratch_data_initializer_##Name; \
-  template <> Threads::ThreadLocalStorage<ConstraintMatrixData<Number>::ScratchData> \
-  ConstraintMatrixData<Number>::scratch_data(scratch_data_initializer_##Name)
-
-  SCRATCH_INITIALIZER(double,double);
-  SCRATCH_INITIALIZER(float,float);
-  SCRATCH_INITIALIZER(std::complex<double>,cdouble);
-  SCRATCH_INITIALIZER(std::complex<float>,cfloat);
+#define SCRATCH_INITIALIZER(MatrixScalar,VectorScalar,Name)                  \
+  ConstraintMatrixData<MatrixScalar,VectorScalar>::ScratchData scratch_data_initializer_##Name; \
+  template <> Threads::ThreadLocalStorage<ConstraintMatrixData<MatrixScalar,VectorScalar>::ScratchData> \
+  ConstraintMatrixData<MatrixScalar,VectorScalar>::scratch_data(scratch_data_initializer_##Name)
+
+  SCRATCH_INITIALIZER(double,double,dd);
+  SCRATCH_INITIALIZER(float,float,ff);
+  SCRATCH_INITIALIZER(std::complex<double>,std::complex<double>,zz);
+  SCRATCH_INITIALIZER(std::complex<float>,std::complex<float>,cc);
+  SCRATCH_INITIALIZER(double,std::complex<double>,dz);
+  SCRATCH_INITIALIZER(float,std::complex<float>,fc);
 #undef SCRATCH_INITIALIZER
 }
 
diff --git a/source/lac/constraint_matrix.inst.in b/source/lac/constraint_matrix.inst.in
index 9244eed2fb..b3c08c4829 100644
--- a/source/lac/constraint_matrix.inst.in
+++ b/source/lac/constraint_matrix.inst.in
@@ -29,21 +29,57 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES)
 {
     template void ConstraintMatrix::condense<LinearAlgebra::distributed::T<S> >(const LinearAlgebra::distributed::T<S> &, LinearAlgebra::distributed::T<S> &) const;
     template void ConstraintMatrix::condense<LinearAlgebra::distributed::T<S> >(LinearAlgebra::distributed::T<S> &vec) const;
-    template void ConstraintMatrix::distribute_local_to_global<LinearAlgebra::distributed::T<S> > (
-        const Vector<S>&, const std::vector<types::global_dof_index> &, LinearAlgebra::distributed::T<S> &, const FullMatrix<S>&) const;
-    template void ConstraintMatrix::distribute_local_to_global<LinearAlgebra::distributed::T<S> > (
-        const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, LinearAlgebra::distributed::T<S> &, const FullMatrix<S>&, bool) const;
-    template void ConstraintMatrix::set_zero<LinearAlgebra::distributed::T<S> >(LinearAlgebra::distributed::T<S> &) const;
-    template void ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> > > (
-        const FullMatrix<S> &, const std::vector< size_type > &, DiagonalMatrix<LinearAlgebra::distributed::T<S> > &) const;
-    template void ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> >, LinearAlgebra::distributed::T<S> > (
-        const FullMatrix<S> &, const Vector<S>&, const std::vector< size_type > &,
-        DiagonalMatrix<LinearAlgebra::distributed::T<S> > &, LinearAlgebra::distributed::T<S>&,
-        bool, std::integral_constant<bool, false>) const;
-    template void ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> >, T<S> > (
-        const FullMatrix<S> &, const Vector<S>&, const std::vector< size_type > &,
-        DiagonalMatrix<LinearAlgebra::distributed::T<S> > &, T<S>&,
-        bool, std::integral_constant<bool, false>) const;
+
+    template
+    void
+    ConstraintMatrix::distribute_local_to_global<LinearAlgebra::distributed::T<S> >
+    (const Vector<S>&,
+     const std::vector<types::global_dof_index> &,
+     LinearAlgebra::distributed::T<S> &,
+     const FullMatrix<S>&) const;
+
+    template
+    void
+    ConstraintMatrix::distribute_local_to_global<LinearAlgebra::distributed::T<S> >
+    (const Vector<S>&,
+     const std::vector<types::global_dof_index> &,
+     const std::vector<types::global_dof_index> &,
+     LinearAlgebra::distributed::T<S> &,
+     const FullMatrix<S>&,
+     bool) const;
+
+    template
+    void
+    ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> > >
+    (const FullMatrix<S> &,
+     const std::vector< size_type > &,
+     DiagonalMatrix<LinearAlgebra::distributed::T<S> > &) const;
+
+    template
+    void
+    ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> >, LinearAlgebra::distributed::T<S> >
+    (const FullMatrix<S> &,
+     const Vector<S>&,
+     const std::vector< size_type > &,
+     DiagonalMatrix<LinearAlgebra::distributed::T<S> > &,
+     LinearAlgebra::distributed::T<S>&,
+     bool,
+     std::integral_constant<bool, false>) const;
+
+    template
+    void
+    ConstraintMatrix::distribute_local_to_global<DiagonalMatrix<LinearAlgebra::distributed::T<S> >, T<S> >
+    (const FullMatrix<S> &,
+     const Vector<S>&,
+     const std::vector< size_type > &,
+     DiagonalMatrix<LinearAlgebra::distributed::T<S> > &,
+     T<S>&,
+     bool,
+     std::integral_constant<bool, false>) const;
+
+    template
+    void
+    ConstraintMatrix::set_zero<LinearAlgebra::distributed::T<S> >(LinearAlgebra::distributed::T<S> &) const;
 }
 
 
diff --git a/source/lac/sparse_matrix.inst.in b/source/lac/sparse_matrix.inst.in
index 1a54ee2587..bdb65baec8 100644
--- a/source/lac/sparse_matrix.inst.in
+++ b/source/lac/sparse_matrix.inst.in
@@ -137,6 +137,19 @@ for (S1, S2, S3 : REAL_SCALARS;
     Tvmult_add (V1<S2> &, const V2<S3> &) const;
 }
 
+for (S1 : REAL_SCALARS; S2, S3 : COMPLEX_SCALARS;
+        V1, V2     : DEAL_II_VEC_TEMPLATES)
+{
+    template void SparseMatrix<S1>::
+    vmult (V1<S2> &, const V2<S3> &) const;
+    template void SparseMatrix<S1>::
+    Tvmult (V1<S2> &, const V2<S3> &) const;
+    template void SparseMatrix<S1>::
+    vmult_add (V1<S2> &, const V2<S3> &) const;
+    template void SparseMatrix<S1>::
+    Tvmult_add (V1<S2> &, const V2<S3> &) const;
+}
+
 for (S1 : REAL_SCALARS)
 {
     template void SparseMatrix<S1>::
diff --git a/source/numerics/matrix_creator.inst.in b/source/numerics/matrix_creator.inst.in
index 2ddbd7e66d..e5cb415035 100644
--- a/source/numerics/matrix_creator.inst.in
+++ b/source/numerics/matrix_creator.inst.in
@@ -14,6 +14,8 @@
 // ---------------------------------------------------------------------
 
 
+
+// MatrixCreator::create_mass_matrix for the matrix only --> only real-valued scalars
 for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
 {
 #if deal_II_dimension <= deal_II_space_dimension
@@ -33,29 +35,30 @@ for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimensi
      SparseMatrix<scalar>     &matrix,
      const Function<deal_II_space_dimension,scalar> * const coefficient,
      const ConstraintMatrix   &constraints);
+
+    // hp versions of functions
     template
     void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const Quadrature<deal_II_dimension>    &q,
+    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
+     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+     const hp::QCollection<deal_II_dimension>    &q,
      SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar>      &rhs,
-     Vector<scalar>           &rhs_vector,
      const Function<deal_II_space_dimension,scalar> * const coefficient,
      const ConstraintMatrix   &constraints);
+
     template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const Quadrature<deal_II_dimension>    &q,
+    void MatrixCreator::create_mass_matrix<deal_II_dimension, deal_II_space_dimension, scalar>
+    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+     const hp::QCollection<deal_II_dimension>    &q,
      SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar>      &rhs,
-     Vector<scalar>           &rhs_vector,
      const Function<deal_II_space_dimension,scalar> * const coefficient,
      const ConstraintMatrix   &constraints);
 #endif
 }
 
-for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+
+// MatrixCreator::create_mass_matrix for matrix + rhs --> real- and complex-valued scalars
+for (scalar: REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
 {
 #if deal_II_dimension <= deal_II_space_dimension
     // non-hp version of create_mass_matrix
@@ -64,103 +67,58 @@ for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dime
     (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
      const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
      const Quadrature<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>     &matrix,
+     const Function<deal_II_space_dimension,scalar>      &rhs,
+     Vector<scalar>           &rhs_vector,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const coefficient,
      const ConstraintMatrix   &constraints);
     template
     void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
     (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
      const Quadrature<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>     &matrix,
+     const Function<deal_II_space_dimension,scalar>      &rhs,
+     Vector<scalar>           &rhs_vector,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const coefficient,
      const ConstraintMatrix   &constraints);
+
+    // hp version
     template
     void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const Quadrature<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
+    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
+     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+     const hp::QCollection<deal_II_dimension>    &q,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>     &matrix,
      const Function<deal_II_space_dimension,scalar>      &rhs,
      Vector<scalar>           &rhs_vector,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const coefficient,
      const ConstraintMatrix   &constraints);
+
     template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const Quadrature<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
+    void MatrixCreator::create_mass_matrix<deal_II_dimension, deal_II_space_dimension, scalar>
+    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+     const hp::QCollection<deal_II_dimension>    &q,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>     &matrix,
      const Function<deal_II_space_dimension,scalar>      &rhs,
      Vector<scalar>           &rhs_vector,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const coefficient,
      const ConstraintMatrix   &constraints);
 #endif
 }
 
 
-for (scalar: REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
-{
-#if deal_II_dimension <= deal_II_space_dimension
-    template
-    void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
-     const Quadrature<deal_II_dimension-1>   &q,
-     SparseMatrix<scalar>      &matrix,
-     const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*> &rhs,
-     Vector<scalar>            &rhs_vector,
-     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-     const Function<deal_II_space_dimension,scalar> * const a,
-     std::vector<unsigned int>);
-
-    template
-    void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension, scalar>
-    (const Mapping<deal_II_dimension,deal_II_space_dimension> &,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
-     const Quadrature<deal_II_dimension-1>   &q,
-     SparseMatrix<scalar>      &matrix,
-     const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*> &rhs,
-     Vector<scalar>            &rhs_vector,
-     std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-     const Function<deal_II_space_dimension,scalar> * const a,
-     std::vector<unsigned int>);
-
-    template
-    void
-    MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&,
-     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
-     const hp::QCollection<deal_II_dimension-1>&,
-     SparseMatrix<scalar>&,
-     const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*>&,
-     Vector<scalar>&,
-     std::vector<types::global_dof_index>&,
-     const Function<deal_II_space_dimension,scalar> * const,
-     std::vector<unsigned int>);
-
-    template
-    void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
-     const hp::QCollection<deal_II_dimension-1>&,
-     SparseMatrix<scalar>&,
-     const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*>&,
-     Vector<scalar>&,
-     std::vector<types::global_dof_index>&,
-     const Function<deal_II_space_dimension,scalar> * const,
-     std::vector<unsigned int>);
-#endif
-}
-
-for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+for (scalar: REAL_AND_COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
 {
 #if deal_II_dimension <= deal_II_space_dimension
     template
     void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
     (const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
      const Quadrature<deal_II_dimension-1>   &q,
-     SparseMatrix<scalar>      &matrix,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>      &matrix,
      const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*> &rhs,
      Vector<scalar>            &rhs_vector,
      std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-     const Function<deal_II_space_dimension,scalar> * const a,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const a,
      std::vector<unsigned int>);
 
     template
@@ -168,11 +126,11 @@ for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dime
     (const Mapping<deal_II_dimension,deal_II_space_dimension> &,
      const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
      const Quadrature<deal_II_dimension-1>   &q,
-     SparseMatrix<scalar>      &matrix,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>      &matrix,
      const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*> &rhs,
      Vector<scalar>            &rhs_vector,
      std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-     const Function<deal_II_space_dimension,scalar> * const a,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const a,
      std::vector<unsigned int>);
 
     template
@@ -181,127 +139,27 @@ for (scalar: COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dime
     (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&,
      const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
      const hp::QCollection<deal_II_dimension-1>&,
-     SparseMatrix<scalar>&,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>&,
      const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*>&,
      Vector<scalar>&,
      std::vector<types::global_dof_index>&,
-     const Function<deal_II_space_dimension,scalar> * const,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const,
      std::vector<unsigned int>);
 
     template
     void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
     (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
      const hp::QCollection<deal_II_dimension-1>&,
-     SparseMatrix<scalar>&,
+     SparseMatrix<numbers::NumberTraits<scalar>::real_type>&,
      const std::map<types::boundary_id, const Function<deal_II_space_dimension,scalar>*>&,
      Vector<scalar>&,
      std::vector<types::global_dof_index>&,
-     const Function<deal_II_space_dimension,scalar> * const,
+     const Function<deal_II_space_dimension,numbers::NumberTraits<scalar>::real_type> * const,
      std::vector<unsigned int>);
 #endif
 }
 
-//TODO[SP]: replace <deal_II_dimension> by <deal_II_dimension, deal_II_space_dimension>
-// where applicable and move to codimension cases above also when applicable
-for (scalar : REAL_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
-{
-// hp versions of functions
-#if deal_II_dimension <= deal_II_space_dimension
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar>      &rhs,
-     Vector<scalar>           &rhs_vector,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-#endif
-
-#if deal_II_dimension == deal_II_space_dimension
-
-
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension, deal_II_space_dimension, scalar>
-    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension, deal_II_space_dimension, scalar>
-    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar>      &rhs,
-     Vector<scalar>           &rhs_vector,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-#endif
-}
-
-for (scalar : COMPLEX_SCALARS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
-{
-// hp versions of functions
-#if deal_II_dimension <= deal_II_space_dimension
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension,deal_II_space_dimension,scalar>
-    (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-     const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar>      &rhs,
-     Vector<scalar>           &rhs_vector,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-#endif
 
-#if deal_II_dimension == deal_II_space_dimension
-
-
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension, deal_II_space_dimension, scalar>
-    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-    template
-    void MatrixCreator::create_mass_matrix<deal_II_dimension, deal_II_space_dimension, scalar>
-    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-     const hp::QCollection<deal_II_dimension>    &q,
-     SparseMatrix<scalar>     &matrix,
-     const Function<deal_II_space_dimension,scalar>      &rhs,
-     Vector<scalar>           &rhs_vector,
-     const Function<deal_II_space_dimension,scalar> * const coefficient,
-     const ConstraintMatrix   &constraints);
-
-#endif
-}
 
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
 {
diff --git a/source/numerics/vector_tools_boundary.inst.in b/source/numerics/vector_tools_boundary.inst.in
index dc07a99915..27488e0150 100644
--- a/source/numerics/vector_tools_boundary.inst.in
+++ b/source/numerics/vector_tools_boundary.inst.in
@@ -109,27 +109,46 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 #endif
 }
 
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; number : REAL_AND_COMPLEX_SCALARS)
 {
-    namespace VectorTools \{
 #if deal_II_dimension == deal_II_space_dimension
+    namespace VectorTools \{
 
     template
     void project_boundary_values<deal_II_dimension>
     (const Mapping<deal_II_dimension>     &,
      const DoFHandler<deal_II_dimension>  &,
-     const std::map<types::boundary_id, const Function<deal_II_dimension,double>*> &,
+     const std::map<types::boundary_id, const Function<deal_II_dimension,number>*> &,
      const Quadrature<deal_II_dimension-1>&,
-     std::map<types::global_dof_index,double>&, std::vector<unsigned int>);
+     std::map<types::global_dof_index,number>&,
+     std::vector<unsigned int>);
 
     template
     void project_boundary_values<deal_II_dimension>
     (const DoFHandler<deal_II_dimension>  &,
-     const std::map<types::boundary_id, const Function<deal_II_dimension,double>*> &,
+     const std::map<types::boundary_id, const Function<deal_II_dimension,number>*> &,
      const Quadrature<deal_II_dimension-1>&,
-     std::map<types::global_dof_index,double>&,
+     std::map<types::global_dof_index,number>&,
      std::vector<unsigned int>);
 
+    template
+    void project_boundary_values<deal_II_dimension,deal_II_space_dimension>
+    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
+     const std::map<types::boundary_id, const Function<deal_II_dimension,number>*> &,
+     const hp::QCollection<deal_II_dimension-1> &,
+     std::map<types::global_dof_index,number> &,
+     std::vector<unsigned int>);
+
+    \}
+#endif
+}
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+    namespace VectorTools \{
+#if deal_II_dimension == deal_II_space_dimension
+
+
     template
     void project_boundary_values<deal_II_dimension>
     (const Mapping<deal_II_dimension>     &,
@@ -146,16 +165,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
      ConstraintMatrix&,
      std::vector<unsigned int>);
 
-    template
-    void project_boundary_values<deal_II_dimension,deal_II_space_dimension>
-    (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const std::map<types::boundary_id, const Function<deal_II_dimension,double>*> &,
-     const hp::QCollection<deal_II_dimension-1> &,
-     std::map<types::global_dof_index,double> &,
-     std::vector<unsigned int>);
-
-
-
 #if deal_II_dimension != 1
     template
     void project_boundary_values_curl_conforming<deal_II_dimension>