]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use row-wise addition into matrices also for PETSc wrappers matrix. When building...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 26 Nov 2008 15:39:21 +0000 (15:39 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 26 Nov 2008 15:39:21 +0000 (15:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@17751 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/source/petsc_matrix_base.cc
deal.II/lac/source/petsc_parallel_sparse_matrix.cc
deal.II/lac/source/petsc_sparse_matrix.cc

index fe9b86d8e2c1f75a71391f66f7fb00cfffffee1e..7436e1b037bbdfd336e7cad0bdc9bc753f94b583 100644 (file)
@@ -16,6 +16,7 @@
 
 #include <base/config.h>
 #include <base/subscriptor.h>
+#include <lac/full_matrix.h>
 #include <lac/exceptions.h>
 
 #ifdef DEAL_II_USE_PETSC
@@ -330,49 +331,275 @@ namespace PETScWrappers
       void clear ();
 
                                        /**
-                                        * Set the element (<i>i,j</i>)
-                                        * to @p value.
+                                        * Set the element (<i>i,j</i>) 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 <tt>value</tt> 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 <tt>value</tt> 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>i,j</i>).
+                                        * Set all elements given in a
+                                        * FullMatrix<double> into the sparse
+                                        * matrix locations given by
+                                        * <tt>indices</tt>. In other words,
+                                        * this function writes the elements
+                                        * in <tt>full_matrix</tt> into the
+                                        * calling matrix, using the
+                                        * local-to-global indexing specified
+                                        * by <tt>indices</tt> 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 <tt>value</tt> 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
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
+                                       */
+      void set (const std::vector<unsigned int> &indices,
+               const FullMatrix<PetscScalar>   &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<unsigned int> &row_indices,
+               const std::vector<unsigned int> &col_indices,
+               const FullMatrix<PetscScalar>   &full_matrix,
+               const bool                       elide_zero_values = false);
+
+                                       /**
+                                        * Set several elements in the
+                                        * specified row of the matrix with
+                                        * column indices as given by
+                                        * <tt>col_indices</tt> 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
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
+                                       */
+      void set (const unsigned int               row,
+               const std::vector<unsigned int> &col_indices,
+               const std::vector<PetscScalar>  &values,
+               const bool                       elide_zero_values = false);
+
+                                       /**
+                                        * Set several elements to values
+                                        * given by <tt>values</tt> 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
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, 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>i,j</i>).
+                                       *
+                                       * 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 <tt>value</tt> 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<double> into sparse
+                                        * matrix locations given by
+                                        * <tt>indices</tt>. In other words,
+                                        * this function adds the elements in
+                                        * <tt>full_matrix</tt> to the
+                                        * respective entries in calling
+                                        * matrix, using the local-to-global
+                                        * indexing specified by
+                                        * <tt>indices</tt> 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
+                                       * <tt>elide_zero_values</tt> 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 <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
+                                       */
+      void add (const std::vector<unsigned int> &indices,
+               const FullMatrix<PetscScalar>   &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<unsigned int> &row_indices,
+               const std::vector<unsigned int> &col_indices,
+               const FullMatrix<PetscScalar>   &full_matrix,
+               const bool                       elide_zero_values = true);
+
+                                       /**
+                                        * Set several elements in the
+                                        * specified row of the matrix with
+                                        * column indices as given by
+                                        * <tt>col_indices</tt> 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
+                                       * <tt>elide_zero_values</tt> 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 <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
+                                       */
+      void add (const unsigned int               row,
+               const std::vector<unsigned int> &col_indices,
+               const std::vector<PetscScalar>  &values,
+               const bool                       elide_zero_values = true);
+
+                                       /**
+                                        * Add an array of values given by
+                                        * <tt>values</tt> 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
+                                       * <tt>elide_zero_values</tt> 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 <tt>true</tt>,
+                                       * 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 <tt>row</tt> 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<unsigned int> 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<PetscScalar> 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<unsigned int> &indices,
+                  const FullMatrix<PetscScalar>   &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<indices.size(); ++i)
+      set (indices[i], indices.size(), &indices[0], &values(i,0),
+          elide_zero_values);
+  }
+
+
+
+  inline
+  void
+  MatrixBase::set (const std::vector<unsigned int> &row_indices,
+                  const std::vector<unsigned int> &col_indices,
+                  const FullMatrix<PetscScalar>   &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<row_indices.size(); ++i)
+      set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
+          elide_zero_values);
+  }
+
+
+
+  inline
+  void
+  MatrixBase::set (const unsigned int               row,
+                  const std::vector<unsigned int> &col_indices,
+                  const std::vector<PetscScalar>  &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<n_cols; ++j)
+         {
+           const double value = values[j];
+           Assert (numbers::is_finite(value),
+                   ExcMessage("The given value is not finite but either "
+                              "infinite or Not A Number (NaN)"));
+           if (value != 0)
+             {
+               column_indices[n_columns] = col_indices[j];
+               column_values[n_columns] = value;
+               n_columns++;
+             }
+         }
+       Assert(n_columns <= (int)n_cols, ExcInternalError());
+
+       col_index_ptr = (int*)&column_indices[0];
+       col_value_ptr = &column_values[0];
+      }
+          
+    const int ierr
+      = MatSetValues (matrix, 1, &petsc_i, n_columns, col_index_ptr,
+                     col_value_ptr, INSERT_VALUES);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+  }
+
+
+
+  inline
+  void
+  MatrixBase::add (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)"));
+
+    if (value == 0)
+      {
+                                 // 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
+                                 // 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 != 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;
+         }
+
+       return;
+      }
+    else
+      add (i, 1, &j, &value, false);
+  }
+
+
+
+  inline
+  void
+  MatrixBase::add (const std::vector<unsigned int> &indices,
+                  const FullMatrix<PetscScalar>   &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<indices.size(); ++i)
+      add (indices[i], indices.size(), &indices[0], &values(i,0),
+          elide_zero_values);
+  }
+
+
+
+  inline
+  void
+  MatrixBase::add (const std::vector<unsigned int> &row_indices,
+                  const std::vector<unsigned int> &col_indices,
+                  const FullMatrix<PetscScalar>   &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<row_indices.size(); ++i)
+      add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
+          elide_zero_values);
+  }
+
+
+
+  inline
+  void
+  MatrixBase::add (const unsigned int               row,
+                  const std::vector<unsigned int> &col_indices,
+                  const std::vector<PetscScalar>  &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<n_cols; ++j)
+         {
+           const double value = values[j];
+           Assert (numbers::is_finite(value),
+                   ExcMessage("The given value is not finite but either "
+                              "infinite or Not A Number (NaN)"));
+           if (value != 0)
+             {
+               column_indices[n_columns] = col_indices[j];
+               column_values[n_columns] = value;
+               n_columns++;
+             }
+         }
+       Assert(n_columns <= (int)n_cols, ExcInternalError());
+
+       col_index_ptr = (int*)&column_indices[0];
+       col_value_ptr = &column_values[0];
+      }
+          
+    const int ierr
+      = MatSetValues (matrix, 1, &petsc_i, n_columns, col_index_ptr,
+                     col_value_ptr, ADD_VALUES);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+  }
+
+
+
+
+
+
   inline
   PetscScalar
   MatrixBase::operator() (const unsigned int i,
index dd46556c999fa07fb0d5a02a4d91f489f5e24a7d..5441490d80d8c4defdf50a09cf4b54e7172badc3 100755 (executable)
@@ -831,14 +831,14 @@ namespace TrilinosWrappers
                                        * transpose.  TODO: Not
                                        * implemented.
                                        */
-
       bool is_hermitian () const;
-                                    /**
-                                     * Determine an estimate for the
-                                     * memory consumption (in bytes)
-                                     * of this object. Currently not
-                                     * implemented for this class.
-                                     */
+
+                                      /**
+                                       * Determine an estimate for the
+                                       * memory consumption (in bytes)
+                                       * of this object. Currently not
+                                       * implemented for this class.
+                                       */
       unsigned int memory_consumption () const;
 
 //@}
@@ -870,8 +870,18 @@ namespace TrilinosWrappers
                                        /**
                                         * Set all elements given in a
                                         * FullMatrix<double> into the sparse
-                                        * matrix, according to row_indices
-                                        * and col_indices.
+                                        * matrix locations given by
+                                        * <tt>indices</tt>. In other words,
+                                        * this function writes the elements
+                                        * in <tt>full_matrix</tt> into the
+                                        * calling matrix, using the
+                                        * local-to-global indexing specified
+                                        * by <tt>indices</tt> 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
                                        * <tt>false</tt>, i.e., even zero
                                        * values are inserted/replaced.
                                        */
+      void set (const std::vector<unsigned int>  &indices,
+               const FullMatrix<TrilinosScalar> &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<unsigned int>  &row_indices,
                const std::vector<unsigned int>  &col_indices,
                const FullMatrix<TrilinosScalar> &full_matrix,
@@ -983,11 +1004,21 @@ namespace TrilinosWrappers
                 const TrilinosScalar value);
 
                                        /**
-                                        * Add all elements in a FullMatrix
-                                        * <tt>full_matrix</tt> (cell matrix)
-                                        * into the sparse matrix, at
-                                        * locations according to row_indices
-                                        * and col_indices.
+                                        * Add all elements given in a
+                                        * FullMatrix<double> into sparse
+                                        * matrix locations given by
+                                        * <tt>indices</tt>. In other words,
+                                        * this function adds the elements in
+                                        * <tt>full_matrix</tt> to the
+                                        * respective entries in calling
+                                        * matrix, using the local-to-global
+                                        * indexing specified by
+                                        * <tt>indices</tt> 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<Number>
@@ -1008,6 +1039,17 @@ namespace TrilinosWrappers
                                        * i.e., zero values won't be added
                                        * into the matrix.
                                        */
+      void add (const std::vector<unsigned int>  &indices,
+               const FullMatrix<TrilinosScalar> &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<unsigned int>  &row_indices,
                const std::vector<unsigned int>  &col_indices,
                const FullMatrix<TrilinosScalar> &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<unsigned int>  &indices,
+                    const FullMatrix<TrilinosScalar> &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<indices.size(); ++i)
+      set (indices[i], indices.size(), &indices[0], &values(i,0),
+          elide_zero_values);
+  }
+
+
+
   inline
   void
   SparseMatrix::set (const std::vector<unsigned int>  &row_indices,
@@ -2109,12 +2175,19 @@ namespace TrilinosWrappers
 
        n_columns = 0;
        for (unsigned int j=0; j<n_cols; ++j)
-         if (values[j] != 0)
-           {
-             column_indices[n_columns] = col_indices[j];
-             column_values[n_columns] = values[j];
-             n_columns++;
-           }
+         {
+           const double value = values[j];
+           Assert (numbers::is_finite(value),
+                   ExcMessage("The given value is not finite but either "
+                              "infinite or Not A Number (NaN)"));
+           if (value != 0)
+             {
+               column_indices[n_columns] = col_indices[j];
+               column_values[n_columns] = value;
+               n_columns++;
+             }
+         }
+
        Assert(n_columns <= (int)n_cols, ExcInternalError());
 
        col_index_ptr = (int*)&column_indices[0];
@@ -2235,6 +2308,23 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  void
+  SparseMatrix::add (const std::vector<unsigned int>  &indices,
+                    const FullMatrix<TrilinosScalar> &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<indices.size(); ++i)
+      add (indices[i], indices.size(), &indices[0], &values(i,0),
+          elide_zero_values);
+  }
+
+
+
   inline
   void
   SparseMatrix::add (const std::vector<unsigned int>  &row_indices,
@@ -2311,12 +2401,19 @@ namespace TrilinosWrappers
 
        n_columns = 0;
        for (unsigned int j=0; j<n_cols; ++j)
-         if (values[j] != 0)
-           {
-             column_indices[n_columns] = col_indices[j];
-             column_values[n_columns] = values[j];
-             n_columns++;
-           }
+         {
+           const double value = values[j];
+           Assert (numbers::is_finite(value),
+                   ExcMessage("The given value is not finite but either "
+                              "infinite or Not A Number (NaN)"));
+           if (value != 0)
+             {
+               column_indices[n_columns] = col_indices[j];
+               column_values[n_columns] = value;
+               n_columns++;
+             }
+         }
+
        Assert(n_columns <= (int)n_cols, ExcInternalError());
 
        col_index_ptr = (int*)&column_indices[0];
index 6c0f3a51c056ae6d93cec3ea2ae20ef42128807b..91c046ec9d37a74b528e570f967803fe47dd3286 100644 (file)
@@ -129,80 +129,6 @@ namespace PETScWrappers
 
     return *this;
   }
-    
-
-
-  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)"));
-
-    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));
-      }
-    
-    const signed int petsc_i = i;
-    const signed int petsc_j = j;
-          
-    const int ierr
-      = MatSetValues (matrix, 1, &petsc_i, 1, &petsc_j,
-                      &value, INSERT_VALUES);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-    last_action = LastAction::insert;
-  }
-
-
-
-  void
-  MatrixBase::add (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)"));
-    
-    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;
-      }
-    
-                                     // we have to do above actions in any
-                                     // case to be consistent with the MPI
-                                     // communication model (see the
-                                     // comments in the documentation of
-                                     // PETScWrappers::MPI::Vector), but we
-                                     // can save some work if the addend is
-                                     // zero
-    if (value == 0)
-      return;
-
-    const signed int petsc_i = i;
-    const signed int petsc_j = j;
-          
-    const int ierr
-      = MatSetValues (matrix, 1, &petsc_i, 1, &petsc_j,
-                      &value, ADD_VALUES);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));
-  }
 
 
 
index 44c65a3d01ed85187e9bcbb6925b65db0ef35414..04c0aba9a19fb3b98b2a4d7de59888cab46418f8 100644 (file)
@@ -449,7 +449,11 @@ namespace PETScWrappers
           
 #endif
 
-//TODO: We should use the appropriate option with MatSetOption here indicating that we do not intend to add any further matrix entries. Maybe this would accelerate some things as well
+                                          // Now we won't insert any
+                                          // further entries, so PETSc can
+                                          // internally optimize some data
+                                          // structures.
+         MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
         }
     }
 
index 95ca455d4d5c07524cda683a61695527dd3f30da..9c2b655bb521f08fd4a75bb0d740f7ce6e222cc7 100644 (file)
@@ -236,6 +236,11 @@ namespace PETScWrappers
           }
         compress ();
       }
+
+                                       // In the end, tell the matrix that
+                                       // it should not expect any new
+                                       // entries.
+    MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.