add BlockMatrixBase::add(factor,matrix) (patch from Jean-Paul Pelteret)

git-svn-id: https://svn.dealii.org/trunk@28234 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/lac/block_matrix_base.h b/deal.II/include/deal.II/lac/block_matrix_base.h
index d7de0a7eb2bf8d591888da827e0560905db93810..7f7940a879bfa0a525e7ffc90ceb43370f836368 100644 (file)
--- a/deal.II/include/deal.II/lac/block_matrix_base.h
+++ b/deal.II/include/deal.II/lac/block_matrix_base.h
@@ -681,6 +681,21 @@ public:
             const bool          elide_zero_values = true,
             const bool          col_indices_are_sorted = false);
 
+  /**
+   * Add <tt>matrix</tt> scaled by
+   * <tt>factor</tt> to this matrix,
+   * i.e. the matrix
+   * <tt>factor*matrix</tt> is added to
+   * <tt>this</tt>. If the sparsity
+   * pattern of the calling matrix does
+   * not contain all the elements in
+   * the sparsity pattern of the input
+   * matrix, this function will throw
+   * an exception.
+   */
+  void add (const value_type                   factor,
+            const BlockMatrixBase<MatrixType> &matrix);
+
   /**
    * Return the value of the entry
    * (i,j).  This may be an
@@ -758,28 +773,6 @@ public:
    */
   BlockMatrixBase &operator /= (const value_type factor);
 
-  /**
-   * Add <tt>matrix</tt> scaled by
-   * <tt>factor</tt> to this matrix,
-   * i.e. the matrix <tt>factor*matrix</tt>
-   * is added to <tt>this</tt>. This
-   * function throws an error if the
-   * sparsity patterns of the two involved
-   * matrices do not point to the same
-   * object, since in this case the
-   * operation is cheaper.
-   *
-   * The source matrix may be a sparse
-   * matrix over an arbitrary underlying
-   * scalar type, as long as its data type
-   * is convertible to the data type of
-   * this matrix.
-   */
-  template <class BlockMatrixType>
-  void add (const value_type       factor,
-            const BlockMatrixType &matrix);
-
-
   /**
    * Adding Matrix-vector
    * multiplication. Add $M*src$ on
@@ -851,6 +844,19 @@ public:
                        const BlockVectorType &x,
                        const BlockVectorType &b) const;
 
+  /**
+   * Print the matrix to the given
+   * stream, using the format
+   * <tt>(line,col) value</tt>, i.e. one
+   * nonzero entry of the matrix per
+   * line. The optional flag outputs the
+   * sparsity pattern in a different style
+   * according to the underlying
+   * sparsematrix type.
+   */
+  void print (std::ostream &out,
+              const bool    alternative_output = false) const;
+
   /**
    * STL-like iterator with the
    * first entry.
@@ -2266,6 +2272,34 @@ BlockMatrixBase<MatrixType>::add (const unsigned int   row,
 
 
 
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::add (const value_type                   factor,
+                                  const BlockMatrixBase<MatrixType> &matrix)
+{
+  Assert (numbers::is_finite(factor), ExcNumberNotFinite());
+
+  prepare_add_operation();
+
+  // save some cycles for zero additions, but
+  // only if it is safe for the matrix we are
+  // working with
+  typedef typename MatrixType::Traits MatrixTraits;
+  if ((MatrixTraits::zero_addition_can_be_elided == true)
+      &&
+      (factor == 0))
+    return;
+
+  for (unsigned int row=0; row<n_block_rows(); ++row)
+    for (unsigned int col=0; col<n_block_cols(); ++col)
+      // This function should throw if the sparsity
+      // patterns of the two blocks differ
+      block(row, col).add(factor, matrix.block(row,col));
+}
+
+
+
 template <class MatrixType>
 inline
 typename BlockMatrixBase<MatrixType>::value_type
@@ -2684,6 +2718,24 @@ residual (BlockVectorType          &dst,
 
 
 
+template <class MatrixType>
+inline
+void
+BlockMatrixBase<MatrixType>::print (std::ostream &out,
+                                    const bool    alternative_output) const
+{
+  for (unsigned int row=0; row<n_block_rows(); ++row)
+    for (unsigned int col=0; col<n_block_cols(); ++col)
+      {
+        if (!alternative_output)
+          out << "Block (" << row << ", " << col << ")" << std::endl;
+
+        block(row, col).print(out, alternative_output);
+      }
+}
+
+
+
 template <class MatrixType>
 inline
 typename BlockMatrixBase<MatrixType>::const_iterator