From ba1bd34ba15cd643770704965bb96fcbe6fe691f Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 16 Apr 2009 11:28:42 +0000 Subject: [PATCH] Make the sparse matrix collective add function a bit smarter. Avoid having some additional arrays. git-svn-id: https://svn.dealii.org/trunk@18619 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_matrix.h | 22 +--- .../lac/include/lac/sparse_matrix.templates.h | 121 +++++++++--------- 2 files changed, 66 insertions(+), 77 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 76ddb7edae..b6d3480613 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -1898,26 +1898,6 @@ class SparseMatrix : public virtual Subscriptor */ unsigned int max_len; - /** - * An internal array of integer - * values that is used to store the - * indices in the global value array - * val when adding/inserting local - * data into the (large) sparse - * matrix. - */ - std::vector global_indices; - - /** - * An internal array of double values - * that is used to store the column - * indices when adding/inserting - * local data into the (large) sparse - * matrix. - */ - std::vector column_values; - - /** * Version of vmult() which only * performs its actions on the @@ -1990,6 +1970,8 @@ class SparseMatrix : public virtual Subscriptor // make all other sparse matrices // friends template friend class SparseMatrix; + template friend class SparseLUDecomposition; + template friend class SparseILU; }; /*@}*/ diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index f89a70a952..70e32b125e 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -342,7 +342,7 @@ SparseMatrix::add (const unsigned int row, // 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 + // doing many binary searches if (col_indices_are_sorted == true) if (n_cols * 8 > cols->row_length(row)) { @@ -359,17 +359,37 @@ SparseMatrix::add (const unsigned int row, &cols->get_column_numbers()[cols->get_rowstart_indices()[row]]; number * val_ptr = &val[cols->get_rowstart_indices()[row]]; Assert (this_cols[0] == row, ExcInternalError()); + + // find diagonal and add it if found + const unsigned int * diag_pos = std::lower_bound (col_indices, + col_indices+n_cols, + row); + const unsigned int diag = diag_pos - col_indices; + unsigned int post_diag = diag; + if (diag != n_cols && *diag_pos == row) + { + val_ptr[0] += *(values+(diag_pos-col_indices)); + ++post_diag; + } + + // add indices before diagonal 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])); + + val_ptr[counter] += values[i]; + } + + // add indices after diagonal + for (unsigned int i=post_diag; i= this_cols[counter], ExcInternalError()); while (this_cols[counter] < col_indices[i]) @@ -380,6 +400,7 @@ SparseMatrix::add (const unsigned int row, val_ptr[counter] += values[i]; } + Assert (counter < cols->row_length(row), ExcInternalError()); } else @@ -405,17 +426,11 @@ SparseMatrix::add (const unsigned int row, 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. + // unsorted case: first, search all the + // indices to find out which values we + // actually need to add. + const unsigned int * const my_cols = cols->get_column_numbers(); + unsigned int index = cols->get_rowstart_indices()[row]; for (unsigned int j=0; j::add (const unsigned int row, if (value == 0 && elide_zero_values == true) continue; - const unsigned int index = cols->operator()(row, col_indices[j]); + // check whether the next index to add is + // the next present index in the sparsity + // pattern (otherwise, do a binary + // search) + if (index != cols->get_rowstart_indices()[row+1] && + my_cols[index] == col_indices[j]) + goto add_value; + + 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 + // value we add is zero if (index == SparsityPattern::invalid_entry) { Assert ((index != SparsityPattern::invalid_entry) || @@ -440,18 +463,10 @@ SparseMatrix::add (const unsigned int row, continue; } - global_indices[n_columns] = index; - column_values[n_columns] = value; - n_columns++; + add_value: + val[index] += value; + ++index; } - - 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::set (const unsigned int row, { 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. + const unsigned int * my_cols = cols->get_column_numbers(); + unsigned int index = cols->get_rowstart_indices()[row]; for (unsigned int j=0; j::set (const unsigned int row, if (value == 0 && elide_zero_values == true) continue; - const unsigned int index = cols->operator()(row, col_indices[j]); + // check whether the next index to set is + // the next present index in the sparsity + // pattern (otherwise, do a binary + // search) + if (index != cols->get_rowstart_indices()[row+1] && + my_cols[index] == col_indices[j]) + goto set_value; + + index = cols->operator()(row, col_indices[j]); // it is allowed to set elements in // the matrix that are not part of @@ -502,18 +519,10 @@ SparseMatrix::set (const unsigned int row, continue; } - global_indices[n_columns] = index; - column_values[n_columns] = value; - n_columns++; + set_value: + val[index] = value; + ++index; } - - 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 unsigned int SparseMatrix::memory_consumption () const { - return sizeof(*this) + max_len*sizeof(number) + - global_indices.size()*sizeof(unsigned int) + - column_values.size()*sizeof(number); + return sizeof(*this) + max_len*sizeof(number); } -- 2.39.5