From: kronbichler Date: Wed, 26 Nov 2008 15:39:21 +0000 (+0000) Subject: Use row-wise addition into matrices also for PETSc wrappers matrix. When building... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a2714d04154da08df3fa54ca9ead56e0e2cff932;p=dealii-svn.git Use row-wise addition into matrices also for PETSc wrappers matrix. When building a parallel PETSc matrix from a sparsity pattern, set the option NO_NEW_ENTRIES. git-svn-id: https://svn.dealii.org/trunk@17751 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/petsc_matrix_base.h b/deal.II/lac/include/lac/petsc_matrix_base.h index fe9b86d8e2..7436e1b037 100644 --- a/deal.II/lac/include/lac/petsc_matrix_base.h +++ b/deal.II/lac/include/lac/petsc_matrix_base.h @@ -16,6 +16,7 @@ #include #include +#include #include #ifdef DEAL_II_USE_PETSC @@ -330,49 +331,275 @@ namespace PETScWrappers void clear (); /** - * Set the element (i,j) - * to @p value. + * Set the element (i,j) to @p + * value. * - * If the present object (from - * a derived class of this one) - * happens to be a sparse - * matrix, then this function - * adds a new entry to the - * matrix if it didn't exist - * before, very much in - * contrast to the SparseMatrix - * class which throws an error - * if the entry does not exist. - * If value is not a - * finite number an exception - * is thrown. + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds a new entry to the + * matrix if it didn't exist before, + * very much in contrast to the + * SparseMatrix class which throws an + * error if the entry does not exist. + * If value is not a finite + * number an exception is thrown. */ void set (const unsigned int i, const unsigned int j, const PetscScalar value); /** - * Add @p value to the - * element (i,j). + * 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. * - * If the present object (from - * a derived class of this one) - * happens to be a sparse - * matrix, then this function - * adds a new entry to the - * matrix if it didn't exist - * before, very much in - * contrast to the SparseMatrix - * class which throws an error - * if the entry does not exist. - * If value is not a - * finite number an exception - * is thrown. + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds some new entries to + * the matrix if they didn't exist + * before, very much in contrast to + * the SparseMatrix class which + * throws an error if the entry does + * not exist. + * + * 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 &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. + * + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds some new entries to + * the matrix if they didn't exist + * before, very much in contrast to + * the SparseMatrix class which + * throws an error if the entry does + * not exist. + * + * 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 a + * given row in columns given by + * col_indices into the sparse + * matrix. + * + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds some new entries to + * the matrix if they didn't exist + * before, very much in contrast to + * the SparseMatrix class which + * throws an error if the entry does + * not exist. + * + * 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 PetscScalar *values, + const bool elide_zero_values = false); + + /** + * Add @p value to the element + * (i,j). + * + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds a new entry to the + * matrix if it didn't exist before, + * very much in contrast to the + * SparseMatrix class which throws an + * error if the entry does not exist. + * If value is not a finite + * number an exception is thrown. */ void add (const unsigned int i, const unsigned int j, const PetscScalar 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. + * + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds some new entries to + * the matrix if they didn't exist + * before, very much in contrast to + * the SparseMatrix class which + * throws an error if the entry does + * not exist. + * + * 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. + * + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds some new entries to + * the matrix if they didn't exist + * before, very much in contrast to + * the SparseMatrix class which + * throws an error if the entry does + * not exist. + * + * 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. + * + * If the present object (from a + * derived class of this one) happens + * to be a sparse matrix, then this + * function adds some new entries to + * the matrix if they didn't exist + * before, very much in contrast to + * the SparseMatrix class which + * throws an error if the entry does + * not exist. + * + * 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 PetscScalar *values, + const bool elide_zero_values = true); + /** * Remove all elements from * this row by setting @@ -887,7 +1114,7 @@ namespace PETScWrappers * Exception */ DeclException0 (ExcSourceEqualsDestination); - + protected: /** * A generic matrix object in @@ -926,7 +1153,27 @@ namespace PETScWrappers * Store whether the last action was a * write or add operation. */ - LastAction::Values last_action; + LastAction::Values last_action; + + private: + /** + * 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; + }; @@ -1082,8 +1329,322 @@ namespace PETScWrappers } } - - + + + + // 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. + inline + void + MatrixBase::set (const unsigned int i, + const unsigned int j, + const PetscScalar 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); + } + + + + inline + void + MatrixBase::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 &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); + } + + + + inline + void + MatrixBase::set (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const PetscScalar *values, + const bool elide_zero_values) + { + + if (last_action != LastAction::insert) + { + int ierr; + ierr = MatAssemblyBegin(matrix,MAT_FLUSH_ASSEMBLY); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = MatAssemblyEnd(matrix,MAT_FLUSH_ASSEMBLY); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + last_action = LastAction::insert; + } + + const signed int petsc_i = row; + int * col_index_ptr; + PetscScalar const* col_value_ptr; + int n_columns; + + // If we don't elide zeros, the pointers + // are already available... + if (elide_zero_values == false) + { + col_index_ptr = (int*)col_indices; + col_value_ptr = values; + n_columns = n_cols; + } + else + { + // Otherwise, extract nonzero values in + // each row and pass on to the other + // function. + column_indices.resize(n_cols); + column_values.resize(n_cols); + + n_columns = 0; + for (unsigned int j=0; j &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 &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 + MatrixBase::add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const PetscScalar *values, + const bool elide_zero_values) + { + if (last_action != LastAction::add) + { + int ierr; + ierr = MatAssemblyBegin(matrix,MAT_FLUSH_ASSEMBLY); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = MatAssemblyEnd(matrix,MAT_FLUSH_ASSEMBLY); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + last_action = LastAction::add; + } + + + const signed int petsc_i = row; + int * col_index_ptr; + PetscScalar const* col_value_ptr; + int n_columns; + + // If we don't elide zeros, the pointers + // are already available... + if (elide_zero_values == false) + { + col_index_ptr = (int*)col_indices; + col_value_ptr = values; + n_columns = n_cols; + } + else + { + // Otherwise, extract nonzero values in + // each row and pass on to the other + // function. + column_indices.resize(n_cols); + column_values.resize(n_cols); + + n_columns = 0; + for (unsigned int j=0; j into the sparse - * matrix, according to row_indices - * and col_indices. + * 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. * * This function is able to insert * new elements into the matrix as @@ -893,6 +903,17 @@ namespace TrilinosWrappers * false, i.e., even zero * values are inserted/replaced. */ + 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, @@ -983,11 +1004,21 @@ namespace TrilinosWrappers 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. + * 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. * * Just as the respective call in * deal.II SparseMatrix @@ -1008,6 +1039,17 @@ namespace TrilinosWrappers * 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, @@ -2017,6 +2059,13 @@ namespace TrilinosWrappers + // 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. inline void SparseMatrix::set (const unsigned int i, @@ -2033,6 +2082,23 @@ namespace TrilinosWrappers + 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 &row_indices, @@ -2109,12 +2175,19 @@ namespace TrilinosWrappers n_columns = 0; for (unsigned int j=0; j &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 &row_indices, @@ -2311,12 +2401,19 @@ namespace TrilinosWrappers n_columns = 0; for (unsigned int j=0; j