extend ConstraintMatrix to cover complex-valued cases

diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h
index c42b67265b3b2f3ee39a71c72490bea71f469ad7..d2549ae3b92992c10d8de94ffca4dcf1df3a45f1 100644 (file)
--- a/include/deal.II/lac/constraint_matrix.h
+++ b/include/deal.II/lac/constraint_matrix.h
@@ -717,12 +717,12 @@ public:
    * site. There is no locking mechanism inside this method to prevent data
    * races.
    */
-  template <typename VectorType>
+  template <typename VectorType, typename LocalType>
   void
-  distribute_local_to_global (const Vector<double>         &local_vector,
+  distribute_local_to_global (const Vector<LocalType>      &local_vector,
                               const std::vector<size_type> &local_dof_indices,
                               VectorType                   &global_vector,
-                              const FullMatrix<double>     &local_matrix) const;
+                              const FullMatrix<LocalType>  &local_matrix) const;
 
   /**
    * Enter a single value into a result vector, obeying constraints.
@@ -818,7 +818,7 @@ public:
    */
   template <typename MatrixType>
   void
-  distribute_local_to_global (const FullMatrix<double>     &local_matrix,
+  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type> &local_matrix,
                               const std::vector<size_type> &local_dof_indices,
                               MatrixType                   &global_matrix) const;
 
@@ -851,7 +851,7 @@ public:
    */
   template <typename MatrixType>
   void
-  distribute_local_to_global (const FullMatrix<double>     &local_matrix,
+  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type> &local_matrix,
                               const std::vector<size_type> &row_indices,
                               const std::vector<size_type> &col_indices,
                               MatrixType                   &global_matrix) const;
@@ -874,8 +874,8 @@ public:
    */
   template <typename MatrixType, typename VectorType>
   void
-  distribute_local_to_global (const FullMatrix<double>      &local_matrix,
-                              const Vector<double>          &local_vector,
+  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type> &local_matrix,
+                              const Vector<typename VectorType::value_type>     &local_vector,
                               const std::vector<size_type>  &local_dof_indices,
                               MatrixType                    &global_matrix,
                               VectorType                    &global_vector,
@@ -1239,8 +1239,8 @@ private:
    */
   template <typename MatrixType, typename VectorType>
   void
-  distribute_local_to_global (const FullMatrix<double>     &local_matrix,
-                              const Vector<double>         &local_vector,
+  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>  &local_matrix,
+                              const Vector<typename VectorType::value_type>      &local_vector,
                               const std::vector<size_type> &local_dof_indices,
                               MatrixType                   &global_matrix,
                               VectorType                   &global_vector,
@@ -1253,8 +1253,8 @@ private:
    */
   template <typename MatrixType, typename VectorType>
   void
-  distribute_local_to_global (const FullMatrix<double>     &local_matrix,
-                              const Vector<double>         &local_vector,
+  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>  &local_matrix,
+                              const Vector<typename VectorType::value_type>      &local_vector,
                               const std::vector<size_type> &local_dof_indices,
                               MatrixType                   &global_matrix,
                               VectorType                   &global_vector,
@@ -1310,12 +1310,13 @@ private:
   /**
    * Internal helper function for distribute_local_to_global function.
    */
-  double
+  template <typename LocalType>
+  LocalType
   resolve_vector_entry (const size_type                       i,
                         const internals::GlobalRowsFromLocal &global_rows,
-                        const Vector<double>                 &local_vector,
+                        const Vector<LocalType>              &local_vector,
                         const std::vector<size_type>         &local_dof_indices,
-                        const FullMatrix<double>             &local_matrix) const;
+                        const FullMatrix<LocalType>          &local_matrix) const;
 
   /**
    * The assignment operator is not implemented for performance reasons. You
@@ -1639,13 +1640,13 @@ template <typename MatrixType>
 inline
 void
 ConstraintMatrix::
-distribute_local_to_global (const FullMatrix<double>     &local_matrix,
+distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>     &local_matrix,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType                   &global_matrix) const
 {
   // create a dummy and hand on to the function actually implementing this
   // feature in the cm.templates.h file.
-  Vector<double> dummy(0);
+  Vector<typename MatrixType::value_type> dummy(0);
   distribute_local_to_global (local_matrix, dummy, local_dof_indices,
                               global_matrix, dummy, false,
                               dealii::internal::bool2type<IsBlockMatrix<MatrixType>::value>());
@@ -1653,12 +1654,13 @@ distribute_local_to_global (const FullMatrix<double>     &local_matrix,
 
 
 
+
 template <typename MatrixType, typename VectorType>
 inline
 void
 ConstraintMatrix::
-distribute_local_to_global (const FullMatrix<double>     &local_matrix,
-                            const Vector<double>         &local_vector,
+distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>     &local_matrix,
+                            const Vector<typename VectorType::value_type>         &local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType                   &global_matrix,
                             VectorType                   &global_vector,
@@ -1673,6 +1675,7 @@ distribute_local_to_global (const FullMatrix<double>     &local_matrix,
 
 
 
+
 template <typename SparsityType>
 inline
 void

diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h
index dad4ca41444b4676694ac76c60a73b25110c7b8e..c0194a7cf47ef0f58e6c1746d127036eaee748bd 100644 (file)
--- a/include/deal.II/lac/constraint_matrix.templates.h
+++ b/include/deal.II/lac/constraint_matrix.templates.h
@@ -601,13 +601,13 @@ ConstraintMatrix::set_zero (VectorType &vec) const
 
 
 
-template <typename VectorType>
+template <typename VectorType, typename LocalType>
 void
 ConstraintMatrix::
-distribute_local_to_global (const Vector<double>            &local_vector,
+distribute_local_to_global (const Vector<LocalType>      &local_vector,
                             const std::vector<size_type> &local_dof_indices,
-                            VectorType                      &global_vector,
-                            const FullMatrix<double>        &local_matrix) const
+                            VectorType                   &global_vector,
+                            const FullMatrix<LocalType>  &local_matrix) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
   AssertDimension (local_vector.size(), local_dof_indices.size());
@@ -649,9 +649,9 @@ distribute_local_to_global (const Vector<double>            &local_vector,
                   continue;
                 }
 
-              const double matrix_entry = local_matrix(j,i);
+              const LocalType matrix_entry = local_matrix(j,i);
 
-              if (matrix_entry == 0)
+              if (matrix_entry == LocalType())
                 continue;
 
               const ConstraintLine &position_j =
@@ -1535,16 +1535,17 @@ namespace internals
   // resolves constraints of one column at the innermost loop. goes through
   // the origin of each global entry and finds out which data we need to
   // collect.
+  template<typename LocalType>
   static inline
-  double resolve_matrix_entry (const GlobalRowsFromLocal &global_rows,
-                               const GlobalRowsFromLocal &global_cols,
-                               const size_type            i,
-                               const size_type            j,
-                               const size_type            loc_row,
-                               const FullMatrix<double> &local_matrix)
+  LocalType resolve_matrix_entry (const GlobalRowsFromLocal   &global_rows,
+                                  const GlobalRowsFromLocal &global_cols,
+                                  const size_type            i,
+                                  const size_type            j,
+                                  const size_type            loc_row,
+                                  const FullMatrix<LocalType> &local_matrix)
   {
     const size_type loc_col = global_cols.local_row(j);
-    double col_val;
+    LocalType col_val;
 
     // case 1: row has direct contribution in local matrix. decide whether col
     // has a direct contribution. if not, set the value to zero.
@@ -1567,8 +1568,8 @@ namespace internals
     // the direct and indirect references in the given column.
     for (size_type q=0; q<global_rows.size(i); ++q)
       {
-        double add_this = (loc_col != numbers::invalid_size_type)
-                          ? local_matrix(global_rows.local_row(i,q), loc_col) : 0;
+        LocalType add_this = (loc_col != numbers::invalid_size_type)
+                             ? local_matrix(global_rows.local_row(i,q), loc_col) : 0;
 
         for (size_type p=0; p<global_cols.size(j); ++p)
           add_this += (local_matrix(global_rows.local_row(i,q),
@@ -1585,7 +1586,7 @@ namespace internals
   // computes all entries that need to be written into global_rows[i]. Lists
   // the resulting values in val_ptr, and the corresponding column indices in
   // col_ptr.
-  template <typename number>
+  template <typename number, typename LocalType>
   inline
   void
   resolve_matrix_row (const GlobalRowsFromLocal &global_rows,
@@ -1593,7 +1594,7 @@ namespace internals
                       const size_type            i,
                       const size_type            column_start,
                       const size_type            column_end,
-                      const FullMatrix<double>  &local_matrix,
+                      const FullMatrix<LocalType> &local_matrix,
                       size_type                *&col_ptr,
                       number                   *&val_ptr)
   {
@@ -1610,14 +1611,14 @@ namespace internals
         global_cols.have_indirect_rows() == false)
       {
         AssertIndexRange(loc_row, local_matrix.m());
-        const double *matrix_ptr = &local_matrix(loc_row, 0);
+        const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
 
         for (size_type j=column_start; j<column_end; ++j)
           {
             const size_type loc_col = global_cols.local_row(j);
             AssertIndexRange(loc_col, local_matrix.n());
-            const double col_val = matrix_ptr[loc_col];
-            if (col_val != 0.)
+            const LocalType col_val = matrix_ptr[loc_col];
+            if (col_val != LocalType ())
               {
                 *val_ptr++ = static_cast<number> (col_val);
                 *col_ptr++ = global_cols.global_row(j);
@@ -1631,12 +1632,12 @@ namespace internals
       {
         for (size_type j=column_start; j<column_end; ++j)
           {
-            double col_val = resolve_matrix_entry (global_rows, global_cols, i, j,
-                                                   loc_row, local_matrix);
+            LocalType col_val = resolve_matrix_entry (global_rows, global_cols, i, j,
+                                                      loc_row, local_matrix);
 
             // if we got some nontrivial value, append it to the array of
             // values.
-            if (col_val != 0.)
+            if (col_val != LocalType ())
               {
                 *val_ptr++ = static_cast<number> (col_val);
                 *col_ptr++ = global_cols.global_row(j);
@@ -1651,15 +1652,15 @@ namespace internals
   // SparseMatrix<number>.
   namespace dealiiSparseMatrix
   {
-    template <typename SparseMatrixIterator>
+    template <typename SparseMatrixIterator, typename LocalType>
     static inline
-    void add_value (const double          value,
+    void add_value (const LocalType       value,
                     const size_type       row,
                     const size_type       column,
                     SparseMatrixIterator &matrix_values)
     {
       (void)row;
-      if (value != 0.)
+      if (value != LocalType ())
         {
           while (matrix_values->column() < column)
             ++matrix_values;
@@ -1674,14 +1675,14 @@ namespace internals
   // similar as before, now with shortcut for deal.II sparse matrices. this
   // lets us avoid using extra arrays, and does all the operations just in
   // place, i.e., in the respective matrix row
-  template <typename number>
+  template <typename number, typename LocalType>
   inline
   void
   resolve_matrix_row (const GlobalRowsFromLocal &global_rows,
                       const size_type            i,
                       const size_type            column_start,
                       const size_type            column_end,
-                      const FullMatrix<double>  &local_matrix,
+                      const FullMatrix<LocalType> &local_matrix,
                       SparseMatrix<number>      *sparse_matrix)
   {
     if (column_end == column_start)
@@ -1709,12 +1710,12 @@ namespace internals
         if (global_rows.have_indirect_rows() == false)
           {
             AssertIndexRange (loc_row, local_matrix.m());
-            const double *matrix_ptr = &local_matrix(loc_row, 0);
+            const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
 
             for (size_type j=column_start; j<column_end; ++j)
               {
                 const size_type loc_col = global_rows.local_row(j);
-                const double col_val = matrix_ptr[loc_col];
+                const LocalType col_val = matrix_ptr[loc_col];
                 dealiiSparseMatrix::add_value (col_val, row,
                                                global_rows.global_row(j),
                                                matrix_values);
@@ -1724,8 +1725,8 @@ namespace internals
           {
             for (size_type j=column_start; j<column_end; ++j)
               {
-                double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
-                                                       loc_row, local_matrix);
+                LocalType col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
+                                                          loc_row, local_matrix);
                 dealiiSparseMatrix::add_value (col_val, row,
                                                global_rows.global_row(j),
                                                matrix_values);
@@ -1738,13 +1739,13 @@ namespace internals
         if (global_rows.have_indirect_rows() == false)
           {
             AssertIndexRange (loc_row, local_matrix.m());
-            const double *matrix_ptr = &local_matrix(loc_row, 0);
+            const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
 
             sparse_matrix->begin(row)->value() += matrix_ptr[loc_row];
             for (size_type j=column_start; j<i; ++j)
               {
                 const size_type loc_col = global_rows.local_row(j);
-                const double col_val = matrix_ptr[loc_col];
+                const LocalType col_val = matrix_ptr[loc_col];
                 dealiiSparseMatrix::add_value(col_val, row,
                                               global_rows.global_row(j),
                                               matrix_values);
@@ -1752,7 +1753,7 @@ namespace internals
             for (size_type j=i+1; j<column_end; ++j)
               {
                 const size_type loc_col = global_rows.local_row(j);
-                const double col_val = matrix_ptr[loc_col];
+                const LocalType col_val = matrix_ptr[loc_col];
                 dealiiSparseMatrix::add_value(col_val, row,
                                               global_rows.global_row(j),
                                               matrix_values);
@@ -1765,16 +1766,16 @@ namespace internals
                                     loc_row, local_matrix);
             for (size_type j=column_start; j<i; ++j)
               {
-                double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
-                                                       loc_row, local_matrix);
+                LocalType col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
+                                                          loc_row, local_matrix);
                 dealiiSparseMatrix::add_value (col_val, row,
                                                global_rows.global_row(j),
                                                matrix_values);
               }
             for (size_type j=i+1; j<column_end; ++j)
               {
-                double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
-                                                       loc_row, local_matrix);
+                LocalType col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
+                                                          loc_row, local_matrix);
                 dealiiSparseMatrix::add_value (col_val, row,
                                                global_rows.global_row(j),
                                                matrix_values);
@@ -1786,12 +1787,12 @@ namespace internals
       {
         ++matrix_values; // jump over diagonal element
         AssertIndexRange (loc_row, local_matrix.m());
-        const double *matrix_ptr = &local_matrix(loc_row, 0);
+        const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
 
         for (size_type j=column_start; j<column_end; ++j)
           {
             const size_type loc_col = global_rows.local_row(j);
-            const double col_val = matrix_ptr[loc_col];
+            const LocalType col_val = matrix_ptr[loc_col];
             if (row==global_rows.global_row(j))
               sparse_matrix->begin(row)->value() += col_val;
             else
@@ -1805,8 +1806,8 @@ namespace internals
         ++matrix_values; // jump over diagonal element
         for (size_type j=column_start; j<column_end; ++j)
           {
-            double col_val = resolve_matrix_entry (global_rows, global_rows, i,
-                                                   j, loc_row, local_matrix);
+            LocalType col_val = resolve_matrix_entry (global_rows, global_rows, i,
+                                                      j, loc_row, local_matrix);
             if (row==global_rows.global_row(j))
               sparse_matrix->begin(row)->value() += col_val;
             else
@@ -1917,7 +1918,7 @@ add_this_index:
   inline void
   set_matrix_diagonals (const internals::GlobalRowsFromLocal &global_rows,
                         const std::vector<size_type>         &local_dof_indices,
-                        const FullMatrix<double>             &local_matrix,
+                        const FullMatrix<typename MatrixType::value_type>  &local_matrix,
                         const ConstraintMatrix               &constraints,
                         MatrixType                           &global_matrix,
                         VectorType                           &global_vector,
@@ -1925,9 +1926,9 @@ add_this_index:
   {
     if (global_rows.n_constraints() > 0)
       {
-        double average_diagonal = 0;
+        typename MatrixType::value_type average_diagonal = typename MatrixType::value_type();
         for (size_type i=0; i<local_matrix.m(); ++i)
-          average_diagonal += std::fabs (local_matrix(i,i));
+          average_diagonal += std::abs (local_matrix(i,i));
         average_diagonal /= static_cast<double>(local_matrix.m());
 
         for (size_type i=0; i<global_rows.n_constraints(); i++)
@@ -1935,8 +1936,8 @@ add_this_index:
             const size_type local_row = global_rows.constraint_origin(i);
             const size_type global_row = local_dof_indices[local_row];
             const typename MatrixType::value_type new_diagonal
-              = (std::fabs(local_matrix(local_row,local_row)) != 0 ?
-                 std::fabs(local_matrix(local_row,local_row)) : average_diagonal);
+              = (std::abs(local_matrix(local_row,local_row)) != 0 ?
+                 std::abs(local_matrix(local_row,local_row)) : average_diagonal);
             global_matrix.add(global_row, global_row, new_diagonal);
 
             // if the use_inhomogeneities_for_rhs flag is set to true, the
@@ -2121,18 +2122,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>
 inline
-double
+LocalType
 ConstraintMatrix::
 resolve_vector_entry (const size_type                       i,
                       const internals::GlobalRowsFromLocal &global_rows,
-                      const Vector<double>                 &local_vector,
+                      const Vector<LocalType>              &local_vector,
                       const std::vector<size_type>         &local_dof_indices,
-                      const FullMatrix<double>             &local_matrix) const
+                      const FullMatrix<LocalType>          &local_matrix) const
 {
   const size_type loc_row = global_rows.local_row(i);
   const size_type n_inhomogeneous_rows = global_rows.n_inhomogeneities();
-  double val = 0;
+  LocalType 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.
@@ -2150,7 +2152,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);
-      double add_this = local_vector (loc_row_q);
+      LocalType add_this = local_vector (loc_row_q);
       for (size_type k=0; k<n_inhomogeneous_rows; ++k)
         add_this -= (lines[lines_cache[calculate_line_index
                                        (local_dof_indices
@@ -2168,8 +2170,8 @@ resolve_vector_entry (const size_type                       i,
 template <typename MatrixType, typename VectorType>
 void
 ConstraintMatrix::distribute_local_to_global (
-  const FullMatrix<double>        &local_matrix,
-  const Vector<double>            &local_vector,
+  const FullMatrix<typename MatrixType::value_type>     &local_matrix,
+  const Vector<typename VectorType::value_type>         &local_vector,
   const std::vector<size_type>    &local_dof_indices,
   MatrixType                      &global_matrix,
   VectorType                      &global_vector,
@@ -2255,12 +2257,12 @@ ConstraintMatrix::distribute_local_to_global (
       // hand side.
       if (use_vectors == true)
         {
-          const double val = resolve_vector_entry (i, global_rows,
+          const number val = resolve_vector_entry (i, global_rows,
                                                    local_vector,
                                                    local_dof_indices,
                                                    local_matrix);
 
-          if (val != 0)
+          if (val != number ())
             global_vector(row) += static_cast<typename VectorType::value_type>(val);
         }
     }
@@ -2275,12 +2277,12 @@ ConstraintMatrix::distribute_local_to_global (
 template <typename MatrixType>
 void
 ConstraintMatrix::distribute_local_to_global (
-  const FullMatrix<double>     &local_matrix,
+  const FullMatrix<typename MatrixType::value_type>  &local_matrix,
   const std::vector<size_type> &row_indices,
   const std::vector<size_type> &col_indices,
   MatrixType                   &global_matrix) const
 {
-  typedef double number;
+  typedef typename MatrixType::value_type number;
 
   AssertDimension (local_matrix.m(), row_indices.size());
   AssertDimension (local_matrix.n(), col_indices.size());
@@ -2331,8 +2333,8 @@ ConstraintMatrix::distribute_local_to_global (
 template <typename MatrixType, typename VectorType>
 void
 ConstraintMatrix::
-distribute_local_to_global (const FullMatrix<double>     &local_matrix,
-                            const Vector<double>         &local_vector,
+distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>  &local_matrix,
+                            const Vector<typename VectorType::value_type>      &local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType                   &global_matrix,
                             VectorType                   &global_vector,
@@ -2429,12 +2431,12 @@ distribute_local_to_global (const FullMatrix<double>     &local_matrix,
 
           if (use_vectors == true)
             {
-              const double val = resolve_vector_entry (i, global_rows,
+              const number val = resolve_vector_entry (i, global_rows,
                                                        local_vector,
                                                        local_dof_indices,
                                                        local_matrix);
 
-              if (val != 0)
+              if (val != number ())
                 global_vector(global_indices[i]) +=
                   static_cast<typename VectorType::value_type>(val);
             }

diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc
index 479d9f18947c94a3e914944d928ec0328ac64c2c..eb8d83675b60def4ff74aac5c6830bf65f27c6f3 100644 (file)
--- a/source/lac/constraint_matrix.cc
+++ b/source/lac/constraint_matrix.cc
@@ -1247,8 +1247,8 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
 
 #define MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \
   template void ConstraintMatrix:: \
-  distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<double>        &, \
-                                                      const Vector<double>            &, \
+  distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<MatrixType::value_type>        &, \
+                                                      const Vector<VectorType::value_type>            &, \
                                                       const std::vector<ConstraintMatrix::size_type> &, \
                                                       MatrixType                      &, \
                                                       VectorType                      &, \
@@ -1256,17 +1256,17 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
                                                       internal::bool2type<false>) const
 #define MATRIX_FUNCTIONS(MatrixType) \
   template void ConstraintMatrix:: \
-  distribute_local_to_global<MatrixType,Vector<double> > (const FullMatrix<double>        &, \
-                                                          const Vector<double>            &, \
-                                                          const std::vector<ConstraintMatrix::size_type> &, \
-                                                          MatrixType                      &, \
-                                                          Vector<double>                  &, \
-                                                          bool                             , \
-                                                          internal::bool2type<false>) const
+  distribute_local_to_global<MatrixType,Vector<MatrixType::value_type> > (const FullMatrix<MatrixType::value_type>        &, \
+      const Vector<MatrixType::value_type>            &, \
+      const std::vector<ConstraintMatrix::size_type> &, \
+      MatrixType                      &, \
+      Vector<MatrixType::value_type>                  &, \
+      bool                             , \
+      internal::bool2type<false>) const
 #define BLOCK_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType)   \
   template void ConstraintMatrix:: \
-  distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<double>        &, \
-                                                      const Vector<double>            &, \
+  distribute_local_to_global<MatrixType,VectorType > (const FullMatrix<MatrixType::value_type>        &, \
+                                                      const Vector<VectorType::value_type>            &, \
                                                       const std::vector<ConstraintMatrix::size_type> &, \
                                                       MatrixType                      &, \
                                                       VectorType                      &, \
@@ -1274,32 +1274,29 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
                                                       internal::bool2type<true>) const
 #define BLOCK_MATRIX_FUNCTIONS(MatrixType)      \
   template void ConstraintMatrix:: \
-  distribute_local_to_global<MatrixType,Vector<double> > (const FullMatrix<double>        &, \
-                                                          const Vector<double>            &, \
-                                                          const std::vector<ConstraintMatrix::size_type> &, \
-                                                          MatrixType                      &, \
-                                                          Vector<double>                  &, \
-                                                          bool                             , \
-                                                          internal::bool2type<true>) const
+  distribute_local_to_global<MatrixType,Vector<MatrixType::value_type> > (const FullMatrix<MatrixType::value_type>        &, \
+      const Vector<MatrixType::value_type>            &, \
+      const std::vector<ConstraintMatrix::size_type> &, \
+      MatrixType                      &, \
+      Vector<MatrixType::value_type>                  &, \
+      bool                             , \
+      internal::bool2type<true>) const
 
 MATRIX_FUNCTIONS(SparseMatrix<double>);
 MATRIX_FUNCTIONS(SparseMatrix<float>);
 MATRIX_FUNCTIONS(FullMatrix<double>);
 MATRIX_FUNCTIONS(FullMatrix<float>);
-MATRIX_VECTOR_FUNCTIONS(SparseMatrix<float>, Vector<float>);
+MATRIX_FUNCTIONS(FullMatrix<std::complex<double> >);
 
 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>);
-BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<float>,  BlockVector<double>);
 
 MATRIX_FUNCTIONS(SparseMatrixEZ<double>);
 MATRIX_FUNCTIONS(SparseMatrixEZ<float>);
 MATRIX_FUNCTIONS(ChunkSparseMatrix<double>);
 MATRIX_FUNCTIONS(ChunkSparseMatrix<float>);
-MATRIX_VECTOR_FUNCTIONS(SparseMatrixEZ<float>,  Vector<float>);
-MATRIX_VECTOR_FUNCTIONS(ChunkSparseMatrix<float>, Vector<float>);
 
 // BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrixEZ<double>);
 // BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ<float>,  Vector<float>);
@@ -1365,7 +1362,7 @@ BLOCK_SPARSITY_FUNCTIONS(TrilinosWrappers::BlockSparsityPattern);
 
 #define ONLY_MATRIX_FUNCTIONS(MatrixType) \
   template void ConstraintMatrix::distribute_local_to_global<MatrixType > (\
-      const FullMatrix<double>        &, \
+      const FullMatrix<MatrixType::value_type>        &, \
       const std::vector<ConstraintMatrix::size_type> &, \
       const std::vector<ConstraintMatrix::size_type> &, \
       MatrixType                      &) const