]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Introduced a few add() and set() functions for Trilinos matrices that allow to set...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 21 Nov 2008 18:33:29 +0000 (18:33 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 21 Nov 2008 18:33:29 +0000 (18:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@17678 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/source/trilinos_block_sparse_matrix.cc
deal.II/lac/source/trilinos_sparse_matrix.cc

index 90133e1bdc68174051d80304cfa6bbdaa2e77fc9..0032a433dcd64d2c650b05add47fa664ac535890 100644 (file)
@@ -20,6 +20,7 @@
 #include <lac/block_sparse_matrix.h>
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_block_vector.h>
+#include <lac/full_matrix.h>
 #include <lac/exceptions.h>
 
 #include <cmath>
@@ -291,6 +292,95 @@ namespace TrilinosWrappers
                                         */
       unsigned int n_nonzero_elements () const;
 
+                                       /**
+                                        * Set the element (<i>i,j</i>)
+                                        * to @p value.
+                                       */
+      void set (const unsigned int i,
+                const unsigned int j,
+                const TrilinosScalar value);
+
+                                       /**
+                                        * Set all elements given in a
+                                        * FullMatrix<double> into the sparse
+                                        * matrix, according to row_indices
+                                        * and col_indices.
+                                       */
+      void set (const std::vector<unsigned int>  &row_indices,
+               const std::vector<unsigned int>  &col_indices,
+               const FullMatrix<TrilinosScalar> &full_matrix);
+
+                                       /**
+                                        * Set several elements in the
+                                        * specified row of the matrix with
+                                        * column indices as given by
+                                        * <tt>col_indices</tt> to the
+                                        * respective value.
+                                       */
+      void set (const unsigned int                row,
+               const std::vector<unsigned int>   &col_indices,
+               const std::vector<TrilinosScalar> &values);
+
+                                       /**
+                                        * Set several elements to values
+                                        * given by <tt>values</tt> at
+                                        * locations specified by row_indices
+                                        * and col_indices in the sparse
+                                        * matrix. The array <tt>values</tt>
+                                        * 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>i,j</i>).
+                                        */
+      void add (const unsigned int i,
+                const unsigned int j,
+                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.
+                                       */
+      void add (const std::vector<unsigned int>  &row_indices,
+               const std::vector<unsigned int>  &col_indices,
+               const FullMatrix<TrilinosScalar> &full_matrix);
+
+                                       /**
+                                        * Set several elements in the
+                                        * specified row of the matrix with
+                                        * column indices as given by
+                                        * <tt>col_indices</tt> to the
+                                        * respective value.
+                                       */
+      void add (const unsigned int                row,
+               const std::vector<unsigned int>   &col_indices,
+               const std::vector<TrilinosScalar> &values);
+
+                                       /**
+                                        * Set several elements to values
+                                        * given by <tt>values</tt> at
+                                        * locations specified by row_indices
+                                        * and col_indices in the sparse
+                                        * matrix. The array <tt>values</tt>
+                                        * 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<unsigned int> block_col_indices, local_row_length, local_col_indices;
   };
 
 
index e06d4ef2682b2ba6f68bdd206933befa28bfc091..6f66f4233fff924600b1af5d394758e480c07c52 100755 (executable)
@@ -17,6 +17,7 @@
 #include <base/config.h>
 #include <base/subscriptor.h>
 #include <lac/sparse_matrix.h>
+#include <lac/full_matrix.h>
 #include <lac/exceptions.h>
 #include <lac/trilinos_vector_base.h>
 
@@ -866,6 +867,75 @@ namespace TrilinosWrappers
                 const unsigned int j,
                 const TrilinosScalar value);
 
+                                       /**
+                                        * Set all elements given in a
+                                        * FullMatrix<double> 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<unsigned int>  &row_indices,
+               const std::vector<unsigned int>  &col_indices,
+               const FullMatrix<TrilinosScalar> &full_matrix);
+
+                                       /**
+                                        * Set several elements in the
+                                        * specified row of the matrix with
+                                        * column indices as given by
+                                        * <tt>col_indices</tt> 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<unsigned int>   &col_indices,
+               const std::vector<TrilinosScalar> &values);
+
+                                       /**
+                                        * Set several elements to values
+                                        * given by <tt>values</tt> at
+                                        * locations specified by row_indices
+                                        * and col_indices in the sparse
+                                        * matrix. The array <tt>values</tt>
+                                        * 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>i,j</i>).
@@ -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
+                                        * <tt>full_matrix</tt> (cell matrix)
+                                        * into the sparse matrix, at
+                                        * locations according to row_indices
+                                        * and col_indices.
+                                       *
+                                       * Just as the respective call in
+                                       * deal.II SparseMatrix<Number>
+                                       * 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<unsigned int>  &row_indices,
+               const std::vector<unsigned int>  &col_indices,
+               const FullMatrix<TrilinosScalar> &full_matrix);
+
+                                       /**
+                                        * Set several elements in the
+                                        * specified row of the matrix with
+                                        * column indices as given by
+                                        * <tt>col_indices</tt> to the
+                                        * respective value.
+                                       *
+                                       * Just as the respective call in
+                                       * deal.II SparseMatrix<Number>
+                                       * 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<unsigned int>   &col_indices,
+               const std::vector<TrilinosScalar> &values);
+
+                                       /**
+                                        * Set several elements to values
+                                        * given by <tt>values</tt> at
+                                        * locations specified by row_indices
+                                        * and col_indices in the sparse
+                                        * matrix. The array <tt>values</tt>
+                                        * is assumed to store data in C
+                                        * style, i.e., row-wise.
+                                       *
+                                       * Just as the respective call in
+                                       * deal.II SparseMatrix<Number>
+                                       * 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.
index 8fc610188087486b351473788b92b61557495b23..55169ba5694f5c6ebfd521b33580ee14517cef67 100644 (file)
@@ -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<unsigned int>  &row_indices,
+                         const std::vector<unsigned int>  &col_indices,
+                         const FullMatrix<TrilinosScalar> &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<unsigned int>   &col_indices,
+                         const std::vector<TrilinosScalar> &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; i<this->n_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<n_cols; ++j, ++row_length)
+       {
+         const std::pair<unsigned int, unsigned int>
+           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; i<this->n_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<n_rows; ++i)
+      {
+       const std::pair<unsigned int,unsigned int> 
+         row_index = this->row_block_indices.global_to_local (row_indices[i]);
+       for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
+         {
+           if (local_row_length[block_col] == 0)
+             continue;
+
+           block(row_index.first, block_col).set 
+             (1, &row_index.second, 
+              local_row_length[block_col], 
+              &local_col_indices[block_col_indices[block_col]],
+              &values[n_cols*i + block_col_indices[block_col]]);
+         }
+      }
+  }
+
+
+
+  void
+  BlockSparseMatrix::add (const unsigned int   i,
+                         const unsigned int   j,
+                         const TrilinosScalar value)
+  {
+                                  // For adding one single element, it is
+                                  // faster to rely on the BaseClass add
+                                  // function, than doing all the strange
+                                  // operations below.
+    BaseClass::add (i, j, value);
+  }
+
+
+
+  void
+  BlockSparseMatrix::add (const std::vector<unsigned int>  &row_indices,
+                         const std::vector<unsigned int>  &col_indices,
+                         const FullMatrix<TrilinosScalar> &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<unsigned int>   &col_indices,
+                         const std::vector<TrilinosScalar> &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; i<this->n_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<n_cols; ++j, ++row_length)
+       {
+         const std::pair<unsigned int, unsigned int>
+           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; i<this->n_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<n_rows; ++i)
+      {
+       const std::pair<unsigned int,unsigned int> 
+         row_index = this->row_block_indices.global_to_local (row_indices[i]);
+       for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
+         {
+           if (local_row_length[block_col] == 0)
+             continue;
+
+           block(row_index.first, block_col).add 
+             (1, &row_index.second, 
+              local_row_length[block_col], 
+              &local_col_indices[block_col_indices[block_col]],
+              &values[n_cols*i + block_col_indices[block_col]]);
+         }
+      }
+  }
+
+
+
   TrilinosScalar
   BlockSparseMatrix::residual (MPI::BlockVector       &dst,
                               const MPI::BlockVector &x,
index 6e80264f5f7d336a0877165c89c7188ed6da8d35..1802a5a9e340781a7f8df33bae2af405fbfcd197 100755 (executable)
@@ -669,147 +669,233 @@ namespace TrilinosWrappers
            ExcMessage("The given value is not finite but either "
                       "infinite or Not A Number (NaN)"));
 
-                                       // If the last thing we did was an
-                                       // addition, we need to do all the
-                                       // data exchange.
+    set (1, &i, 1, &j, &value);
+  }
+
+
+
+  void
+  SparseMatrix::set (const std::vector<unsigned int>  &row_indices,
+                    const std::vector<unsigned int>  &col_indices,
+                    const FullMatrix<TrilinosScalar> &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<unsigned int>   &col_indices,
+                    const std::vector<TrilinosScalar> &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; i<n_rows; ++i)
       {
-                                       // 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)
+                                  // 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<TrilinosScalar*>(&value),
-                                          (int*)&j);
+           if (matrix->Filled() == false)
+             {
+               ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
+                                  (int)row_indices[i], (int)n_cols,
+                                  const_cast<TrilinosScalar*>(&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<TrilinosScalar*>(&values[i*n_cols]),
+                                  (int*)&col_indices[0]);
          }
        else
          {
-           ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues ( 
-                                          (int)i, 1,
-                                          const_cast<TrilinosScalar*>(&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<unsigned int>  &row_indices,
+                    const std::vector<unsigned int>  &col_indices,
+                    const FullMatrix<TrilinosScalar> &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<unsigned int>   &col_indices,
+                    const std::vector<TrilinosScalar> &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; i<n_rows; ++i)
+      {
                                   // If the calling matrix owns the row to
-                                  // which we want to add this value, we
+                                  // 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::SumIntoGlobalValues(
-                                          (int)i, 1,
-                                          const_cast<TrilinosScalar*>(&value),
-                                          (int*)&j);
-      }
-    else
-      {
+       if (row_map.MyGID(i) == true)
+         {
+           ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(
+                                  (int)row_indices[i], (int)n_cols,
+                                  const_cast<TrilinosScalar*>(&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);

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.