]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
The deal.II SparseMatrix class can now set a FullMatrix of data at once, as can all...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 26 Nov 2008 19:21:53 +0000 (19:21 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 26 Nov 2008 19:21:53 +0000 (19:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@17754 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/block_matrix_base.h
deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h

index fe842cbff735c233adf2170d5f335066976bb986..b408f40742ac46105c41b3558c30783bc648b2a5 100644 (file)
@@ -19,6 +19,7 @@
 #include <base/smartpointer.h>
 #include <lac/block_indices.h>
 #include <lac/exceptions.h>
+#include <lac/full_matrix.h>
 #include <lac/matrix_iterator.h>
 
 #include <cmath>
@@ -477,21 +478,196 @@ class BlockMatrixBase : public Subscriptor
     void set (const unsigned int i,
              const unsigned int j,
              const value_type value);
-    
+
+                                     /**
+                                     * Set all elements given in a
+                                     * FullMatrix 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.
+                                     *
+                                     * The optional parameter
+                                     * <tt>elide_zero_values</tt> can be
+                                     * used to specify whether zero
+                                     * values should be set anyway or
+                                     * they should be filtered away (and
+                                     * not change the previous content in
+                                     * the respective element if it
+                                     * exists). The default value is
+                                     * <tt>false</tt>, i.e., even zero
+                                     * values are treated.
+                                     */
+    void set (const std::vector<unsigned int> &indices,
+             const FullMatrix<value_type>    &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<value_type>    &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.
+                                     *
+                                     * The optional parameter
+                                     * <tt>elide_zero_values</tt> can be
+                                     * used to specify whether zero
+                                     * values should be set anyway or
+                                     * they should be filtered away (and
+                                     * not change the previous content in
+                                     * the respective element if it
+                                     * exists). The default value is
+                                     * <tt>false</tt>, i.e., even zero
+                                     * values are treated.
+                                     */
+    void set (const unsigned int               row,
+             const std::vector<unsigned int> &col_indices,
+             const std::vector<value_type>   &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.
+                                     *
+                                     * 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 value_type   *values,
+             const bool          elide_zero_values = false);
+
                                     /**
                                      * Add <tt>value</tt> to the
-                                     * element <tt>(i,j)</tt>.
-                                     * Throws an error if the entry
-                                     * does not exist or if
-                                     * <tt>value</tt> is not a finite
-                                     * number. Still, it is allowed
-                                     * to store zero values in
-                                     * non-existent fields.
+                                     * element (<i>i,j</i>).  Throws
+                                     * an error if the entry does not
+                                     * exist or if <tt>value</tt> is
+                                     * not a finite number. Still, it
+                                     * is allowed to store zero
+                                     * values in non-existent fields.
                                      */
     void add (const unsigned int i,
               const unsigned int j,
              const value_type 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.
+                                     *
+                                     * 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<value_type>    &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<value_type>    &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.
+                                     *
+                                     * 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<value_type>   &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.
+                                     *
+                                     * 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 value_type   *values,
+             const bool          elide_zero_values = true);
+
                                     /**
                                      * Return the value of the entry
                                      * (i,j).  This may be an
@@ -984,7 +1160,30 @@ class BlockMatrixBase : public Subscriptor
     friend class ConstAccessor;
     
     friend class const_iterator;
-#endif    
+#endif
+
+
+  private:
+                                       /**
+                                       * Temporary vector for block indices
+                                       * when writing local to global data.
+                                       */
+    std::vector<unsigned int> block_col_indices; 
+
+                                       /**
+                                       * Temporary vector for column
+                                       * indices when writing local to
+                                       * global data and passing the call
+                                       * down to a SparseMatrix.
+                                       */
+    std::vector<unsigned int> local_col_indices;
+
+                                       /**
+                                       * Temporary vector for storing the
+                                       * local lengths of the row data.
+                                       */
+    std::vector<unsigned int> local_row_length;
+
 };
 
 
@@ -1554,6 +1753,9 @@ BlockMatrixBase<MatrixType>::n_block_rows () const
 
 
 
+                                   // Write the single set manually,
+                                   // since the other function has a lot
+                                   // of overhead in that case.
 template <class MatrixType>
 inline
 void
@@ -1575,6 +1777,144 @@ BlockMatrixBase<MatrixType>::set (const unsigned int i,
 
 
 
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::set (const std::vector<unsigned int> &row_indices,
+                                 const std::vector<unsigned int> &col_indices,
+                                 const FullMatrix<value_type>    &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);
+}
+
+
+
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::set (const std::vector<unsigned int> &indices,
+                                 const FullMatrix<value_type>    &values,
+                                 const bool                       elide_zero_values)
+{
+  Assert (indices.size() == values.m(),
+         ExcDimensionMismatch(indices.size(), values.m()));
+  Assert (values.n() == values.m(), ExcNotQuadratic());
+
+  for (unsigned int i=0; i<indices.size(); ++i)
+    set (indices[i], indices.size(), &indices[0], &values(i,0),
+        elide_zero_values);
+}
+
+
+
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::set (const unsigned int               row,
+                                 const std::vector<unsigned int> &col_indices,
+                                 const std::vector<value_type>   &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);
+}
+
+
+
+                                  // This is a very messy function, since
+                                  // we need to calculate to each position 
+                                  // the location in the global array.
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::set (const unsigned int   row,
+                                 const unsigned int   n_cols,
+                                 const unsigned int  *col_indices,
+                                 const value_type    *values,
+                                 const bool           elide_zero_values)
+{
+                                  // Resize scratch arrays
+  local_row_length.resize (this->n_block_cols());
+  block_col_indices.resize (this->n_block_cols());
+  local_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.
+  {
+    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!
+  const std::pair<unsigned int,unsigned int> 
+    row_index = this->row_block_indices.global_to_local (row);
+  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 
+       (row_index.second, 
+        local_row_length[block_col], 
+        &local_col_indices[block_col_indices[block_col]],
+        &values[block_col_indices[block_col]],
+        elide_zero_values);
+    }
+}
+
+
+
 template <class MatrixType>
 inline
 void
@@ -1605,6 +1945,148 @@ BlockMatrixBase<MatrixType>::add (const unsigned int i,
 
 
 
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::add (const std::vector<unsigned int> &row_indices,
+                                 const std::vector<unsigned int> &col_indices,
+                                 const FullMatrix<value_type>    &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);
+}
+
+
+
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::add (const std::vector<unsigned int> &indices,
+                                 const FullMatrix<value_type>    &values,
+                                 const bool                       elide_zero_values)
+{
+  Assert (indices.size() == values.m(),
+         ExcDimensionMismatch(indices.size(), values.m()));
+  Assert (values.n() == values.m(), ExcNotQuadratic());
+
+  for (unsigned int i=0; i<indices.size(); ++i)
+    add (indices[i], indices.size(), &indices[0], &values(i,0),
+        elide_zero_values);
+}
+
+
+
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::add (const unsigned int               row,
+                                 const std::vector<unsigned int> &col_indices,
+                                 const std::vector<value_type>   &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);
+}
+
+
+
+                                  // This is a very messy function, since
+                                  // we need to calculate to each position 
+                                  // the location in the global array.
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::add (const unsigned int   row,
+                                 const unsigned int   n_cols,
+                                 const unsigned int  *col_indices,
+                                 const value_type    *values,
+                                 const bool           elide_zero_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.
+  {
+    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!
+  const std::pair<unsigned int,unsigned int> 
+    row_index = this->row_block_indices.global_to_local (row);
+  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 
+       (row_index.second, 
+        local_row_length[block_col], 
+        &local_col_indices[block_col_indices[block_col]],
+        &values[block_col_indices[block_col]],
+        elide_zero_values);
+    }
+}
+
+
+
 template <class MatrixType>
 inline
 typename BlockMatrixBase<MatrixType>::value_type
index 7436e1b037bbdfd336e7cad0bdc9bc753f94b583..14c92050805dffc01b81f25603c81f7e532ddd15 100644 (file)
@@ -1444,8 +1444,7 @@ namespace PETScWrappers
     else
       {
                                   // Otherwise, extract nonzero values in
-                                  // each row and pass on to the other 
-                                  // function.
+                                  // each row and get the respective index.
        column_indices.resize(n_cols);
        column_values.resize(n_cols);
 
@@ -1609,8 +1608,7 @@ namespace PETScWrappers
     else
       {
                                   // Otherwise, extract nonzero values in
-                                  // each row and pass on to the other 
-                                  // function.
+                                  // each row and get the respective index.
        column_indices.resize(n_cols);
        column_values.resize(n_cols);
 
index da0a9c31475dd23d65a94c0441deb6425b72c43f..16f04091b73a96d60b90c20037d9d118eed0a028 100644 (file)
@@ -765,7 +765,95 @@ class SparseMatrix : public virtual Subscriptor
     void set (const unsigned int i,
               const unsigned int j,
              const number value);
-    
+
+                                       /**
+                                        * Set all elements given in a
+                                        * FullMatrix 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.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be set anyway or
+                                       * they should be filtered away (and
+                                       * not change the previous content in
+                                       * the respective element if it
+                                       * exists). The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are treated.
+                                       */
+    void set (const std::vector<unsigned int> &indices,
+             const FullMatrix<number>        &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<number>        &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.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be set anyway or
+                                       * they should be filtered away (and
+                                       * not change the previous content in
+                                       * the respective element if it
+                                       * exists). The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are treated.
+                                       */
+    void set (const unsigned int               row,
+             const std::vector<unsigned int> &col_indices,
+             const std::vector<number>       &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.
+                                       *
+                                       * 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 number       *values,
+             const bool          elide_zero_values = false);
+
                                     /**
                                      * Add <tt>value</tt> to the
                                      * element (<i>i,j</i>).  Throws
@@ -779,6 +867,94 @@ class SparseMatrix : public virtual Subscriptor
               const unsigned int j,
              const number 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.
+                                       *
+                                       * 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<number>        &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<number>        &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.
+                                       *
+                                       * 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<number>       &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.
+                                       *
+                                       * 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 number       *values,
+             const bool          elide_zero_values = true);
+
                                     /**
                                      * Multiply the entire matrix by a
                                      * fixed factor.
@@ -1713,6 +1889,26 @@ class SparseMatrix : public virtual Subscriptor
                                      */
     unsigned int max_len;
 
+                                      /**
+                                       * An internal array of integer
+                                       * values that is used to store the
+                                       * indices in the global value array
+                                       * val when adding/inserting local
+                                       * data into the (large) sparse
+                                       * matrix.
+                                       */
+    std::vector<unsigned int> global_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<number> column_values;
+
+
                                     /**
                                      * Version of vmult() which only
                                      * performs its actions on the
@@ -1812,51 +2008,290 @@ unsigned int SparseMatrix<number>::n () const
 }
 
 
-template <typename number>
+                                       // 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.
+template <typename number>  
+inline
+void
+SparseMatrix<number>::set (const unsigned int i,
+                          const unsigned int j,
+                          const number       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);
+}
+
+
+
+template <typename number>  
+inline
+void
+SparseMatrix<number>::set (const std::vector<unsigned int> &indices,
+                          const FullMatrix<number>        &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);
+}
+
+
+
+template <typename number>  
+inline
+void
+SparseMatrix<number>::set (const std::vector<unsigned int> &row_indices,
+                          const std::vector<unsigned int> &col_indices,
+                          const FullMatrix<number>        &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);
+}
+
+
+
+template <typename number>  
 inline
-void SparseMatrix<number>::set (const unsigned int i,
-                               const unsigned int j,
-                               const number value)
+void
+SparseMatrix<number>::set (const unsigned int               row,
+                          const std::vector<unsigned int> &col_indices,
+                          const std::vector<number>       &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);
+}
+
 
-  Assert (numbers::is_finite(value), 
-          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
+template <typename number>  
+inline
+void
+SparseMatrix<number>::set (const unsigned int  row,
+                          const unsigned int  n_cols,
+                          const unsigned int *col_indices,
+                          const number       *values,
+                          const bool          elide_zero_values)
+{
   Assert (cols != 0, ExcNotInitialized());
-                                  // it is allowed to set elements of
+
+  unsigned int n_columns = 0;
+
+                                  // Otherwise, extract nonzero values in
+                                  // each row and pass on to the other 
+                                  // function.
+  global_indices.resize(n_cols);
+  column_values.resize(n_cols);
+
+  const unsigned int * index_ptr = &col_indices[0];
+  const number * value_ptr = &values[0];
+
+                                  // First, search all the indices to find
+                                  // out which values we actually need to
+                                  // set.
+                                  // TODO: Could probably made more
+                                  // efficient.
+  for (unsigned int j=0; j<n_cols; ++j)
+    {
+      const number value = *value_ptr++;
+      Assert (numbers::is_finite(value),
+             ExcMessage("The given value is not finite but either "
+                        "infinite or Not A Number (NaN)"));
+
+      if (value == 0 && elide_zero_values == true)
+       {
+         continue;
+         ++index_ptr;
+       }
+
+      const unsigned int index = cols->operator()(row, *index_ptr++);
+
+                                  // it is allowed to set elements in
                                   // the matrix that are not part of
                                   // the sparsity pattern, if the
                                   // value to which we set it is zero
-  const unsigned int index = cols->operator()(i,j);
-  Assert ((index != SparsityPattern::invalid_entry) ||
-         (value == 0.),
-         ExcInvalidIndex(i,j));
+      if (index == SparsityPattern::invalid_entry)
+       {
+         Assert ((index != SparsityPattern::invalid_entry) ||
+                 (value == 0.),
+                 ExcInvalidIndex(row,col_indices[j]));
+         continue;
+       }
+
+      global_indices[n_columns] = index;
+      column_values[n_columns] = value;
+      n_columns++;
+    }
 
-  if (index != SparsityPattern::invalid_entry)
-    val[index] = value;
+  index_ptr = &global_indices[0];
+  value_ptr = &column_values[0];
+  
+                                   // Finally, go through the index list 
+                                   // and set the elements one by one.
+  for (unsigned int j=0; j<n_columns; ++j)
+    val[*index_ptr++] = *value_ptr++;
 }
 
 
 
-template <typename number>
+template <typename number>  
+inline
+void
+SparseMatrix<number>::add (const unsigned int i,
+                          const unsigned int j,
+                          const number       value)
+{
+
+  Assert (numbers::is_finite(value),
+         ExcMessage("The given value is not finite but either "
+                    "infinite or Not A Number (NaN)"));
+
+  add (i, 1, &j, &value, false);
+}
+
+
+
+template <typename number>  
+inline
+void
+SparseMatrix<number>::add (const std::vector<unsigned int> &indices,
+                          const FullMatrix<number>        &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);
+}
+
+
+
+template <typename number>  
+inline
+void
+SparseMatrix<number>::add (const std::vector<unsigned int> &row_indices,
+                          const std::vector<unsigned int> &col_indices,
+                          const FullMatrix<number>        &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);
+}
+
+
+
+template <typename number>  
 inline
-void SparseMatrix<number>::add (const unsigned int i,
-                               const unsigned int j,
-                               const number value)
+void
+SparseMatrix<number>::add (const unsigned int               row,
+                          const std::vector<unsigned int> &col_indices,
+                          const std::vector<number>       &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);
+}
+
 
-  Assert (numbers::is_finite(value), 
-          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
+template <typename number>  
+inline
+void
+SparseMatrix<number>::add (const unsigned int  row,
+                          const unsigned int  n_cols,
+                          const unsigned int *col_indices,
+                          const number       *values,
+                          const bool          elide_zero_values)
+{
   Assert (cols != 0, ExcNotInitialized());
 
-  const unsigned int index = cols->operator()(i,j);
-  Assert ((index != SparsityPattern::invalid_entry) ||
-         (value == 0.),
-         ExcInvalidIndex(i,j));
+  unsigned int n_columns = 0;
 
-  if (value != 0.)
-    val[index] += value;
+  global_indices.resize(n_cols);
+  column_values.resize(n_cols);
+
+  const unsigned int * index_ptr = &col_indices[0];
+  const number * value_ptr = &values[0];
+
+                                  // First, search all the indices to find
+                                  // out which values we actually need to
+                                  // add.
+                                  // TODO: Could probably made more
+                                  // efficient.
+  for (unsigned int j=0; j<n_cols; ++j)
+    {
+      const number value = *value_ptr++;
+      Assert (numbers::is_finite(value),
+             ExcMessage("The given value is not finite but either "
+                        "infinite or Not A Number (NaN)"));
+
+      if (value == 0 && elide_zero_values == true)
+       {
+         continue;
+         ++index_ptr;
+       }
+
+      const unsigned int index = cols->operator()(row, *index_ptr++);
+
+                                  // it is allowed to add elements to
+                                  // the matrix that are not part of
+                                  // the sparsity pattern, if the
+                                  // value to which we set it is zero
+      if (index == SparsityPattern::invalid_entry)
+       {
+         Assert ((index != SparsityPattern::invalid_entry) ||
+                 (value == 0.),
+                 ExcInvalidIndex(row,col_indices[j]));
+         continue;
+       }
+
+      global_indices[n_columns] = index;
+      column_values[n_columns] = value;
+      n_columns++;
+    }
+
+  index_ptr = &global_indices[0];
+  value_ptr = &column_values[0];
+  
+                                   // Finally, go through the index list 
+                                   // and add the elements one by one.
+  for (unsigned int j=0; j<n_columns; ++j)
+    val[*index_ptr++] += *value_ptr++;
 }
 
 
index 6aaea7c2ad5093891ded60cda34633943584a969..973a204fea7574864c89c05da25cece975766304 100644 (file)
@@ -292,152 +292,6 @@ 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.
-                                       *
-                                       * 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>  &row_indices,
-               const std::vector<unsigned int>  &col_indices,
-               const FullMatrix<TrilinosScalar> &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.
-                                       *
-                                       * 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<TrilinosScalar> &values,
-               const bool                         elide_zero_values = false);
-
-                                       /**
-                                        * Set several elements to values
-                                        * given by <tt>values</tt> in the
-                                        * given row at column locations
-                                        * specified by col_indices in the
-                                        * sparse matrix.
-                                       *
-                                       * 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 TrilinosScalar *values,
-               const bool            elide_zero_values = false);
-
-                                       /**
-                                        * 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.
-                                       *
-                                       * 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>  &row_indices,
-               const std::vector<unsigned int>  &col_indices,
-               const FullMatrix<TrilinosScalar> &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.
-                                       *
-                                       * 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<TrilinosScalar> &values,
-               const bool                         elide_zero_values = true);
-
-                                       /**
-                                        * Set several elements to values
-                                        * given by <tt>values</tt> in the
-                                        * given row at column locations
-                                        * specified by col_indices in the
-                                        * sparse matrix.
-                                       *
-                                       * 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 TrilinosScalar *values,
-               const bool            elide_zero_values = true);
-
                                        /**
                                         * Matrix-vector multiplication:
                                         * let $dst = M*src$ with $M$
@@ -770,14 +624,6 @@ 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_col_indices, local_row_length;
   };
 
 
@@ -804,267 +650,6 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  void
-  BlockSparseMatrix::set (const unsigned int   i,
-                         const unsigned int   j,
-                         const TrilinosScalar value)
-  {
-    BaseClass::set (i, j, value);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::set (const std::vector<unsigned int>  &row_indices,
-                         const std::vector<unsigned int>  &col_indices,
-                         const FullMatrix<TrilinosScalar> &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
-  BlockSparseMatrix::set (const unsigned int                 row,
-                         const std::vector<unsigned int>   &col_indices,
-                         const std::vector<TrilinosScalar> &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);
-  }
-
-
-
-                                  // This is a very messy function, since
-                                  // we need to calculate to each position 
-                                  // the location in the global array.
-  inline
-  void
-  BlockSparseMatrix::set (const unsigned int    row,
-                         const unsigned int    n_cols,
-                         const unsigned int   *col_indices,
-                         const TrilinosScalar *values,
-                         const bool            elide_zero_values)
-  {
-                                  // Resize scratch arrays
-    local_row_length.resize (this->n_block_cols());
-    block_col_indices.resize (this->n_block_cols());
-    local_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.
-    {
-      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!
-    const std::pair<unsigned int,unsigned int> 
-      row_index = this->row_block_indices.global_to_local (row);
-    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 
-             (row_index.second, 
-              local_row_length[block_col], 
-              &local_col_indices[block_col_indices[block_col]],
-              &values[block_col_indices[block_col]],
-              elide_zero_values);
-      }
-  }
-
-
-
-  inline
-  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);
-  }
-
-
-
-  inline
-  void
-  BlockSparseMatrix::add (const std::vector<unsigned int>  &row_indices,
-                         const std::vector<unsigned int>  &col_indices,
-                         const FullMatrix<TrilinosScalar> &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
-  BlockSparseMatrix::add (const unsigned int                 row,
-                         const std::vector<unsigned int>   &col_indices,
-                         const std::vector<TrilinosScalar> &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
-  BlockSparseMatrix::add (const unsigned int    row,
-                         const unsigned int    n_cols,
-                         const unsigned int   *col_indices,
-                         const TrilinosScalar *values,
-                         const bool            elide_zero_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.
-    {
-      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!
-    const std::pair<unsigned int,unsigned int> 
-      row_index = this->row_block_indices.global_to_local (row);
-    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 
-             (row_index.second, 
-              local_row_length[block_col], 
-              &local_col_indices[block_col_indices[block_col]],
-              &values[block_col_indices[block_col]],
-              elide_zero_values);
-      }
-  }
-
-
-
   inline
   void
   BlockSparseMatrix::vmult (MPI::BlockVector       &dst,
index 5441490d80d8c4defdf50a09cf4b54e7172badc3..5ac56828d500044496652add7b7712a75de37d31 100755 (executable)
@@ -2168,8 +2168,8 @@ namespace TrilinosWrappers
     else
       {
                                   // Otherwise, extract nonzero values in
-                                  // each row and pass on to the other 
-                                  // function.
+                                  // each row and get the respective
+                                  // indices.
        column_indices.resize(n_cols);
        column_values.resize(n_cols);
 
@@ -2394,8 +2394,7 @@ namespace TrilinosWrappers
     else
       {
                                   // Otherwise, extract nonzero values in
-                                  // each row and pass on to the other 
-                                  // function.
+                                  // each row and the corresponding index.
        column_indices.resize(n_cols);
        column_values.resize(n_cols);
 

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.