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
-template <typename ForwardIterator>
-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<unsigned int>(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<unsigned int>::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<unsigned int>::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<unsigned int>::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
// 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 <typename number>
inline
void
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;
}
-template <typename number>
-template <typename number2>
-inline
-void
-SparseMatrix<number>::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; j<n_cols; ++j)
- {
- const number value = values[j];
- Assert (numbers::is_finite(value),
- ExcMessage("The given value is not finite but either "
- "infinite or Not A Number (NaN)"));
-
- if (value == 0 && elide_zero_values == true)
- continue;
-
- const unsigned int index = cols->operator()(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<n_columns; ++j)
- val[*index_ptr++] = *value_ptr++;
-}
-
-
-
template <typename number>
inline
void
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;
}
-template <typename number>
-template <typename number2>
-inline
-void
-SparseMatrix<number>::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<n_cols; ++i)
- Assert (col_indices[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<n_cols; ++i)
- {
- // diagonal
- if (col_indices[i] == row)
- {
- if (values[i] != 0)
- val_ptr[0] += values[i];
- continue;
- }
-
- Assert (col_indices[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<n_cols; ++i)
- {
- Assert (col_indices[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; j<n_cols; ++j)
- {
- const number value = values[j];
- Assert (numbers::is_finite(value),
- ExcMessage("The given value is not finite but either "
- "infinite or Not A Number (NaN)"));
-
- if (value == 0 && elide_zero_values == true)
- continue;
-
- const unsigned int index = cols->operator()(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<n_columns; ++j)
- val[*index_ptr++] += *value_ptr++;
-}
-
-
-
template <typename number>
inline
SparseMatrix<number> &
}
+
+template <typename number>
+template <typename number2>
+void
+SparseMatrix<number>::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<n_cols; ++i)
+ Assert (col_indices[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<n_cols; ++i)
+ {
+ // diagonal
+ if (col_indices[i] == row)
+ {
+ if (values[i] != 0)
+ val_ptr[0] += values[i];
+ continue;
+ }
+
+ Assert (col_indices[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<n_cols; ++i)
+ {
+ Assert (col_indices[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; j<n_cols; ++j)
+ {
+ const number value = values[j];
+ Assert (numbers::is_finite(value),
+ ExcMessage("The given value is not finite but either "
+ "infinite or Not A Number (NaN)"));
+
+ if (value == 0 && elide_zero_values == true)
+ continue;
+
+ const unsigned int index = cols->operator()(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<n_columns; ++j)
+ val[*index_ptr++] += *value_ptr++;
+}
+
+
+
+template <typename number>
+template <typename number2>
+void
+SparseMatrix<number>::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; j<n_cols; ++j)
+ {
+ const number value = values[j];
+ Assert (numbers::is_finite(value),
+ ExcMessage("The given value is not finite but either "
+ "infinite or Not A Number (NaN)"));
+
+ if (value == 0 && elide_zero_values == true)
+ continue;
+
+ const unsigned int index = cols->operator()(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<n_columns; ++j)
+ val[*index_ptr++] = *value_ptr++;
+}
+
+
+
template <typename number>
template <class OutVector, class InVector>
void
DEAL_II_NAMESPACE_OPEN
+
+template <typename ForwardIterator>
+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<unsigned int>(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<unsigned int>::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<unsigned int>::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<unsigned int>::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),
return n;
}
+
+// explicit instantiations
+template void CompressedSimpleSparsityPattern::Line::add_entries(unsigned int *,
+ unsigned int *,
+ const bool);
+template void CompressedSimpleSparsityPattern::Line::
+add_entries(std::vector<unsigned int>::iterator,
+ std::vector<unsigned int>::iterator,
+ const bool);
+
+
DEAL_II_NAMESPACE_CLOSE
template void SparseMatrix<S1>::add<S2> (const S1,
const SparseMatrix<S2> &);
+
+ template void SparseMatrix<S1>::add<S2> (const unsigned int,
+ const unsigned int,
+ const unsigned int *,
+ const S2 *,
+ const bool,
+ const bool);
+
+ template void SparseMatrix<S1>::set<S2> (const unsigned int,
+ const unsigned int,
+ const unsigned int *,
+ const S2 *,
+ const bool);
}
projects / dealii-svn.git / commitdiff