From 2dafc0aa45767f6f13dd032c5b1a50439b589741 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 5 Dec 2013 11:35:19 +0000 Subject: [PATCH] Cleanup append_index call. Do not make the scratch data a member of ConstraintMatrix but make it a static variable that lives for the whole program. This simplifies the calls singificantly. git-svn-id: https://svn.dealii.org/trunk@31898 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/lac/constraint_matrix.h | 76 +---- .../deal.II/lac/constraint_matrix.templates.h | 289 +++++++++++------- deal.II/source/lac/constraint_matrix.cc | 21 ++ 3 files changed, 206 insertions(+), 180 deletions(-) diff --git a/deal.II/include/deal.II/lac/constraint_matrix.h b/deal.II/include/deal.II/lac/constraint_matrix.h index e6cecd0ceb..5377861801 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.h @@ -1290,76 +1290,6 @@ private: */ bool sorted; - /** - * Scratch data that is used during calls to distribute_local_to_global and - * add_entries_local_to_global. In order to avoid frequent memory - * allocation, we keep the data alive from one call to the next. - */ - struct ScratchData - { - /** - * Constructor, does nothing. - */ - ScratchData () : - in_use (false) - {} - - /** - * Copy constructor, does nothing - */ - ScratchData (const ScratchData &) : - in_use (false) - {} - - /** - * Stores whether the data is currently in use. - */ - bool in_use; - - /** - * Temporary array for column indices - */ - std::vector columns; - - /** - * Temporary array for column values - */ - std::vector values; - - /** - * Temporary array for block start indices - */ - std::vector block_starts; - - /** - * Temporary array for vector indices - */ - std::vector vector_indices; - - /** - * Data array for reorder row/column indices. Use a shared ptr to - * global_rows to avoid defining in the .h file - */ - std_cxx1x::shared_ptr global_rows; - - /** - * Data array for reorder row/column indices. Use a shared ptr to - * global_rows to avoid defining in the .h file - */ - std_cxx1x::shared_ptr global_columns; - }; - - /** - * Here comes the actual data structure for the scratch data. It is made - * mutable since it is modified in a const function. Since only one thread - * can access it at a time, no conflicting access can occur. For this to be - * valid, we need to make sure that no call within - * distribute_local_to_global is made that by itself can spawn - * tasks. Otherwise, we might end up in a situation where several threads - * fight for the data. - */ - mutable Threads::ThreadLocalStorage scratch_data; - /** * Internal function to calculate the index of line @p line in the vector * lines_cache using local_lines. @@ -1472,8 +1402,7 @@ ConstraintMatrix::ConstraintMatrix (const IndexSet &local_constraints) : lines (), local_lines (local_constraints), - sorted (false), - scratch_data (ScratchData()) + sorted (false) { // make sure the IndexSet is compressed. Otherwise this can lead to crashes // that are hard to find (only happen in release mode). @@ -1490,8 +1419,7 @@ ConstraintMatrix::ConstraintMatrix (const ConstraintMatrix &constraint_matrix) lines (constraint_matrix.lines), lines_cache (constraint_matrix.lines_cache), local_lines (constraint_matrix.local_lines), - sorted (constraint_matrix.sorted), - scratch_data (ScratchData()) + sorted (constraint_matrix.sorted) {} diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h index d2e4dce533..a047165fde 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h @@ -1309,15 +1309,14 @@ namespace internals { AssertDimension(data.size(), individual_size.size()*row_length); // no space left in this row, need to double row_length and - // rearrange the data items + // rearrange the data items. Move all items to the right except the + // first one, starting at the back. Since individual_size contains + // at least one element when we get here, subtracting 1 works fine. data.resize(2*data.size()); - for (size_type i=individual_size.size(); i>0; ) - { - --i; - std::memmove(&data[i*row_length*2], &data[i*row_length], - individual_size[i]* - sizeof(std::pair)); - } + for (size_type i=individual_size.size()-1; i>0; --i) + std::memmove(&data[i*row_length*2], &data[i*row_length], + individual_size[i]* + sizeof(std::pair)); row_length *= 2; } data[index*row_length+my_length] = pair; @@ -1633,6 +1632,142 @@ namespace internals } } + + + /** + * Scratch data that is used during calls to distribute_local_to_global and + * add_entries_local_to_global. In order to avoid frequent memory + * allocation, we keep the data alive from one call to the next in a static + * variable. Since we want to allow for different number types in matrices, + * this is a template. + * + * Since each thread gets its private version of scratch data out of the + * ThreadLocalStorage, no conflicting access can occur. For this to be + * valid, we need to make sure that no call within + * distribute_local_to_global is made that by itself can spawn + * tasks. Otherwise, we might end up in a situation where several threads + * fight for the data. + * + * Access to the scratch data is only through the accessor class which + * handles the access as well as marking the data as used. + */ + template + class ConstraintMatrixData + { + public: + struct ScratchData + { + /** + * Constructor, does nothing. + */ + ScratchData () + : + in_use (false) + {} + + /** + * Copy constructor, does nothing + */ + ScratchData (const ScratchData &) + : + in_use (false) + {} + + /** + * Stores whether the data is currently in use. + */ + bool in_use; + + /** + * Temporary array for column indices + */ + std::vector columns; + + /** + * Temporary array for column values + */ + std::vector values; + + /** + * Temporary array for block start indices + */ + std::vector block_starts; + + /** + * Temporary array for vector indices + */ + std::vector vector_indices; + + /** + * Data array for reorder row/column indices. Use a shared ptr to + * global_rows to avoid defining in the .h file + */ + GlobalRowsFromLocal global_rows; + + /** + * Data array for reorder row/column indices. Use a shared ptr to + * global_rows to avoid defining in the .h file + */ + GlobalRowsFromLocal global_columns; + }; + + /** + * Accessor class to guard access to scratch_data + */ + class ScratchDataAccessor + { + public: + /** + * Constructor. Grabs a scratch data object on the current thread and + * mark it as used + */ + ScratchDataAccessor() + : + my_scratch_data(&ConstraintMatrixData::scratch_data.get()) + { + Assert(my_scratch_data->in_use == false, + ExcMessage("Access to thread-local scratch data tried, but it is already " + "in use")); + my_scratch_data->in_use = true; + } + + /** + * Destructor. Mark scratch data as available again. + */ + ~ScratchDataAccessor() + { + my_scratch_data->in_use = false; + } + + /** + * Dereferencing operator. + */ + ScratchData & operator* () + { + return *my_scratch_data; + } + + /** + * Dereferencing operator. + */ + ScratchData * operator-> () + { + return my_scratch_data; + } + + private: + ScratchData* my_scratch_data; + }; + + private: + /** + * The actual data object that contains a scratch data for each thread. + */ + static Threads::ThreadLocalStorage scratch_data; + }; + + + // function for block matrices: Find out where in the list of local dofs // (sorted according to global ids) the individual blocks start. Transform // the global indices to block-local indices in order to be able to use @@ -2373,17 +2508,10 @@ ConstraintMatrix::distribute_local_to_global ( const size_type n_local_dofs = local_dof_indices.size(); - ScratchData &my_scratch_data = scratch_data.get(); - Assert(my_scratch_data.in_use == false, - ExcMessage("Access to thread-local scratch data tried, but it is already " - "in use")); - // TODO: might want to have a scoped variable for in_use here and in the - // methods below - my_scratch_data.in_use = true; - - if (my_scratch_data.global_rows.get() == 0) - my_scratch_data.global_rows.reset(new internals::GlobalRowsFromLocal()); - internals::GlobalRowsFromLocal &global_rows = *my_scratch_data.global_rows; + typename internals::ConstraintMatrixData::ScratchDataAccessor + scratch_data; + + internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows; global_rows.reinit(n_local_dofs); make_sorted_row_list (local_dof_indices, global_rows); @@ -2395,18 +2523,14 @@ ConstraintMatrix::distribute_local_to_global ( // an array in any case since we cannot know about the actual data type in // the ConstraintMatrix class (unless we do cast). This involves a little // bit of logic to determine the type of the matrix value. - std::vector & cols = my_scratch_data.columns; - std::vector & vals = my_scratch_data.values; - std::vector values_non_double; + std::vector & cols = scratch_data->columns; + std::vector & vals = scratch_data->values; SparseMatrix *sparse_matrix = dynamic_cast *>(&global_matrix); if (use_dealii_matrix == false) { cols.resize (n_actual_dofs); - if (types_are_equal::value == false) - values_non_double.resize(n_actual_dofs); - else - vals.resize (n_actual_dofs); + vals.resize (n_actual_dofs); } else Assert (sparse_matrix != 0, ExcInternalError()); @@ -2423,14 +2547,13 @@ ConstraintMatrix::distribute_local_to_global ( size_type *col_ptr = &cols[0]; // cast is uncritical here and only used to avoid compiler // warnings. We never access a non-double array - number *val_ptr = types_are_equal::value ? - reinterpret_cast(&vals[0]) : &values_non_double[0]; + number *val_ptr = &vals[0]; internals::resolve_matrix_row (global_rows, global_rows, i, 0, n_actual_dofs, local_matrix, col_ptr, val_ptr); const size_type n_values = col_ptr - &cols[0]; if (n_values > 0) - global_matrix.add(row, n_values, &cols[0], val_ptr-n_values, false, + global_matrix.add(row, n_values, &cols[0], &vals[0], false, true); } else @@ -2457,7 +2580,6 @@ ConstraintMatrix::distribute_local_to_global ( internals::set_matrix_diagonals (global_rows, local_dof_indices, local_matrix, *this, global_matrix, global_vector, use_inhomogeneities_for_rhs); - my_scratch_data.in_use = false; } @@ -2479,19 +2601,11 @@ 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(); - ScratchData &my_scratch_data = scratch_data.get(); - Assert(my_scratch_data.in_use == false, - ExcMessage("Access to thread-local scratch data tried, but it is already " - "in use")); - my_scratch_data.in_use = true; - - if (my_scratch_data.global_rows.get() == 0) - my_scratch_data.global_rows.reset(new internals::GlobalRowsFromLocal()); - if (my_scratch_data.global_columns.get() == 0) - my_scratch_data.global_columns.reset(new internals::GlobalRowsFromLocal()); - internals::GlobalRowsFromLocal &global_rows = *my_scratch_data.global_rows; + typename internals::ConstraintMatrixData::ScratchDataAccessor + scratch_data; + internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows; global_rows.reinit(n_local_row_dofs); - internals::GlobalRowsFromLocal &global_cols = *my_scratch_data.global_columns; + internals::GlobalRowsFromLocal &global_cols = scratch_data->global_columns; global_cols.reinit(n_local_col_dofs); make_sorted_row_list (row_indices, global_rows); make_sorted_row_list (col_indices, global_cols); @@ -2501,14 +2615,10 @@ ConstraintMatrix::distribute_local_to_global ( // create arrays for the column data (indices and values) that will then be // written into the matrix. Shortcut for deal.II sparse matrix - std::vector & cols = my_scratch_data.columns; - std::vector & vals = my_scratch_data.values; - std::vector values_non_double; + std::vector & cols = scratch_data->columns; + std::vector & vals = scratch_data->values; cols.resize(n_actual_col_dofs); - if (types_are_equal::value == true) - vals.resize(n_actual_col_dofs); - else - values_non_double.resize(n_actual_col_dofs); + vals.resize(n_actual_col_dofs); // now do the actual job. for (size_type i=0; i::value ? - reinterpret_cast(&vals[0]) : &values_non_double[0]; + number *val_ptr = &vals[0]; internals::resolve_matrix_row (global_rows, global_cols, i, 0, n_actual_col_dofs, local_matrix, col_ptr, val_ptr); const size_type n_values = col_ptr - &cols[0]; if (n_values > 0) - global_matrix.add(row, n_values, &cols[0], val_ptr-n_values, - false, true); + global_matrix.add(row, n_values, &cols[0], &vals[0], false, true); } - - my_scratch_data.in_use = false; } @@ -2563,22 +2669,17 @@ distribute_local_to_global (const FullMatrix &local_matrix, } Assert (sorted == true, ExcMatrixNotClosed()); - ScratchData &my_scratch_data = scratch_data.get(); - Assert(my_scratch_data.in_use == false, - ExcMessage("Access to thread-local scratch data tried, but it is already " - "in use")); - my_scratch_data.in_use = true; + typename internals::ConstraintMatrixData::ScratchDataAccessor + scratch_data; const size_type n_local_dofs = local_dof_indices.size(); - if (my_scratch_data.global_rows.get() == 0) - my_scratch_data.global_rows.reset(new internals::GlobalRowsFromLocal()); - internals::GlobalRowsFromLocal &global_rows = *my_scratch_data.global_rows; + internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows; global_rows.reinit(n_local_dofs); make_sorted_row_list (local_dof_indices, global_rows); const size_type n_actual_dofs = global_rows.size(); - std::vector &global_indices = my_scratch_data.vector_indices; + std::vector &global_indices = scratch_data->vector_indices; if (use_vectors == true) { global_indices.resize(n_actual_dofs); @@ -2588,20 +2689,16 @@ distribute_local_to_global (const FullMatrix &local_matrix, // additional construct that also takes care of block indices. const size_type num_blocks = global_matrix.n_block_rows(); - std::vector &block_starts = my_scratch_data.block_starts; + std::vector &block_starts = scratch_data->block_starts; block_starts.resize(num_blocks+1); internals::make_block_starts (global_matrix, global_rows, block_starts); - std::vector & cols = my_scratch_data.columns; - std::vector & vals = my_scratch_data.values; - std::vector values_non_double; + std::vector & cols = scratch_data->columns; + std::vector & vals = scratch_data->values; if (use_dealii_matrix == false) { cols.resize (n_actual_dofs); - if (types_are_equal::value == true) - vals.resize(n_actual_dofs); - else - values_non_double.resize(n_actual_dofs); + vals.resize (n_actual_dofs); } // the basic difference to the non-block variant from now onwards is that we @@ -2621,16 +2718,14 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (use_dealii_matrix == false) { size_type *col_ptr = &cols[0]; - number *val_ptr = types_are_equal::value ? - reinterpret_cast(&vals[0]) : &values_non_double[0]; + number *val_ptr = &vals[0]; internals::resolve_matrix_row (global_rows, global_rows, i, start_block, end_block, local_matrix, col_ptr, val_ptr); const size_type n_values = col_ptr - &cols[0]; if (n_values > 0) global_matrix.block(block, block_col).add(row, n_values, - &cols[0], - val_ptr-n_values, + &cols[0], &vals[0], false, true); } else @@ -2661,8 +2756,6 @@ distribute_local_to_global (const FullMatrix &local_matrix, internals::set_matrix_diagonals (global_rows, local_dof_indices, local_matrix, *this, global_matrix, global_vector, use_inhomogeneities_for_rhs); - - my_scratch_data.in_use = false; } @@ -2686,11 +2779,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, AssertDimension (dof_mask.n_cols(), n_local_dofs); } - ScratchData &my_scratch_data = scratch_data.get(); - Assert(my_scratch_data.in_use == false, - ExcMessage("Access to thread-local scratch data tried, but it is already " - "in use")); - my_scratch_data.in_use = true; + internals::ConstraintMatrixData::ScratchDataAccessor scratch_data; // if the dof mask is not active, all we have to do is to add some indices // in a matrix format. To do this, we first create an array of all the @@ -2698,7 +2787,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, // plus some indices that come from constraints. if (dof_mask_is_active == false) { - std::vector & actual_dof_indices = my_scratch_data.columns; + std::vector & actual_dof_indices = scratch_data->columns; actual_dof_indices.resize(n_local_dofs); make_sorted_row_list (local_dof_indices, actual_dof_indices); const size_type n_actual_dofs = actual_dof_indices.size(); @@ -2726,25 +2815,22 @@ add_entries_local_to_global (const std::vector &local_dof_indices, else sparsity_pattern.add (local_dof_indices[i], local_dof_indices[i]); } - - my_scratch_data.in_use = false; - return; + + return; } // complicated case: we need to filter out some indices. then the function // gets similar to the function for distributing matrix entries, see there // for additional comments. - if (my_scratch_data.global_rows.get() == 0) - my_scratch_data.global_rows.reset(new internals::GlobalRowsFromLocal()); - internals::GlobalRowsFromLocal &global_rows = *my_scratch_data.global_rows; + internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows; global_rows.reinit(n_local_dofs); make_sorted_row_list (local_dof_indices, global_rows); const size_type n_actual_dofs = global_rows.size(); // create arrays for the column indices that will then be written into the // sparsity pattern. - std::vector & cols = my_scratch_data.columns; + std::vector & cols = scratch_data->columns; cols.resize(n_actual_dofs); for (size_type i=0; i &local_dof_indices, internals::set_sparsity_diagonals (global_rows, local_dof_indices, dof_mask, keep_constrained_entries, sparsity_pattern); - my_scratch_data.in_use = false; } @@ -2846,11 +2931,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, const size_type n_local_dofs = local_dof_indices.size(); const size_type num_blocks = sparsity_pattern.n_block_rows(); - ScratchData &my_scratch_data = scratch_data.get(); - Assert(my_scratch_data.in_use == false, - ExcMessage("Access to thread-local scratch data tried, but it is already " - "in use")); - my_scratch_data.in_use = true; + internals::ConstraintMatrixData::ScratchDataAccessor scratch_data; bool dof_mask_is_active = false; if (dof_mask.n_rows() == n_local_dofs) @@ -2861,13 +2942,13 @@ add_entries_local_to_global (const std::vector &local_dof_indices, if (dof_mask_is_active == false) { - std::vector & actual_dof_indices = my_scratch_data.columns; + std::vector & actual_dof_indices = scratch_data->columns; actual_dof_indices.resize(n_local_dofs); make_sorted_row_list (local_dof_indices, actual_dof_indices); const size_type n_actual_dofs = actual_dof_indices.size(); // additional construct that also takes care of block indices. - std::vector &block_starts = my_scratch_data.block_starts; + std::vector &block_starts = scratch_data->block_starts; block_starts.resize(num_blocks+1); internals::make_block_starts (sparsity_pattern, actual_dof_indices, block_starts); @@ -2908,25 +2989,22 @@ add_entries_local_to_global (const std::vector &local_dof_indices, sparsity_pattern.add (local_dof_indices[i], local_dof_indices[i]); } - my_scratch_data.in_use = false; return; } // difficult case with dof_mask, similar to the distribute_local_to_global // function for block matrices - if (my_scratch_data.global_rows.get() == 0) - my_scratch_data.global_rows.reset(new internals::GlobalRowsFromLocal()); - internals::GlobalRowsFromLocal &global_rows = *my_scratch_data.global_rows; + internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows; global_rows.reinit(n_local_dofs); make_sorted_row_list (local_dof_indices, global_rows); const size_type n_actual_dofs = global_rows.size(); // additional construct that also takes care of block indices. - std::vector & block_starts = my_scratch_data.block_starts; + std::vector & block_starts = scratch_data->block_starts; block_starts.resize(num_blocks+1); internals::make_block_starts(sparsity_pattern, global_rows, block_starts); - std::vector &cols = my_scratch_data.columns; + std::vector &cols = scratch_data->columns; cols.resize(n_actual_dofs); // the basic difference to the non-block variant from now onwards is that we @@ -2956,7 +3034,6 @@ add_entries_local_to_global (const std::vector &local_dof_indices, internals::set_sparsity_diagonals (global_rows, local_dof_indices, dof_mask, keep_constrained_entries, sparsity_pattern); - my_scratch_data.in_use = false; } diff --git a/deal.II/source/lac/constraint_matrix.cc b/deal.II/source/lac/constraint_matrix.cc index 5ace029bd2..8c8c86cb3c 100644 --- a/deal.II/source/lac/constraint_matrix.cc +++ b/deal.II/source/lac/constraint_matrix.cc @@ -1900,4 +1900,25 @@ ONLY_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix); #include "constraint_matrix.inst" +// allocate scratch data. Cannot use the generic template instantiation +// because we need to provide an initializer object of type +// internals::ConstraintMatrixData that can be passed to the +// constructor of scratch_data (it won't allow one to be constructed in place). +namespace internals +{ +#define SCRATCH_INITIALIZER(Number,Name) \ + ConstraintMatrixData::ScratchData scratch_data_initializer_##Name; \ + template<> Threads::ThreadLocalStorage::ScratchData> \ + ConstraintMatrixData::scratch_data(scratch_data_initializer_##Name) + + SCRATCH_INITIALIZER(double,double); + SCRATCH_INITIALIZER(float,float); + SCRATCH_INITIALIZER(long double,ldouble); + SCRATCH_INITIALIZER(std::complex,cdouble); + SCRATCH_INITIALIZER(std::complex,cfloat); + SCRATCH_INITIALIZER(std::complex,cldouble); +#undef SCRATCH_INITIALIZER +} + + DEAL_II_NAMESPACE_CLOSE -- 2.39.5