From a989d56d7be4ba1ba398589356fd7c2388efb492 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 24 Nov 2008 16:25:46 +0000 Subject: [PATCH] The function TrilinosWrappers::SparseMatrix::add() that sets more than one element can now elide zero elements in the array. git-svn-id: https://svn.dealii.org/trunk@17697 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/trilinos_block_sparse_matrix.h | 202 ++++++---- .../lac/include/lac/trilinos_sparse_matrix.h | 352 ++++++++++++------ deal.II/lac/source/trilinos_sparse_matrix.cc | 22 +- 3 files changed, 376 insertions(+), 200 deletions(-) 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 2b3fb37901..92db2a99aa 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -305,10 +305,20 @@ namespace TrilinosWrappers * 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 FullMatrix &full_matrix, + const bool elide_zero_values = false); /** * Set several elements in the @@ -316,25 +326,42 @@ namespace TrilinosWrappers * 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, + void set (const unsigned int row, const std::vector &col_indices, - const std::vector &values); + const std::vector &values, + const bool elide_zero_values = false); /** * 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. + * 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 n_rows, - const unsigned int *row_indices, + void set (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values); + const TrilinosScalar *values, + const bool elide_zero_values = false); /** * Add @p value to the element @@ -350,10 +377,21 @@ namespace TrilinosWrappers * 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 FullMatrix &full_matrix, + const bool elide_zero_values = true); /** * Set several elements in the @@ -361,25 +399,44 @@ namespace TrilinosWrappers * 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, + void add (const unsigned int row, const std::vector &col_indices, - const std::vector &values); + const std::vector &values, + const bool elide_zero_values = true); /** * 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. + * 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 n_rows, - const unsigned int *row_indices, + void add (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values); + const TrilinosScalar *values, + const bool elide_zero_values = true); /** * Matrix-vector multiplication: @@ -720,7 +777,7 @@ namespace TrilinosWrappers * writing data from local to global * matrices. */ - std::vector block_col_indices, local_row_length, local_col_indices; + std::vector block_col_indices, local_col_indices, local_row_length; }; @@ -762,15 +819,17 @@ namespace TrilinosWrappers void BlockSparseMatrix::set (const std::vector &row_indices, const std::vector &col_indices, - const FullMatrix &values) + 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())); - set (row_indices.size(), &row_indices[0], - col_indices.size(), &col_indices[0], &values(0,0)); + for (unsigned int i=0; i &col_indices, - const std::vector &values) + const std::vector &values, + const bool elide_zero_values) { Assert (col_indices.size() == values.size(), ExcDimensionMismatch(col_indices.size(), values.size())); - set (1, &row, col_indices.size(), &col_indices[0], &values[0]); + 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 n_rows, - const unsigned int *row_indices, + BlockSparseMatrix::set (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values) + 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_row_length.resize (this->n_block_cols()); local_col_indices.resize (n_cols); // Clear the content in local_row_length @@ -809,10 +873,7 @@ namespace TrilinosWrappers // 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. + // matrix. { unsigned int current_block = 0, row_length = 0; block_col_indices[0] = 0; @@ -853,21 +914,19 @@ namespace TrilinosWrappers // 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); + for (unsigned int block_col=0; block_col - 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) + 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())); - add (row_indices.size(), &row_indices[0], - col_indices.size(), &col_indices[0], &values(0,0)); + for (unsigned int i=0; i &col_indices, - const std::vector &values) + const std::vector &values, + const bool elide_zero_values) { Assert (col_indices.size() == values.size(), ExcDimensionMismatch(col_indices.size(), values.size())); - add (1, &row, col_indices.size(), &col_indices[0], &values[0]); + add (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); } inline void - BlockSparseMatrix::add (const unsigned int n_rows, - const unsigned int *row_indices, + BlockSparseMatrix::add (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values) + const TrilinosScalar *values, + const bool elide_zero_values) { // TODO: Look over this to find out // whether we can do that more @@ -943,10 +1006,7 @@ namespace TrilinosWrappers // 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. + // matrix. { unsigned int current_block = 0, row_length = 0; block_col_indices[0] = 0; @@ -987,21 +1047,19 @@ namespace TrilinosWrappers // 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); + for (unsigned int block_col=0; block_col - row_index = this->row_block_indices.global_to_local (row_indices[i]); - for (unsigned int block_col=0; block_colelide_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 FullMatrix &full_matrix, + const bool elide_zero_values = false); /** * Set several elements in the @@ -905,19 +915,27 @@ namespace TrilinosWrappers * insertion of elements at positions * which have not been initialized * will throw an exception. + * + * 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 std::vector &values, + const bool elide_zero_values = false); /** * 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. + * given by values in a + * given row in columns given by + * col_indices into the sparse + * matrix. * * This function is able to insert * new elements into the matrix as @@ -929,12 +947,21 @@ namespace TrilinosWrappers * insertion of elements at positions * which have not been initialized * will throw an exception. + * + * 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 n_rows, - const unsigned int *row_indices, + void set (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values); + const TrilinosScalar *values, + const bool elide_zero_values = false); /** * Add @p value to the element @@ -970,10 +997,21 @@ namespace TrilinosWrappers * throws an exception if an * entry does not exist in the * sparsity pattern. + * + * 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 FullMatrix &full_matrix, + const bool elide_zero_values = true); /** * Set several elements in the @@ -990,34 +1028,52 @@ namespace TrilinosWrappers * throws an exception if an * entry does not exist in the * sparsity pattern. + * + * 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 std::vector &values, + const bool elide_zero_values = true); /** - * 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. + * Add an array of values given by + * values in the given + * global matrix row at columns + * specified by col_indices in the + * sparse matrix. * * 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 + * 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. + * + * 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 n_rows, - const unsigned int *row_indices, + void add (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values); + const TrilinosScalar *values, + const bool elide_zero_values = false); /** * Multiply the entire matrix @@ -1683,6 +1739,24 @@ namespace TrilinosWrappers */ mutable VectorBase temp_vector; + /** + * An internal array of integer + * values that is used to store the + * column indices when + * adding/inserting local data into + * the (large) sparse matrix. + */ + std::vector column_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; + public: /** * A sparse matrix object in @@ -1954,7 +2028,7 @@ namespace TrilinosWrappers ExcMessage("The given value is not finite but either " "infinite or Not A Number (NaN)")); - set (1, &i, 1, &j, &value); + set (i, 1, &j, &value, false); } @@ -1963,15 +2037,17 @@ namespace TrilinosWrappers void SparseMatrix::set (const std::vector &row_indices, const std::vector &col_indices, - const FullMatrix &values) + 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())); - set (row_indices.size(), &row_indices[0], - col_indices.size(), &col_indices[0], &values(0,0)); + for (unsigned int i=0; i &col_indices, - const std::vector &values) + const std::vector &values, + const bool elide_zero_values) { Assert (col_indices.size() == values.size(), ExcDimensionMismatch(col_indices.size(), values.size())); - set (1, &row, col_indices.size(), &col_indices[0], &values[0]); + set (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); } inline void - SparseMatrix::set (const unsigned int n_rows, - const unsigned int *row_indices, + SparseMatrix::set (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values) + const TrilinosScalar *values, + const bool elide_zero_values) { int ierr; if (last_action == Add) @@ -2009,11 +2087,38 @@ namespace TrilinosWrappers if (last_action != Insert) last_action = Insert; - // Now go through all rows that are - // present in the input data. - for (unsigned int i=0; iFilled() == false) { - if (matrix->Filled() == false) - { - ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues( - row, (int)n_cols, - const_cast(&values[i*n_cols]), - (int*)&col_indices[0]); + ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_cols, + const_cast(col_value_ptr), + col_index_ptr); // When inserting elements, we do // not want to create exceptions in // the case when inserting non-local // data (since that's what we want // to do right now). - if (ierr > 0) - ierr = 0; - } - else - ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues( - row, (int)n_cols, - const_cast(&values[i*n_cols]), - (int*)&col_indices[0]); + if (ierr > 0) + ierr = 0; } else - { + ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_cols, + const_cast(col_value_ptr), + col_index_ptr); + } + else + { // When we're at off-processor data, we // have to stick with the standard // Insert/ReplaceGlobalValues @@ -2064,30 +2167,26 @@ namespace TrilinosWrappers // call the function we already use, // which is very unefficient if writing // one element at a time). - compressed = false; + compressed = false; - const TrilinosScalar* value_ptr = &values[i*n_cols]; - - if (matrix->Filled() == false) - { - ierr = matrix->InsertGlobalValues (1, &row, - (int)n_cols, (int*)&col_indices[0], - &value_ptr, - Epetra_FECrsMatrix::ROW_MAJOR); - if (ierr > 0) - ierr = 0; - } - else - ierr = matrix->ReplaceGlobalValues (1, &row, - (int)n_cols, (int*)&col_indices[0], - &value_ptr, - Epetra_FECrsMatrix::ROW_MAJOR); + if (matrix->Filled() == false) + { + ierr = matrix->InsertGlobalValues (1, (int*)&row, + (int)n_cols, col_index_ptr, + &col_value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); + if (ierr > 0) + ierr = 0; } - - Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i], - col_indices[0])); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + else + ierr = matrix->ReplaceGlobalValues (1, (int*)&row, + (int)n_cols, col_index_ptr, + &col_value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); } + + Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0])); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -2105,7 +2204,8 @@ namespace TrilinosWrappers if (value == 0) { - // we have to do above actions in any case + // we have to do checkings on Insert/Add + // in any case // to be consistent with the MPI // communication model (see the comments // in the documentation of @@ -2128,7 +2228,7 @@ namespace TrilinosWrappers return; } else - add (1, &i, 1, &j, &value); + add (i, 1, &j, &value, false); } @@ -2137,15 +2237,17 @@ namespace TrilinosWrappers void SparseMatrix::add (const std::vector &row_indices, const std::vector &col_indices, - const FullMatrix &values) + 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())); - add (row_indices.size(), &row_indices[0], - col_indices.size(), &col_indices[0], &values(0,0)); + for (unsigned int i=0; i &col_indices, - const std::vector &values) + const std::vector &values, + const bool elide_zero_values) { Assert (col_indices.size() == values.size(), ExcDimensionMismatch(col_indices.size(), values.size())); - add (1, &row, col_indices.size(), &col_indices[0], &values[0]); + add (row, col_indices.size(), &col_indices[0], &values[0], + elide_zero_values); } inline void - SparseMatrix::add (const unsigned int n_rows, - const unsigned int *row_indices, + SparseMatrix::add (const unsigned int row, const unsigned int n_cols, const unsigned int *col_indices, - const TrilinosScalar *values) + const TrilinosScalar *values, + const bool elide_zero_values) { int ierr; if (last_action == Insert) @@ -2183,25 +2287,51 @@ namespace TrilinosWrappers if (last_action != Add) last_action = Add; - // Now go through all rows that are - // present in the input data. - for (unsigned int i=0; iEpetra_CrsMatrix::SumIntoGlobalValues( - row, (int)n_cols, - const_cast(&values[i*n_cols]), - (int*)&col_indices[0]); - } - else - { + if (row_map.MyGID(row) == true) + { + ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_cols, + const_cast(col_value_ptr), + col_index_ptr); + } + else + { // When we're at off-processor data, we // have to stick with the standard // SumIntoGlobalValues @@ -2212,20 +2342,16 @@ namespace TrilinosWrappers // 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]; - - ierr = matrix->SumIntoGlobalValues (1, &row, - (int)n_cols, (int*)&col_indices[0], - &value_ptr, - Epetra_FECrsMatrix::ROW_MAJOR); - } + compressed = false; - Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i], - col_indices[0])); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + ierr = matrix->SumIntoGlobalValues (1, (int*)&row, (int)n_cols, + col_index_ptr, + &col_value_ptr, + Epetra_FECrsMatrix::ROW_MAJOR); } + + Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0])); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 18570f3f9d..d45ebca374 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -316,7 +316,7 @@ namespace TrilinosWrappers // sparsity_pattern.n_cols())); std::vector values; - std::vector row_indices; + std::vector row_indices; for (unsigned int row=0; rowEpetra_CrsMatrix::InsertGlobalValues ((int)row, row_length, - &values[0], - &row_indices[0]); + + set (row, row_length, &row_indices[0], &values[0], false); } last_action = Zero; @@ -416,7 +414,7 @@ namespace TrilinosWrappers // sparsity_pattern.n_cols())); std::vector values; - std::vector row_indices; + std::vector row_indices; for (unsigned int row=0; rowEpetra_CrsMatrix::InsertGlobalValues ((int)row, row_length, - &values[0], - &row_indices[0]); + set (row, row_length, &row_indices[0], &values[0], false); } last_action = Zero; @@ -536,7 +532,7 @@ namespace TrilinosWrappers &n_entries_per_row[0], false)); std::vector values; - std::vector row_indices; + std::vector row_indices; for (unsigned int row=0; rowEpetra_CrsMatrix::InsertGlobalValues ((int)row, row_length, - &values[0], - &row_indices[0]); + set (row, index, &row_indices[0], &values[0], false); } compress(); -- 2.39.5