From c0c1b903e2bd7446ae5e83b1f70002faf1cc955a Mon Sep 17 00:00:00 2001
From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 20 Jan 2014 20:01:30 +0000
Subject: [PATCH] Add extract_row_copy method to chunk sparse matrix since
 iterators are very slow

git-svn-id: https://svn.dealii.org/trunk@32267 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../include/deal.II/lac/chunk_sparse_matrix.h | 22 ++++++++++---
 .../lac/chunk_sparse_matrix.templates.h       | 32 +++++++++++++++++++
 2 files changed, 49 insertions(+), 5 deletions(-)

diff --git a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h
index e9ef398417..f44a0318d0 100644
--- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h
@@ -791,11 +791,8 @@ public:
              const size_type j) const;
 
   /**
-   * Return the main diagonal
-   * element in the <i>i</i>th
-   * row. This function throws an
-   * error if the matrix is not
-   * quadratic.
+   * Return the main diagonal element in the <i>i</i>th row. This function
+   * throws an error if the matrix is not quadratic.
    *
    * This function is considerably faster than the operator()(), since for
    * quadratic matrices, the diagonal entry may be the first to be stored in
@@ -810,6 +807,21 @@ public:
    */
   number &diag_element (const size_type i);
 
+  /**
+   * Extracts a copy of the values and indices in the given matrix row.
+   *
+   * The user is expected to pass the length of the arrays column_indices and
+   * values, which gives a means for checking that we do not write to
+   * unallocated memory. This method is motivated by a similar method in
+   * Trilinos row matrices and gives faster access to entries in the matrix as
+   * compared to iterators which are quite slow for this matrix type.
+   */
+  void extract_row_copy (const size_type row,
+                         const size_type array_length,
+                         size_type      &row_length,
+                         size_type      *column_indices,
+                         number         *values) const;
+
 //@}
   /**
    * @name Matrix vector multiplications
diff --git a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
index 62f8409c8f..b493ed3f26 100644
--- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
@@ -596,6 +596,38 @@ ChunkSparseMatrix<number>::add (const number factor,
 }
 
 
+
+template <typename number>
+void
+ChunkSparseMatrix<number>::extract_row_copy (const size_type row,
+                                             const size_type array_length,
+                                             size_type &row_length,
+                                             size_type *column_indices,
+                                             number *values) const
+{
+  row_length = cols->row_length(row);
+  const unsigned int chunk_size = cols->get_chunk_size();
+  const size_type reduced_row = row/chunk_size;
+  AssertIndexRange(cols->sparsity_pattern.row_length(reduced_row)*chunk_size,
+                   array_length+1);
+
+  SparsityPattern::iterator it = cols->sparsity_pattern.begin(reduced_row),
+    itend = cols->sparsity_pattern.end(reduced_row);
+  const number *val_ptr = &val[(it-cols->sparsity_pattern.begin(0))*chunk_size*chunk_size
+                               +(row%chunk_size)*chunk_size];
+  for ( ; it < itend; ++it)
+    {
+      for (unsigned int c=0; c<chunk_size; ++c)
+        {
+          *values++ = val_ptr[c];
+          *column_indices++ = it->column()*chunk_size+c;
+        }
+      val_ptr += chunk_size*chunk_size;
+    }
+}
+
+
+
 template <typename number>
 template <class OutVector, class InVector>
 void
-- 
2.39.5