From b28c73331e4ae3fc381ee2eb8aad3f85a94ee68d Mon Sep 17 00:00:00 2001
From: Karl Ljungkvist <k.ljungkvist@gmail.com>
Date: Mon, 21 Nov 2016 14:53:25 +0100
Subject: [PATCH] extract vectorized data structure from common code

---
 .../deal.II/multigrid/mg_transfer_internal.h  |  2 +-
 source/multigrid/mg_transfer_internal.cc      | 26 +++++++----------
 source/multigrid/mg_transfer_internal.inst.in |  2 +-
 source/multigrid/mg_transfer_matrix_free.cc   | 28 ++++++++++++++++++-
 4 files changed, 39 insertions(+), 19 deletions(-)

diff --git a/include/deal.II/multigrid/mg_transfer_internal.h b/include/deal.II/multigrid/mg_transfer_internal.h
index dec3ea3b11..80f9d291bc 100644
--- a/include/deal.II/multigrid/mg_transfer_internal.h
+++ b/include/deal.II/multigrid/mg_transfer_internal.h
@@ -106,7 +106,7 @@ namespace internal
                         std::vector<std::vector<std::pair<unsigned int,unsigned int> > >    &parent_child_connect,
                         std::vector<unsigned int>                                           &n_owned_level_cells,
                         std::vector<std::vector<std::vector<unsigned short> > >             &dirichlet_indices,
-                        std::vector<AlignedVector<VectorizedArray<Number> > >               &weights_on_refined,
+                        std::vector<std::vector<Number> >                                   &weights_on_refined,
                         std::vector<std::vector<std::pair<unsigned int, unsigned int> > >   &copy_indices_global_mine,
                         MGLevelObject<LinearAlgebra::distributed::Vector<Number> >          &ghosted_level_vector);
 
diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc
index 253292a582..c560bb5e59 100644
--- a/source/multigrid/mg_transfer_internal.cc
+++ b/source/multigrid/mg_transfer_internal.cc
@@ -456,7 +456,7 @@ namespace internal
                         std::vector<std::vector<std::pair<unsigned int,unsigned int> > >    &parent_child_connect,
                         std::vector<unsigned int>                                           &n_owned_level_cells,
                         std::vector<std::vector<std::vector<unsigned short> > >             &dirichlet_indices,
-                        std::vector<AlignedVector<VectorizedArray<Number> > >               &weights_on_refined,
+                        std::vector<std::vector<Number> >                                   &weights_on_refined,
                         std::vector<std::vector<std::pair<unsigned int, unsigned int> > >   &copy_indices_global_mine,
                         MGLevelObject<LinearAlgebra::distributed::Vector<Number> >          &ghosted_level_vector)
     {
@@ -709,7 +709,6 @@ namespace internal
           ghosted_level_vector[level].compress(VectorOperation::add);
           ghosted_level_vector[level].update_ghost_values();
 
-          const unsigned int vec_size = VectorizedArray<Number>::n_array_elements;
           std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
           degree_to_3[0] = 0;
           for (unsigned int i=1; i<n_child_dofs_1d-1; ++i)
@@ -718,21 +717,16 @@ namespace internal
 
           // we only store 3^dim weights because all dofs on a line have the
           // same valence, and all dofs on a quad have the same valence.
-          weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*Utilities::fixed_power<dim>(3));
+          weights_on_refined[level-1].resize(n_owned_level_cells[level-1]*Utilities::fixed_power<dim>(3));
           for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
-            {
-              const unsigned int comp = c/vec_size;
-              const unsigned int v = c%vec_size;
-
-              for (unsigned int k=0, m=0; k<(dim>2 ? n_child_dofs_1d : 1); ++k)
-                for (unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j)
-                  {
-                    unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
-                    for (unsigned int i=0; i<n_child_dofs_1d; ++i, ++m)
-                      weights_on_refined[level-1][comp*Utilities::fixed_power<dim>(3)+shift+degree_to_3[i]][v] = Number(1.)/
-                          ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+m]);
-                  }
-            }
+            for (unsigned int k=0, m=0; k<(dim>2 ? n_child_dofs_1d : 1); ++k)
+              for (unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j)
+                {
+                  unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
+                  for (unsigned int i=0; i<n_child_dofs_1d; ++i, ++m)
+                    weights_on_refined[level-1][c*Utilities::fixed_power<dim>(3)+shift+degree_to_3[i]] = Number(1.)/
+                        ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+m]);
+                }
         }
 
     }
diff --git a/source/multigrid/mg_transfer_internal.inst.in b/source/multigrid/mg_transfer_internal.inst.in
index 96c231acbd..d1b848dd9b 100644
--- a/source/multigrid/mg_transfer_internal.inst.in
+++ b/source/multigrid/mg_transfer_internal.inst.in
@@ -81,7 +81,7 @@ for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS)
                                            std::vector<std::vector<std::pair<unsigned int,unsigned int> > >&,
                                            std::vector<unsigned int>&,
                                            std::vector<std::vector<std::vector<unsigned short> > >&,
-                                           std::vector<AlignedVector<VectorizedArray<S> > >&,
+                                           std::vector<std::vector<S> >&,
                                            std::vector<std::vector<std::pair<unsigned int, unsigned int> > >&,
                                            MGLevelObject<LinearAlgebra::distributed::Vector<S> >&);
     \}
diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc
index e23881f1b3..fefe0ec73b 100644
--- a/source/multigrid/mg_transfer_matrix_free.cc
+++ b/source/multigrid/mg_transfer_matrix_free.cc
@@ -97,7 +97,10 @@ void MGTransferMatrixFree<dim,Number>::build
 {
   this->fill_and_communicate_copy_indices(mg_dof);
 
+  std::vector<std::vector<Number> > weights_unvectorized;
+
   internal::MGTransfer::ElementInfo<Number> elem_info;
+
   internal::MGTransfer::setup_transfer<dim,Number>(mg_dof,
                                                    (const MGConstrainedDoFs *)this->mg_constrained_dofs,
                                                    elem_info,
@@ -105,7 +108,7 @@ void MGTransferMatrixFree<dim,Number>::build
                                                    parent_child_connect,
                                                    n_owned_level_cells,
                                                    dirichlet_indices,
-                                                   weights_on_refined,
+                                                   weights_unvectorized,
                                                    this->copy_indices_global_mine,
                                                    this->ghosted_level_vector);
   // unpack element info data
@@ -115,6 +118,29 @@ void MGTransferMatrixFree<dim,Number>::build
   n_child_cell_dofs        = elem_info.n_child_cell_dofs;
   shape_info               = elem_info.shape_info;
 
+
+  // reshuffle into aligned vector of vectorized arrays
+  const unsigned int vec_size = VectorizedArray<Number>::n_array_elements;
+  const unsigned int n_levels = mg_dof.get_triangulation().n_global_levels();
+
+  const unsigned int n_weights_per_cell = Utilities::fixed_power<dim>(3);
+  weights_on_refined.resize(n_levels-1);
+  for (unsigned int level = 1; level<n_levels; ++level)
+    {
+      weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*n_weights_per_cell);
+
+      for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
+        {
+          const unsigned int comp = c/vec_size;
+          const unsigned int v = c%vec_size;
+          for (unsigned int i = 0; i<n_weights_per_cell; ++i)
+            {
+
+              weights_on_refined[level-1][comp*n_weights_per_cell+i][v] = weights_unvectorized[level-1][c*n_weights_per_cell+i];
+            }
+        }
+    }
+
   evaluation_data.resize(3*n_child_cell_dofs);
 }
 
-- 
2.39.5