From: kronbichler Date: Wed, 26 Nov 2008 19:21:53 +0000 (+0000) Subject: The deal.II SparseMatrix class can now set a FullMatrix of data at once, as can all... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=78a5a78c44730d0ad555930b36f2e516df4b66b2;p=dealii-svn.git The deal.II SparseMatrix class can now set a FullMatrix of data at once, as can all block matrices derived from BlockMatrixBase. Other arrays can be set at once as well. git-svn-id: https://svn.dealii.org/trunk@17754 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/block_matrix_base.h b/deal.II/lac/include/lac/block_matrix_base.h index fe842cbff7..b408f40742 100644 --- a/deal.II/lac/include/lac/block_matrix_base.h +++ b/deal.II/lac/include/lac/block_matrix_base.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -477,21 +478,196 @@ class BlockMatrixBase : public Subscriptor void set (const unsigned int i, const unsigned int j, const value_type value); - + + /** + * Set all elements given in a + * FullMatrix into the sparse matrix + * locations given by + * indices. In other words, + * this function writes the elements + * in full_matrix into the + * calling matrix, using the + * local-to-global indexing specified + * by indices for both the + * rows and the columns of the + * matrix. This function assumes a + * quadratic sparse matrix and a + * quadratic full_matrix, the usual + * situation in FE calculations. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be set anyway or + * they should be filtered away (and + * not change the previous content in + * the respective element if it + * exists). The default value is + * false, i.e., even zero + * values are treated. + */ + void set (const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Same function as before, but now + * including the possibility to use + * rectangular full_matrices and + * different local-to-global indexing + * on rows and columns, respectively. + */ + void set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be set anyway or + * they should be filtered away (and + * not change the previous content in + * the respective element if it + * exists). The default value is + * false, i.e., even zero + * values are treated. + */ + void set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = false); + + /** + * Set several elements to values + * given by values in a + * given row in columns given by + * col_indices into the sparse + * matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be inserted anyway + * or they should be filtered + * away. The default value is + * false, i.e., even zero + * values are inserted/replaced. + */ + void set (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const value_type *values, + const bool elide_zero_values = false); + /** * Add value to the - * element (i,j). - * Throws an error if the entry - * does not exist or if - * value is not a finite - * number. Still, it is allowed - * to store zero values in - * non-existent fields. + * element (i,j). Throws + * an error if the entry does not + * exist or if value is + * not a finite number. Still, it + * is allowed to store zero + * values in non-existent fields. */ void add (const unsigned int i, const unsigned int j, const value_type value); + /** + * Add all elements given in a + * FullMatrix into sparse + * matrix locations given by + * indices. In other words, + * this function adds the elements in + * full_matrix to the + * respective entries in calling + * matrix, using the local-to-global + * indexing specified by + * indices for both the rows + * and the columns of the + * matrix. This function assumes a + * quadratic sparse matrix and a + * quadratic full_matrix, the usual + * situation in FE calculations. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + void add (const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = true); + + /** + * Same function as before, but now + * including the possibility to use + * rectangular full_matrices and + * different local-to-global indexing + * on rows and columns, respectively. + */ + void add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = true); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + void add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = true); + + /** + * Add an array of values given by + * values in the given + * global matrix row at columns + * specified by col_indices in the + * sparse matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + void add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const value_type *values, + const bool elide_zero_values = true); + /** * Return the value of the entry * (i,j). This may be an @@ -984,7 +1160,30 @@ class BlockMatrixBase : public Subscriptor friend class ConstAccessor; friend class const_iterator; -#endif +#endif + + + private: + /** + * Temporary vector for block indices + * when writing local to global data. + */ + std::vector block_col_indices; + + /** + * Temporary vector for column + * indices when writing local to + * global data and passing the call + * down to a SparseMatrix. + */ + std::vector local_col_indices; + + /** + * Temporary vector for storing the + * local lengths of the row data. + */ + std::vector local_row_length; + }; @@ -1554,6 +1753,9 @@ BlockMatrixBase::n_block_rows () const + // Write the single set manually, + // since the other function has a lot + // of overhead in that case. template inline void @@ -1575,6 +1777,144 @@ BlockMatrixBase::set (const unsigned int i, +template +inline +void +BlockMatrixBase::set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (row_indices.size() == values.m(), + ExcDimensionMismatch(row_indices.size(), values.m())); + Assert (col_indices.size() == values.n(), + ExcDimensionMismatch(col_indices.size(), values.n())); + + for (unsigned int i=0; i +inline +void +BlockMatrixBase::set (const std::vector &indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (indices.size() == values.m(), + ExcDimensionMismatch(indices.size(), values.m())); + Assert (values.n() == values.m(), ExcNotQuadratic()); + + for (unsigned int i=0; i +inline +void +BlockMatrixBase::set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) +{ + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + set (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); +} + + + + // This is a very messy function, since + // we need to calculate to each position + // the location in the global array. +template +inline +void +BlockMatrixBase::set (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const value_type *values, + const bool elide_zero_values) +{ + // Resize scratch arrays + local_row_length.resize (this->n_block_cols()); + block_col_indices.resize (this->n_block_cols()); + local_col_indices.resize (n_cols); + + // Clear the content in local_row_length + for (unsigned int i=0; in_block_cols(); ++i) + local_row_length[i] = 0; + + // Go through the column indices to find + // out which portions of the values + // should be written into which block + // matrix. + { + unsigned int current_block = 0, row_length = 0; + block_col_indices[0] = 0; + for (unsigned int j=0; j + col_index = this->column_block_indices.global_to_local(col_indices[j]); + local_col_indices[j] = col_index.second; + if (col_index.first > current_block) + { + local_row_length[current_block] = row_length; + row_length = 0; + while (col_index.first > current_block) + current_block++; + block_col_indices[current_block] = j; + } + + Assert (col_index.first == current_block, + ExcInternalError()); + } + local_row_length[current_block] = row_length; + Assert (current_block < this->n_block_cols(), + ExcInternalError()); + +#ifdef DEBUG + // If in debug mode, do a check whether + // the right length has been obtained. + unsigned int length = 0; + for (unsigned int i=0; in_block_cols(); ++i) + length += local_row_length[i]; + Assert (length == n_cols, + ExcDimensionMismatch(length, n_cols)); +#endif + } + + // Now we found out about where the + // individual columns should start and + // where we should start reading out + // data. Now let's write the data into + // the individual blocks! + const std::pair + row_index = this->row_block_indices.global_to_local (row); + for (unsigned int block_col=0; block_col inline void @@ -1605,6 +1945,148 @@ BlockMatrixBase::add (const unsigned int i, +template +inline +void +BlockMatrixBase::add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (row_indices.size() == values.m(), + ExcDimensionMismatch(row_indices.size(), values.m())); + Assert (col_indices.size() == values.n(), + ExcDimensionMismatch(col_indices.size(), values.n())); + + for (unsigned int i=0; i +inline +void +BlockMatrixBase::add (const std::vector &indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (indices.size() == values.m(), + ExcDimensionMismatch(indices.size(), values.m())); + Assert (values.n() == values.m(), ExcNotQuadratic()); + + for (unsigned int i=0; i +inline +void +BlockMatrixBase::add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) +{ + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + add (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); +} + + + + // This is a very messy function, since + // we need to calculate to each position + // the location in the global array. +template +inline +void +BlockMatrixBase::add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const value_type *values, + const bool elide_zero_values) +{ + // TODO: Look over this to find out + // whether we can do that more + // efficiently. + + // Resize scratch arrays + block_col_indices.resize (this->n_block_cols()); + local_row_length.resize (this->n_block_cols()); + local_col_indices.resize (n_cols); + + // Clear the content in local_row_length + for (unsigned int i=0; in_block_cols(); ++i) + local_row_length[i] = 0; + + // Go through the column indices to find + // out which portions of the values + // should be written into which block + // matrix. + { + unsigned int current_block = 0, row_length = 0; + block_col_indices[0] = 0; + for (unsigned int j=0; j + col_index = this->column_block_indices.global_to_local(col_indices[j]); + local_col_indices[j] = col_index.second; + if (col_index.first > current_block) + { + local_row_length[current_block] = row_length; + row_length = 0; + while (col_index.first > current_block) + current_block++; + block_col_indices[current_block] = j; + } + + Assert (col_index.first == current_block, + ExcInternalError()); + } + local_row_length[current_block] = row_length; + Assert (current_block < this->n_block_cols(), + ExcInternalError()); + +#ifdef DEBUG + // If in debug mode, do a check whether + // the right length has been obtained. + unsigned int length = 0; + for (unsigned int i=0; in_block_cols(); ++i) + length += local_row_length[i]; + Assert (length == n_cols, + ExcDimensionMismatch(length, n_cols)); +#endif + } + + // Now we found out about where the + // individual columns should start and + // where we should start reading out + // data. Now let's write the data into + // the individual blocks! + const std::pair + row_index = this->row_block_indices.global_to_local (row); + for (unsigned int block_col=0; block_col inline typename BlockMatrixBase::value_type diff --git a/deal.II/lac/include/lac/petsc_matrix_base.h b/deal.II/lac/include/lac/petsc_matrix_base.h index 7436e1b037..14c9205080 100644 --- a/deal.II/lac/include/lac/petsc_matrix_base.h +++ b/deal.II/lac/include/lac/petsc_matrix_base.h @@ -1444,8 +1444,7 @@ namespace PETScWrappers else { // Otherwise, extract nonzero values in - // each row and pass on to the other - // function. + // each row and get the respective index. column_indices.resize(n_cols); column_values.resize(n_cols); @@ -1609,8 +1608,7 @@ namespace PETScWrappers else { // Otherwise, extract nonzero values in - // each row and pass on to the other - // function. + // each row and get the respective index. column_indices.resize(n_cols); column_values.resize(n_cols); diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index da0a9c3147..16f04091b7 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -765,7 +765,95 @@ class SparseMatrix : public virtual Subscriptor void set (const unsigned int i, const unsigned int j, const number value); - + + /** + * Set all elements given in a + * FullMatrix into the sparse matrix + * locations given by + * indices. In other words, + * this function writes the elements + * in full_matrix into the + * calling matrix, using the + * local-to-global indexing specified + * by indices for both the + * rows and the columns of the + * matrix. This function assumes a + * quadratic sparse matrix and a + * quadratic full_matrix, the usual + * situation in FE calculations. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be set anyway or + * they should be filtered away (and + * not change the previous content in + * the respective element if it + * exists). The default value is + * false, i.e., even zero + * values are treated. + */ + void set (const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Same function as before, but now + * including the possibility to use + * rectangular full_matrices and + * different local-to-global indexing + * on rows and columns, respectively. + */ + void set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be set anyway or + * they should be filtered away (and + * not change the previous content in + * the respective element if it + * exists). The default value is + * false, i.e., even zero + * values are treated. + */ + void set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = false); + + /** + * Set several elements to values + * given by values in a + * given row in columns given by + * col_indices into the sparse + * matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be inserted anyway + * or they should be filtered + * away. The default value is + * false, i.e., even zero + * values are inserted/replaced. + */ + void set (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool elide_zero_values = false); + /** * Add value to the * element (i,j). Throws @@ -779,6 +867,94 @@ class SparseMatrix : public virtual Subscriptor const unsigned int j, const number value); + /** + * Add all elements given in a + * FullMatrix into sparse + * matrix locations given by + * indices. In other words, + * this function adds the elements in + * full_matrix to the + * respective entries in calling + * matrix, using the local-to-global + * indexing specified by + * indices for both the rows + * and the columns of the + * matrix. This function assumes a + * quadratic sparse matrix and a + * quadratic full_matrix, the usual + * situation in FE calculations. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + void add (const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = true); + + /** + * Same function as before, but now + * including the possibility to use + * rectangular full_matrices and + * different local-to-global indexing + * on rows and columns, respectively. + */ + void add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = true); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + void add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = true); + + /** + * Add an array of values given by + * values in the given + * global matrix row at columns + * specified by col_indices in the + * sparse matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + void add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool elide_zero_values = true); + /** * Multiply the entire matrix by a * fixed factor. @@ -1713,6 +1889,26 @@ 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 @@ -1812,51 +2008,290 @@ unsigned int SparseMatrix::n () const } -template + // 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. +template +inline +void +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); +} + + + +template +inline +void +SparseMatrix::set (const std::vector &indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (indices.size() == values.m(), + ExcDimensionMismatch(indices.size(), values.m())); + Assert (values.m() == values.n(), ExcNotQuadratic()); + + for (unsigned int i=0; i +inline +void +SparseMatrix::set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (row_indices.size() == values.m(), + ExcDimensionMismatch(row_indices.size(), values.m())); + Assert (col_indices.size() == values.n(), + ExcDimensionMismatch(col_indices.size(), values.n())); + + for (unsigned int i=0; i inline -void SparseMatrix::set (const unsigned int i, - const unsigned int j, - const number value) +void +SparseMatrix::set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) { + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + set (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); +} + - Assert (numbers::is_finite(value), - ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); +template +inline +void +SparseMatrix::set (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool elide_zero_values) +{ Assert (cols != 0, ExcNotInitialized()); - // it is allowed to set elements of + + 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); + + const unsigned int * index_ptr = &col_indices[0]; + const number * value_ptr = &values[0]; + + // 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, *index_ptr++); + + // 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 - const unsigned int index = cols->operator()(i,j); - Assert ((index != SparsityPattern::invalid_entry) || - (value == 0.), - ExcInvalidIndex(i,j)); + 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++; + } - if (index != SparsityPattern::invalid_entry) - val[index] = value; + index_ptr = &global_indices[0]; + 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 +inline +void +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); +} + + + +template +inline +void +SparseMatrix::add (const std::vector &indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (indices.size() == values.m(), + ExcDimensionMismatch(indices.size(), values.m())); + Assert (values.m() == values.n(), ExcNotQuadratic()); + + for (unsigned int i=0; i +inline +void +SparseMatrix::add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &values, + const bool elide_zero_values) +{ + Assert (row_indices.size() == values.m(), + ExcDimensionMismatch(row_indices.size(), values.m())); + Assert (col_indices.size() == values.n(), + ExcDimensionMismatch(col_indices.size(), values.n())); + + for (unsigned int i=0; i inline -void SparseMatrix::add (const unsigned int i, - const unsigned int j, - const number value) +void +SparseMatrix::add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) { + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + add (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); +} + - Assert (numbers::is_finite(value), - ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); +template +inline +void +SparseMatrix::add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number *values, + const bool elide_zero_values) +{ Assert (cols != 0, ExcNotInitialized()); - const unsigned int index = cols->operator()(i,j); - Assert ((index != SparsityPattern::invalid_entry) || - (value == 0.), - ExcInvalidIndex(i,j)); + unsigned int n_columns = 0; - if (value != 0.) - val[index] += value; + global_indices.resize(n_cols); + column_values.resize(n_cols); + + const unsigned int * index_ptr = &col_indices[0]; + const number * value_ptr = &values[0]; + + // First, search all the indices to find + // out which values we actually need to + // add. + // TODO: Could probably made more + // efficient. + for (unsigned int j=0; joperator()(row, *index_ptr++); + + // 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++; + } + + index_ptr = &global_indices[0]; + value_ptr = &column_values[0]; + + // Finally, go through the index list + // and add the elements one by one. + for (unsigned int j=0; ji,j) - * to @p value. - */ - void set (const unsigned int i, - const unsigned int j, - const TrilinosScalar value); - - /** - * Set all elements given in a - * FullMatrix into the sparse - * matrix, according to row_indices - * and col_indices. - * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be inserted anyway - * or they should be filtered - * away. The default value is - * false, i.e., even zero - * values are inserted/replaced. - */ - void set (const std::vector &row_indices, - const std::vector &col_indices, - const FullMatrix &full_matrix, - const bool elide_zero_values = false); - - /** - * Set several elements in the - * specified row of the matrix with - * column indices as given by - * col_indices to the - * respective value. - * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be inserted anyway - * or they should be filtered - * away. The default value is - * false, i.e., even zero - * values are inserted/replaced. - */ - void set (const unsigned int row, - const std::vector &col_indices, - const std::vector &values, - const bool elide_zero_values = false); - - /** - * Set several elements to values - * given by values in the - * given row at column locations - * specified by col_indices in the - * sparse matrix. - * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be inserted anyway - * or they should be filtered - * away. The default value is - * false, i.e., even zero - * values are inserted/replaced. - */ - void set (const unsigned int row, - const unsigned int n_cols, - const unsigned int *col_indices, - const TrilinosScalar *values, - const bool elide_zero_values = false); - - /** - * Add @p value to the element - * (i,j). - */ - void add (const unsigned int i, - const unsigned int j, - const TrilinosScalar value); - - /** - * Add all elements in a FullMatrix - * full_matrix (cell matrix) - * into the sparse matrix, at - * locations according to row_indices - * and col_indices. - * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be added anyway or - * these should be filtered away and - * only non-zero data is added. The - * default value is true, - * i.e., zero values won't be added - * into the matrix. - */ - void add (const std::vector &row_indices, - const std::vector &col_indices, - const FullMatrix &full_matrix, - const bool elide_zero_values = true); - - /** - * Set several elements in the - * specified row of the matrix with - * column indices as given by - * col_indices to the - * respective value. - * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be added anyway or - * these should be filtered away and - * only non-zero data is added. The - * default value is true, - * i.e., zero values won't be added - * into the matrix. - */ - void add (const unsigned int row, - const std::vector &col_indices, - const std::vector &values, - const bool elide_zero_values = true); - - /** - * Set several elements to values - * given by values in the - * given row at column locations - * specified by col_indices in the - * sparse matrix. - * - * The optional parameter - * elide_zero_values can be - * used to specify whether zero - * values should be added anyway or - * these should be filtered away and - * only non-zero data is added. The - * default value is true, - * i.e., zero values won't be added - * into the matrix. - */ - void add (const unsigned int row, - const unsigned int n_cols, - const unsigned int *col_indices, - const TrilinosScalar *values, - const bool elide_zero_values = true); - /** * Matrix-vector multiplication: * let $dst = M*src$ with $M$ @@ -770,14 +624,6 @@ namespace TrilinosWrappers << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3 << ',' << arg4 << "] have differing column numbers."); ///@} - - private: - /** - * A few temporary vectors for - * writing data from local to global - * matrices. - */ - std::vector block_col_indices, local_col_indices, local_row_length; }; @@ -804,267 +650,6 @@ namespace TrilinosWrappers - inline - void - BlockSparseMatrix::set (const unsigned int i, - const unsigned int j, - const TrilinosScalar value) - { - BaseClass::set (i, j, value); - } - - - - inline - void - BlockSparseMatrix::set (const std::vector &row_indices, - const std::vector &col_indices, - const FullMatrix &values, - const bool elide_zero_values) - { - Assert (row_indices.size() == values.m(), - ExcDimensionMismatch(row_indices.size(), values.m())); - Assert (col_indices.size() == values.n(), - ExcDimensionMismatch(col_indices.size(), values.n())); - - for (unsigned int i=0; i &col_indices, - const std::vector &values, - const bool elide_zero_values) - { - Assert (col_indices.size() == values.size(), - ExcDimensionMismatch(col_indices.size(), values.size())); - - set (row, col_indices.size(), &col_indices[0], &values[0], - elide_zero_values); - } - - - - // This is a very messy function, since - // we need to calculate to each position - // the location in the global array. - inline - void - BlockSparseMatrix::set (const unsigned int row, - const unsigned int n_cols, - const unsigned int *col_indices, - const TrilinosScalar *values, - const bool elide_zero_values) - { - // Resize scratch arrays - local_row_length.resize (this->n_block_cols()); - block_col_indices.resize (this->n_block_cols()); - local_col_indices.resize (n_cols); - - // Clear the content in local_row_length - for (unsigned int i=0; in_block_cols(); ++i) - local_row_length[i] = 0; - - // Go through the column indices to find - // out which portions of the values - // should be written into which block - // matrix. - { - unsigned int current_block = 0, row_length = 0; - block_col_indices[0] = 0; - for (unsigned int j=0; j - col_index = this->column_block_indices.global_to_local(col_indices[j]); - local_col_indices[j] = col_index.second; - if (col_index.first > current_block) - { - local_row_length[current_block] = row_length; - row_length = 0; - while (col_index.first > current_block) - current_block++; - block_col_indices[current_block] = j; - } - - Assert (col_index.first == current_block, - ExcInternalError()); - } - local_row_length[current_block] = row_length; - Assert (current_block < this->n_block_cols(), - ExcInternalError()); - -#ifdef DEBUG - // If in debug mode, do a check whether - // the right length has been obtained. - unsigned int length = 0; - for (unsigned int i=0; in_block_cols(); ++i) - length += local_row_length[i]; - Assert (length == n_cols, - ExcDimensionMismatch(length, n_cols)); -#endif - } - - // Now we found out about where the - // individual columns should start and - // where we should start reading out - // data. Now let's write the data into - // the individual blocks! - const std::pair - row_index = this->row_block_indices.global_to_local (row); - for (unsigned int block_col=0; block_col &row_indices, - const std::vector &col_indices, - const FullMatrix &values, - const bool elide_zero_values) - { - Assert (row_indices.size() == values.m(), - ExcDimensionMismatch(row_indices.size(), values.m())); - Assert (col_indices.size() == values.n(), - ExcDimensionMismatch(col_indices.size(), values.n())); - - for (unsigned int i=0; i &col_indices, - const std::vector &values, - const bool elide_zero_values) - { - Assert (col_indices.size() == values.size(), - ExcDimensionMismatch(col_indices.size(), values.size())); - - add (row, col_indices.size(), &col_indices[0], &values[0], - elide_zero_values); - } - - - - inline - void - BlockSparseMatrix::add (const unsigned int row, - const unsigned int n_cols, - const unsigned int *col_indices, - const TrilinosScalar *values, - const bool elide_zero_values) - { - // TODO: Look over this to find out - // whether we can do that more - // efficiently. - - // Resize scratch arrays - block_col_indices.resize (this->n_block_cols()); - local_row_length.resize (this->n_block_cols()); - local_col_indices.resize (n_cols); - - // Clear the content in local_row_length - for (unsigned int i=0; in_block_cols(); ++i) - local_row_length[i] = 0; - - // Go through the column indices to find - // out which portions of the values - // should be written into which block - // matrix. - { - unsigned int current_block = 0, row_length = 0; - block_col_indices[0] = 0; - for (unsigned int j=0; j - col_index = this->column_block_indices.global_to_local(col_indices[j]); - local_col_indices[j] = col_index.second; - if (col_index.first > current_block) - { - local_row_length[current_block] = row_length; - row_length = 0; - while (col_index.first > current_block) - current_block++; - block_col_indices[current_block] = j; - } - - Assert (col_index.first == current_block, - ExcInternalError()); - } - local_row_length[current_block] = row_length; - Assert (current_block < this->n_block_cols(), - ExcInternalError()); - -#ifdef DEBUG - // If in debug mode, do a check whether - // the right length has been obtained. - unsigned int length = 0; - for (unsigned int i=0; in_block_cols(); ++i) - length += local_row_length[i]; - Assert (length == n_cols, - ExcDimensionMismatch(length, n_cols)); -#endif - } - - // Now we found out about where the - // individual columns should start and - // where we should start reading out - // data. Now let's write the data into - // the individual blocks! - const std::pair - row_index = this->row_block_indices.global_to_local (row); - for (unsigned int block_col=0; block_col