From 8b64a4ae18d27c8fff5aab7547adc36a1f07c489 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 6 Apr 2009 15:47:02 +0000 Subject: [PATCH] Move some collective add/set functions into the cc file to not make the header file too big. git-svn-id: https://svn.dealii.org/trunk@18558 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/compressed_simple_sparsity_pattern.h | 129 +--------- deal.II/lac/include/lac/sparse_matrix.h | 243 +++--------------- .../lac/include/lac/sparse_matrix.templates.h | 195 ++++++++++++++ .../compressed_simple_sparsity_pattern.cc | 160 ++++++++++++ deal.II/lac/source/sparse_matrix.inst.in | 13 + 5 files changed, 403 insertions(+), 337 deletions(-) diff --git a/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h index 0d25f1dddc..d588f72f7b 100644 --- a/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h +++ b/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h @@ -197,7 +197,7 @@ class CompressedSimpleSparsityPattern : public Subscriptor void add_entries (const unsigned int row, ForwardIterator begin, ForwardIterator end, - const bool indices_are_sorted = false); + const bool indices_are_unique_and_sorted = false); /** * Check if a value at a certain @@ -415,133 +415,6 @@ CompressedSimpleSparsityPattern::Line::add (const unsigned int j) -template -inline -void -CompressedSimpleSparsityPattern::Line::add_entries (ForwardIterator begin, - ForwardIterator end, - const bool indices_are_sorted) -{ - const int n_elements = end - begin; - if (n_elements <= 0) - return; - - const unsigned int n_cols = static_cast(n_elements); - const unsigned int stop_size = entries.size() + n_cols; - - if (indices_are_sorted == true && n_elements > 3) - { - if (entries.size() == 0 || entries.back() < *begin) - { - entries.insert(entries.end(), begin, end); - return; - } - - // resize vector by just inserting the - // list - const unsigned int col = *begin; - std::vector::iterator it = - std::lower_bound(entries.begin(), entries.end(), col); - const unsigned int pos1 = it - entries.begin(); - entries.insert (it, begin, end); - it = entries.begin() + pos1; - - // now merge the two lists. - ForwardIterator my_it = begin; - std::vector::iterator it2 = it + n_cols; - - // as long as there are indices both in - // the end of the entries list and in the - // input list - while (my_it != end && it2 != entries.end()) - { - if (*my_it < *it2) - *it++ = *my_it++; - else if (*my_it == *it2) - { - *it++ = *it2++; - ++my_it; - } - else - *it++ = *it2++; - } - // in case there are indices left in the - // input list - while (my_it != end) - *it++ = *my_it++; - - // in case there are indices left in the - // end of entries - while (it2 != entries.end()) - *it++ = *it2++; - - // resize - const unsigned int new_size = it - entries.begin(); - Assert (new_size <= stop_size, ExcInternalError()); - entries.resize (new_size); - return; - } - - ForwardIterator my_it = begin; - - // If necessary, increase the size of the array. In order to avoid - // allocating just a few entries every time, use five elements at a time - // at least. - if (stop_size > entries.capacity()) - entries.reserve (stop_size); - - unsigned int col = *my_it; - std::vector::iterator it, it2; - // insert the first element as for one entry only first check the last - // element (or if line is still empty) - if ( (entries.size()==0) || ( entries.back() < col) ) { - entries.push_back(col); - it = entries.end()-1; - } - else { - // do a binary search to find the place where to insert: - it2 = std::lower_bound(entries.begin(), entries.end(), col); - - // If this entry is a duplicate, continue immediately Insert at the - // right place in the vector. Vector grows automatically to fit - // elements. Always doubles its size. - if (*it2 != col) - it = entries.insert(it2, col); - else - it = it2; - } - - ++my_it; - // Now try to be smart and insert with bias in the direction we are - // walking. This has the advantage that for sorted lists, we always search - // in the right direction, what should decrease the work needed in here. - for ( ; my_it != end; ++my_it) - { - col = *my_it; - // need a special insertion command when we're at the end of the list - if (col > entries.back()) { - entries.push_back(col); - it = entries.end()-1; - } - // search to the right (preferred search direction) - else if (col > *it) { - it2 = std::lower_bound(it++, entries.end(), col); - if (*it2 != col) - it = entries.insert(it2, col); - } - // search to the left - else if (col < *it) { - it2 = std::lower_bound(entries.begin(), it, col); - if (*it2 != col) - it = entries.insert(it2, col); - } - // if we're neither larger nor smaller, then this was a duplicate and - // we can just continue. - } -} - - - inline unsigned int CompressedSimpleSparsityPattern::n_rows () const diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index fdf0d34e96..76ddb7edae 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -2019,11 +2019,7 @@ unsigned int SparseMatrix::n () const // Inline the set() and add() // functions, since they will be - // called frequently, and the - // compiler can optimize away - // some unnecessary loops when - // the sizes are given at - // compile time. + // called frequently. template inline void @@ -2031,12 +2027,25 @@ SparseMatrix::set (const unsigned int i, const unsigned int j, const number value) { - Assert (numbers::is_finite(value), ExcMessage("The given value is not finite but either " "infinite or Not A Number (NaN)")); - set (i, 1, &j, &value, false); + const unsigned int index = cols->operator()(i, j); + + // it is allowed to set elements of + // the matrix that are not part of + // the sparsity pattern, if the + // value to which we set it is zero + if (index == SparsityPattern::invalid_entry) + { + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(i, j)); + return; + } + + val[index] = value; } @@ -2099,71 +2108,6 @@ SparseMatrix::set (const unsigned int row, -template -template -inline -void -SparseMatrix::set (const unsigned int row, - const unsigned int n_cols, - const unsigned int *col_indices, - const number2 *values, - const bool elide_zero_values) -{ - Assert (cols != 0, ExcNotInitialized()); - - unsigned int n_columns = 0; - - // Otherwise, extract nonzero values in - // each row and pass on to the other - // function. - global_indices.resize(n_cols); - column_values.resize(n_cols); - - // First, search all the indices to find - // out which values we actually need to - // set. - // TODO: Could probably made more - // efficient. - for (unsigned int j=0; joperator()(row, col_indices[j]); - - // it is allowed to set elements in - // the matrix that are not part of - // the sparsity pattern, if the - // value to which we set it is zero - if (index == SparsityPattern::invalid_entry) - { - Assert ((index != SparsityPattern::invalid_entry) || - (value == 0.), - ExcInvalidIndex(row,col_indices[j])); - continue; - } - - global_indices[n_columns] = index; - column_values[n_columns] = value; - n_columns++; - } - - const unsigned int * index_ptr = &global_indices[0]; - const number * value_ptr = &column_values[0]; - - // Finally, go through the index list - // and set the elements one by one. - for (unsigned int j=0; j inline void @@ -2171,12 +2115,28 @@ SparseMatrix::add (const unsigned int i, const unsigned int j, const number value) { - Assert (numbers::is_finite(value), ExcMessage("The given value is not finite but either " "infinite or Not A Number (NaN)")); - add (i, 1, &j, &value, false); + if (value == 0) + return; + + const unsigned int index = cols->operator()(i, j); + + // it is allowed to add elements to + // the matrix that are not part of + // the sparsity pattern, if the + // value to which we set it is zero + if (index == SparsityPattern::invalid_entry) + { + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(i, j)); + return; + } + + val[index] += value; } @@ -2239,141 +2199,6 @@ SparseMatrix::add (const unsigned int row, -template -template -inline -void -SparseMatrix::add (const unsigned int row, - const unsigned int n_cols, - const unsigned int *col_indices, - const number2 *values, - const bool elide_zero_values, - const bool col_indices_are_sorted) -{ - Assert (cols != 0, ExcNotInitialized()); - - // if we have sufficiently many columns - // and sorted indices it is faster to - // just go through the column indices and - // look whether we found one, rather than - // doing a binary search for every column - if (col_indices_are_sorted == true) - if (n_cols * 8 > cols->row_length(row)) - { - // check whether the given indices are - // really sorted -#ifdef DEBUG - for (unsigned int i=1; i col_indices[i-1], - ExcMessage("Indices are not sorted.")); -#endif - if (cols->optimize_diagonal() == true) - { - const unsigned int * this_cols = - &cols->get_column_numbers()[cols->get_rowstart_indices()[row]]; - number * val_ptr = &val[cols->get_rowstart_indices()[row]]; - Assert (this_cols[0] == row, ExcInternalError()); - unsigned int counter = 1; - for (unsigned int i=0; i= this_cols[counter], ExcInternalError()); - - while (this_cols[counter] < col_indices[i]) - ++counter; - - Assert (this_cols[counter] == col_indices[i] || values[i] == 0, - ExcInvalidIndex(row,col_indices[i])); - - if (values[i] != 0) - val_ptr[counter] += values[i]; - } - Assert (counter < cols->row_length(row), ExcInternalError()); - } - else - { - const unsigned int * this_cols = - &cols->get_column_numbers()[cols->get_rowstart_indices()[row]]; - number * val_ptr = &val[cols->get_rowstart_indices()[row]]; - unsigned int counter = 0; - for (unsigned int i=0; i= this_cols[counter], ExcInternalError()); - - while (this_cols[counter] < col_indices[i]) - ++counter; - - Assert (this_cols[counter] == col_indices[i] || values[i] == 0, - ExcInvalidIndex(row,col_indices[i])); - - if (values[i] != 0) - val_ptr[counter] += values[i]; - } - Assert (counter < cols->row_length(row), ExcInternalError()); - } - return; - } - - unsigned int n_columns = 0; - - if (global_indices.size() < n_cols) - { - global_indices.resize(n_cols); - column_values.resize(n_cols); - } - - // First, search all the indices to find - // out which values we actually need to - // add. - // TODO: Could probably be made more - // efficient. - for (unsigned int j=0; joperator()(row, col_indices[j]); - - // it is allowed to add elements to - // the matrix that are not part of - // the sparsity pattern, if the - // value to which we set it is zero - if (index == SparsityPattern::invalid_entry) - { - Assert ((index != SparsityPattern::invalid_entry) || - (value == 0.), - ExcInvalidIndex(row,col_indices[j])); - continue; - } - - global_indices[n_columns] = index; - column_values[n_columns] = value; - n_columns++; - } - - const unsigned int * index_ptr = &global_indices[0]; - const number * value_ptr = &column_values[0]; - - // Finally, go through the index list - // and add the elements one by one. - for (unsigned int j=0; j inline SparseMatrix & diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 0816f00ec3..21639c9c12 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -325,6 +325,201 @@ SparseMatrix::add (const number factor, } + +template +template +void +SparseMatrix::add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number2 *values, + const bool elide_zero_values, + const bool col_indices_are_sorted) +{ + Assert (cols != 0, ExcNotInitialized()); + + // if we have sufficiently many columns + // and sorted indices it is faster to + // just go through the column indices and + // look whether we found one, rather than + // doing a binary search for every column + if (col_indices_are_sorted == true) + if (n_cols * 8 > cols->row_length(row)) + { + // check whether the given indices are + // really sorted +#ifdef DEBUG + for (unsigned int i=1; i col_indices[i-1], + ExcMessage("Indices are not sorted.")); +#endif + if (cols->optimize_diagonal() == true) + { + const unsigned int * this_cols = + &cols->get_column_numbers()[cols->get_rowstart_indices()[row]]; + number * val_ptr = &val[cols->get_rowstart_indices()[row]]; + Assert (this_cols[0] == row, ExcInternalError()); + unsigned int counter = 1; + for (unsigned int i=0; i= this_cols[counter], ExcInternalError()); + + while (this_cols[counter] < col_indices[i]) + ++counter; + + Assert (this_cols[counter] == col_indices[i] || values[i] == 0, + ExcInvalidIndex(row,col_indices[i])); + + if (values[i] != 0) + val_ptr[counter] += values[i]; + } + Assert (counter < cols->row_length(row), ExcInternalError()); + } + else + { + const unsigned int * this_cols = + &cols->get_column_numbers()[cols->get_rowstart_indices()[row]]; + number * val_ptr = &val[cols->get_rowstart_indices()[row]]; + unsigned int counter = 0; + for (unsigned int i=0; i= this_cols[counter], ExcInternalError()); + + while (this_cols[counter] < col_indices[i]) + ++counter; + + Assert (this_cols[counter] == col_indices[i] || values[i] == 0, + ExcInvalidIndex(row,col_indices[i])); + + if (values[i] != 0) + val_ptr[counter] += values[i]; + } + Assert (counter < cols->row_length(row), ExcInternalError()); + } + return; + } + + unsigned int n_columns = 0; + + if (global_indices.size() < n_cols) + { + global_indices.resize(n_cols); + column_values.resize(n_cols); + } + + // First, search all the indices to find + // out which values we actually need to + // add. + for (unsigned int j=0; joperator()(row, col_indices[j]); + + // it is allowed to add elements to + // the matrix that are not part of + // the sparsity pattern, if the + // value to which we set it is zero + if (index == SparsityPattern::invalid_entry) + { + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(row,col_indices[j])); + continue; + } + + global_indices[n_columns] = index; + column_values[n_columns] = value; + n_columns++; + } + + const unsigned int * index_ptr = &global_indices[0]; + const number * value_ptr = &column_values[0]; + + // Finally, go through the index list + // and add the elements one by one. + for (unsigned int j=0; j +template +void +SparseMatrix::set (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number2 *values, + const bool elide_zero_values) +{ + Assert (cols != 0, ExcNotInitialized()); + + unsigned int n_columns = 0; + + // Otherwise, extract nonzero values in + // each row and pass on to the other + // function. + global_indices.resize(n_cols); + column_values.resize(n_cols); + + // First, search all the indices to find + // out which values we actually need to + // set. + for (unsigned int j=0; joperator()(row, col_indices[j]); + + // it is allowed to set elements in + // the matrix that are not part of + // the sparsity pattern, if the + // value to which we set it is zero + if (index == SparsityPattern::invalid_entry) + { + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(row,col_indices[j])); + continue; + } + + global_indices[n_columns] = index; + column_values[n_columns] = value; + n_columns++; + } + + const unsigned int * index_ptr = &global_indices[0]; + const number * value_ptr = &column_values[0]; + + // Finally, go through the index list + // and set the elements one by one. + for (unsigned int j=0; j template void diff --git a/deal.II/lac/source/compressed_simple_sparsity_pattern.cc b/deal.II/lac/source/compressed_simple_sparsity_pattern.cc index 4a3c702434..c3ae666ded 100644 --- a/deal.II/lac/source/compressed_simple_sparsity_pattern.cc +++ b/deal.II/lac/source/compressed_simple_sparsity_pattern.cc @@ -24,6 +24,155 @@ DEAL_II_NAMESPACE_OPEN + +template +void +CompressedSimpleSparsityPattern::Line::add_entries (ForwardIterator begin, + ForwardIterator end, + const bool indices_are_sorted) +{ + const int n_elements = end - begin; + if (n_elements <= 0) + return; + + const unsigned int n_cols = static_cast(n_elements); + const unsigned int stop_size = entries.size() + n_cols; + + if (indices_are_sorted == true && n_elements > 3) + { + // in debug mode, check whether the + // indices really are sorted. +#ifdef DEBUG + { + ForwardIterator test = begin, test1 = begin; + ++test1; + for ( ; test1 != end; ++test, ++test1) + Assert (*test1 > *test, ExcInternalError()); + } +#endif + + if (entries.size() == 0 || entries.back() < *begin) + { + entries.insert(entries.end(), begin, end); + return; + } + + // resize vector by just inserting the + // list + const unsigned int col = *begin; + std::vector::iterator it = + std::lower_bound(entries.begin(), entries.end(), col); + const unsigned int pos1 = it - entries.begin(); + entries.insert (it, begin, end); + it = entries.begin() + pos1; + + // now merge the two lists. + ForwardIterator my_it = begin; + std::vector::iterator it2 = it + n_cols; + + // as long as there are indices both in + // the end of the entries list and in the + // input list + while (my_it != end && it2 != entries.end()) + { + if (*my_it < *it2) + *it++ = *my_it++; + else if (*my_it == *it2) + { + *it++ = *it2++; + ++my_it; + } + else + *it++ = *it2++; + } + // in case there are indices left in the + // input list + while (my_it != end) + *it++ = *my_it++; + + // in case there are indices left in the + // end of entries + while (it2 != entries.end()) + *it++ = *it2++; + + // resize and return + const unsigned int new_size = it - entries.begin(); + Assert (new_size <= stop_size, ExcInternalError()); + entries.resize (new_size); + return; + } + + // unsorted case or case with too few + // elements + ForwardIterator my_it = begin; + + // If necessary, increase the size of the + // array. + if (stop_size > entries.capacity()) + entries.reserve (stop_size); + + unsigned int col = *my_it; + std::vector::iterator it, it2; + // insert the first element as for one + // entry only first check the last + // element (or if line is still empty) + if ( (entries.size()==0) || ( entries.back() < col) ) { + entries.push_back(col); + it = entries.end()-1; + } + else { + // do a binary search to find the place + // where to insert: + it2 = std::lower_bound(entries.begin(), entries.end(), col); + + // If this entry is a duplicate, continue + // immediately Insert at the right place + // in the vector. Vector grows + // automatically to fit elements. Always + // doubles its size. + if (*it2 != col) + it = entries.insert(it2, col); + else + it = it2; + } + + ++my_it; + // Now try to be smart and insert with + // bias in the direction we are + // walking. This has the advantage that + // for sorted lists, we always search in + // the right direction, what should + // decrease the work needed in here. + for ( ; my_it != end; ++my_it) + { + col = *my_it; + // need a special insertion command when + // we're at the end of the list + if (col > entries.back()) { + entries.push_back(col); + it = entries.end()-1; + } + // search to the right (preferred search + // direction) + else if (col > *it) { + it2 = std::lower_bound(it++, entries.end(), col); + if (*it2 != col) + it = entries.insert(it2, col); + } + // search to the left + else if (col < *it) { + it2 = std::lower_bound(entries.begin(), it, col); + if (*it2 != col) + it = entries.insert(it2, col); + } + // if we're neither larger nor smaller, + // then this was a duplicate and we can + // just continue. + } +} + + + CompressedSimpleSparsityPattern::CompressedSimpleSparsityPattern () : rows(0), @@ -238,4 +387,15 @@ CompressedSimpleSparsityPattern::n_nonzero_elements () const return n; } + +// explicit instantiations +template void CompressedSimpleSparsityPattern::Line::add_entries(unsigned int *, + unsigned int *, + const bool); +template void CompressedSimpleSparsityPattern::Line:: +add_entries(std::vector::iterator, + std::vector::iterator, + const bool); + + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/source/sparse_matrix.inst.in b/deal.II/lac/source/sparse_matrix.inst.in index 9db1b9faee..7766b8c87c 100644 --- a/deal.II/lac/source/sparse_matrix.inst.in +++ b/deal.II/lac/source/sparse_matrix.inst.in @@ -29,6 +29,19 @@ for (S1, S2 : REAL_SCALARS) template void SparseMatrix::add (const S1, const SparseMatrix &); + + template void SparseMatrix::add (const unsigned int, + const unsigned int, + const unsigned int *, + const S2 *, + const bool, + const bool); + + template void SparseMatrix::set (const unsigned int, + const unsigned int, + const unsigned int *, + const S2 *, + const bool); } -- 2.39.5