From: kronbichler Date: Fri, 21 Nov 2008 18:33:29 +0000 (+0000) Subject: Introduced a few add() and set() functions for Trilinos matrices that allow to set... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=16d7e171809622fcf51f3f90681fb1a6fc8f0fd1;p=dealii-svn.git Introduced a few add() and set() functions for Trilinos matrices that allow to set several elements at once. We will need to rewrite ConstraintMatrix::distribute_local_to_global in order to use this new capability. git-svn-id: https://svn.dealii.org/trunk@17678 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h index 90133e1bdc..0032a433dc 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -291,6 +292,95 @@ namespace TrilinosWrappers */ unsigned int n_nonzero_elements () const; + /** + * Set the element (i,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. + */ + void set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + */ + void set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values); + + /** + * Set several elements to values + * given by values at + * locations specified by row_indices + * and col_indices in the sparse + * matrix. The array values + * is assumed to store data in C + * style, i.e., row-wise. + */ + void set (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values); + + /** + * 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. + */ + void add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + */ + void add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values); + + /** + * Set several elements to values + * given by values at + * locations specified by row_indices + * and col_indices in the sparse + * matrix. The array values + * is assumed to store data in C + * style, i.e., row-wise. + */ + void add (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values); + /** * Matrix-vector multiplication: * let $dst = M*src$ with $M$ @@ -623,6 +713,14 @@ 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_row_length, local_col_indices; }; diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index e06d4ef268..6f66f4233f 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -866,6 +867,75 @@ namespace TrilinosWrappers 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. + * + * This function is able to insert + * new elements into the matrix as + * long as compress() has not been + * called, so the sparsity pattern + * will be extended. When compress() + * is called for the first time, then + * this is no longer possible and an + * insertion of elements at positions + * which have not been initialized + * will throw an exception. + */ + void set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * This function is able to insert + * new elements into the matrix as + * long as compress() has not been + * called, so the sparsity pattern + * will be extended. When compress() + * is called for the first time, then + * this is no longer possible and an + * insertion of elements at positions + * which have not been initialized + * will throw an exception. + */ + void set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values); + + /** + * Set several elements to values + * given by values at + * locations specified by row_indices + * and col_indices in the sparse + * matrix. The array values + * is assumed to store data in C + * style, i.e., row-wise. + * + * This function is able to insert + * new elements into the matrix as + * long as compress() has not been + * called, so the sparsity pattern + * will be extended. When compress() + * is called for the first time, then + * this is no longer possible and an + * insertion of elements at positions + * which have not been initialized + * will throw an exception. + */ + void set (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values); + /** * Add @p value to the element * (i,j). @@ -884,7 +954,71 @@ namespace TrilinosWrappers 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. + * + * Just as the respective call in + * deal.II SparseMatrix + * class (but in contrast to the + * situation for PETSc based + * matrices), this function + * throws an exception if an + * entry does not exist in the + * sparsity pattern. + */ + void add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix); + + /** + * Set several elements in the + * specified row of the matrix with + * column indices as given by + * col_indices to the + * respective value. + * + * Just as the respective call in + * deal.II SparseMatrix + * class (but in contrast to the + * situation for PETSc based + * matrices), this function + * throws an exception if an + * entry does not exist in the + * sparsity pattern. + */ + void add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values); + + /** + * Set several elements to values + * given by values at + * locations specified by row_indices + * and col_indices in the sparse + * matrix. The array values + * is assumed to store data in C + * style, i.e., row-wise. + * + * Just as the respective call in + * deal.II SparseMatrix + * class (but in contrast to the + * situation for PETSc based + * matrices), this function + * throws an exception if an + * entry does not exist in the + * sparsity pattern. + */ + void add (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values); + /** * Multiply the entire matrix * by a fixed factor. diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc index 8fc6101880..55169ba569 100644 --- a/deal.II/lac/source/trilinos_block_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_block_sparse_matrix.cc @@ -262,6 +262,258 @@ namespace TrilinosWrappers + void + BlockSparseMatrix::set (const unsigned int i, + const unsigned int j, + const TrilinosScalar value) + { + BaseClass::set (i, j, value); + } + + + + void + BlockSparseMatrix::set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &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())); + + set (row_indices.size(), &row_indices[0], + col_indices.size(), &col_indices[0], &values(0,0)); + } + + + + void + BlockSparseMatrix::set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values) + { + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + set (1, &row, col_indices.size(), &col_indices[0], &values[0]); + } + + + + void + BlockSparseMatrix::set (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values) + { + // 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. This can be done before + // starting the loop over all the rows, + // since we assume a rectangular set of + // matrix data. + { + 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! + for (unsigned int i=0; i + row_index = this->row_block_indices.global_to_local (row_indices[i]); + for (unsigned int block_col=0; block_col &row_indices, + const std::vector &col_indices, + const FullMatrix &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())); + + add (row_indices.size(), &row_indices[0], + col_indices.size(), &col_indices[0], &values(0,0)); + } + + + + void + BlockSparseMatrix::add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values) + { + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + add (1, &row, col_indices.size(), &col_indices[0], &values[0]); + } + + + + void + BlockSparseMatrix::add (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *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. This can be done before + // starting the loop over all the rows, + // since we assume a rectangular set of + // matrix data. + { + 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! + for (unsigned int i=0; i + row_index = this->row_block_indices.global_to_local (row_indices[i]); + for (unsigned int block_col=0; block_col &row_indices, + const std::vector &col_indices, + const FullMatrix &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())); + + set (row_indices.size(), &row_indices[0], + col_indices.size(), &col_indices[0], &values(0,0)); + } + + + + void + SparseMatrix::set (const unsigned int row, + const std::vector &col_indices, + const std::vector &values) + { + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + set (1, &row, col_indices.size(), &col_indices[0], &values[0]); + } + + + + void + SparseMatrix::set (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values) + { + int ierr; if (last_action == Add) { - int ierr; ierr = matrix->GlobalAssemble(col_map, row_map, false); - + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } - // If our last action was not an - // addition, we tell that at this - // point here. if (last_action != Insert) last_action = Insert; - compressed = false; - - int ierr; - - // In case this matrix owns the - // actual row, we can save some - // computing time by directly - // accessing the Epetra_CrsMatrix - // set functions. - if (row_map.MyGID(i) == true) + // Now go through all rows that are + // present in the input data. + for (unsigned int i=0; iFilled() == false) + // If the calling matrix owns the row to + // which we want to add values, we + // can directly call the Epetra_CrsMatrix + // input function, which is much faster + // than the Epetra_FECrsMatrix function. + if (row_map.MyGID(i) == true) { - ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues ( - (int)i, 1, - const_cast(&value), - (int*)&j); + if (matrix->Filled() == false) + { + ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues( + (int)row_indices[i], (int)n_cols, + const_cast(&values[i*n_cols]), + (int*)&col_indices[0]); // When adding up elements, we do // not want to create exceptions in // the case when adding elements. - if (ierr > 0) - ierr = 0; + if (ierr > 0) + ierr = 0; + } + else + ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues( + (int)row_indices[i], (int)n_cols, + const_cast(&values[i*n_cols]), + (int*)&col_indices[0]); } else { - ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues ( - (int)i, 1, - const_cast(&value), - (int*)&j); + // When we're at off-processor data, we + // have to stick with the standard + // SumIntoGlobalValues + // function. Nevertheless, the way we + // call it is the fastest one (any other + // will lead to repeated allocation and + // deallocation of memory in order to + // call the function we already use, + // which is very unefficient if writing + // one element at a time). + compressed = false; + + const TrilinosScalar* value_ptr = &values[i*n_cols]; + + if (matrix->Filled() == false) + { + ierr = matrix->InsertGlobalValues (1, (int*)&i, + (int)n_cols, (int*)&col_indices[0], + &value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); + if (ierr > 0) + ierr = 0; + } + else + ierr = matrix->ReplaceGlobalValues (1, (int*)&i, + (int)n_cols, (int*)&col_indices[0], + &value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); } + + Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i], + col_indices[0])); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } - else - { - int trilinos_i = i; - int trilinos_j = j; - const TrilinosScalar* value_ptr = &value; + } - // If the matrix is not yet filled, - // we can insert new entries into - // the matrix using this - // command. Otherwise, we're just - // able to replace existing entries. - if (matrix->Filled() == false) - { - ierr = matrix->InsertGlobalValues (1, &trilinos_i, - 1, &trilinos_j, - &value_ptr, - Epetra_FECrsMatrix::ROW_MAJOR); - // When adding up elements, we do - // not want to create exceptions in - // the case when adding elements. - if (ierr > 0) - ierr = 0; - } - else + void + SparseMatrix::add (const unsigned int i, + const unsigned int j, + const TrilinosScalar value) + { + + Assert (numbers::is_finite(value), + ExcMessage("The given value is not finite but either " + "infinite or Not A Number (NaN)")); + + if (value == 0) + { + // we have to do above actions in any case + // to be consistent with the MPI + // communication model (see the comments + // in the documentation of + // TrilinosWrappers::Vector), but we can + // save some work if the addend is + // zero. However, these actions are done + // in case we pass on to the other + // function. + if (last_action == Insert) { - ierr = matrix->ReplaceGlobalValues (1, &trilinos_i, - 1, &trilinos_j, - &value_ptr, - Epetra_FECrsMatrix::ROW_MAJOR); + int ierr; + ierr = matrix->GlobalAssemble(col_map, row_map, false); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } + + if (last_action != Add) + last_action = Add; + + return; } + else + add (1, &i, 1, &j, &value); + } - Assert (ierr <= 0, ExcAccessToNonPresentElement(i,j)); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + void + SparseMatrix::add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &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())); + + add (row_indices.size(), &row_indices[0], + col_indices.size(), &col_indices[0], &values(0,0)); } void - SparseMatrix::add (const unsigned int i, - const unsigned int j, - const TrilinosScalar value) + SparseMatrix::add (const unsigned int row, + const std::vector &col_indices, + const std::vector &values) { + Assert (col_indices.size() == values.size(), + ExcDimensionMismatch(col_indices.size(), values.size())); + + add (1, &row, col_indices.size(), &col_indices[0], &values[0]); + } + - Assert (numbers::is_finite(value), - ExcMessage("The given value is not finite but either " - "infinite or Not A Number (NaN)")); + void + SparseMatrix::add (const unsigned int n_rows, + const unsigned int *row_indices, + const unsigned int n_cols, + const unsigned int *col_indices, + const TrilinosScalar *values) + { + int ierr; if (last_action == Insert) { - int ierr; ierr = matrix->GlobalAssemble(col_map, row_map, false); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - } + } if (last_action != Add) last_action = Add; - // we have to do above actions in any - // case to be consistent with the MPI - // communication model (see the - // comments in the documentation of - // TrilinosWrappers::Vector), but we - // can save some work if the addend is - // zero - if (value == 0) - return; - - int ierr; + // Now go through all rows that are + // present in the input data. + for (unsigned int i=0; iEpetra_CrsMatrix::SumIntoGlobalValues( - (int)i, 1, - const_cast(&value), - (int*)&j); - } - else - { + if (row_map.MyGID(i) == true) + { + ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues( + (int)row_indices[i], (int)n_cols, + const_cast(&values[i*n_cols]), + (int*)&col_indices[0]); + } + else + { // When we're at off-processor data, we // have to stick with the standard // SumIntoGlobalValues @@ -820,21 +906,20 @@ namespace TrilinosWrappers // call the function we already use, // which is very unefficient if writing // one element at a time). - compressed = false; + compressed = false; - int trilinos_i = i; - int trilinos_j = j; + const TrilinosScalar* value_ptr = &values[i*n_cols]; - const TrilinosScalar* value_ptr = &value; + ierr = matrix->SumIntoGlobalValues (1, (int*)&i, + (int)n_cols, (int*)&col_indices[0], + &value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); + } - ierr = matrix->SumIntoGlobalValues (1, &trilinos_i, - 1, &trilinos_j, - &value_ptr, - Epetra_FECrsMatrix::ROW_MAJOR); + Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i], + col_indices[0])); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } - - AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -1243,6 +1328,11 @@ namespace TrilinosWrappers { Assert (&src != &dst, ExcSourceEqualsDestination()); + // Choose to reinit the vector with fast + // argument set, which does not overwrite + // the content -- this is what we need + // since we're going to overwrite that + // anyway in the vmult operation. temp_vector.reinit(dst, true); vmult (temp_vector, src);