From b1de01d4db2ae9d908b766eb171a082f065f4170 Mon Sep 17 00:00:00 2001
From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 5 Feb 2010 11:59:37 +0000
Subject: [PATCH] When initializing a sparse decomposition (ILU etc) with an
 additional sparsity pattern, this sparsity pattern does not need to be a
 superset of the matrix sparsity pattern. This is allows e.g. to initialize an
 ILU with neglecting some zeros (Dirichlet boundary nodes) in the
 decomposition sparsity pattern without copying the matrix.

git-svn-id: https://svn.dealii.org/trunk@20504 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../lac/include/lac/sparse_decomposition.h    | 16 +++-----
 .../lac/sparse_decomposition.templates.h      | 39 +++++++++++++++----
 2 files changed, 37 insertions(+), 18 deletions(-)

diff --git a/deal.II/lac/include/lac/sparse_decomposition.h b/deal.II/lac/include/lac/sparse_decomposition.h
index 32ec0de416..a946d17481 100644
--- a/deal.II/lac/include/lac/sparse_decomposition.h
+++ b/deal.II/lac/include/lac/sparse_decomposition.h
@@ -256,17 +256,13 @@ class SparseLUDecomposition : protected SparseMatrix<number>,
 					  * causing this sparsity to
 					  * be used.
 					  *
-					  * Note that the sparsity
-					  * structures of
-					  * <tt>*use_this_sparsity</tt> and
+					  * Note that the sparsity structures
+					  * of <tt>*use_this_sparsity</tt> and
 					  * the matrix passed to the
-					  * initialize function need
-					  * not be equal, but that the
-					  * pattern used by this
-					  * matrix needs to contain
-					  * all elements used by the
-					  * matrix to be decomposed.
-					  * Fill-in is thus allowed.
+					  * initialize function need not be
+					  * equal. Fill-in is allowed, as well
+					  * as filtering out some elements in
+					  * the matrix.
 					  */
 	const SparsityPattern *use_this_sparsity;
     };
diff --git a/deal.II/lac/include/lac/sparse_decomposition.templates.h b/deal.II/lac/include/lac/sparse_decomposition.templates.h
index 9ed89e7ecd..f9937f1115 100644
--- a/deal.II/lac/include/lac/sparse_decomposition.templates.h
+++ b/deal.II/lac/include/lac/sparse_decomposition.templates.h
@@ -12,9 +12,12 @@
 #ifndef __deal2__sparse_decomposition_templates_h
 #define __deal2__sparse_decomposition_templates_h
 
+
 #include <base/memory_consumption.h>
+#include <base/template_constraints.h>
 #include <lac/sparse_decomposition.h>
 #include <algorithm>
+#include <cstring>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -185,8 +188,11 @@ SparseLUDecomposition<number>::copy_from (const SparseMatrix<somenumber>& matrix
       const somenumber * input_ptr = matrix.val;
       number * this_ptr = this->val;
       const number * const end_ptr = this_ptr + this->n_nonzero_elements();
-      for ( ; this_ptr != end_ptr; ++input_ptr, ++this_ptr)
-	*this_ptr = *input_ptr;
+      if (types_are_equal<somenumber, number>::value == true)
+	std::memcpy (this_ptr, input_ptr, this->n_nonzero_elements()*sizeof(number));
+      else
+	for ( ; this_ptr != end_ptr; ++input_ptr, ++this_ptr)
+	  *this_ptr = *input_ptr;
       return;
     }
 
@@ -197,16 +203,33 @@ SparseLUDecomposition<number>::copy_from (const SparseMatrix<somenumber>& matrix
 
                                    // note: pointers to the sparsity
                                    // pattern of the old matrix!
-  const std::size_t * const rowstart_indices
+  const std::size_t * const in_rowstart_indices
     = matrix.get_sparsity_pattern().get_rowstart_indices();
-
-  const unsigned int * const column_numbers
+  const unsigned int * const in_cols
     = matrix.get_sparsity_pattern().get_column_numbers();
+  const unsigned int * cols = this->get_sparsity_pattern().get_column_numbers();
+  const std::size_t * rowstart_indices = 
+    this->get_sparsity_pattern().get_rowstart_indices();
 
+				   // both allow more and less entries 
+				   // in the new matrix
+  std::size_t in_index = in_rowstart_indices[0], index;
   for (unsigned int row=0; row<this->m(); ++row)
-    this->set (row, rowstart_indices[row+1]-rowstart_indices[row],
-	       &column_numbers[rowstart_indices[row]],
-	       &matrix.val[rowstart_indices[row]], false);
+    {
+      index = rowstart_indices[row];
+      in_index = in_rowstart_indices[row];
+      this->val[index++] = matrix.val[in_index++];
+      while (in_index < in_rowstart_indices[row+1] && 
+	     index < rowstart_indices[row+1])
+	{
+	  while (cols[index] < in_cols[in_index] && index < rowstart_indices[row+1])
+	    ++index;
+	  while (in_cols[in_index] < cols[index] && in_index < in_rowstart_indices[row+1])
+	    ++in_index;
+
+	  this->val[index++] = matrix.val[in_index++];
+	}
+    }
 }
 
 
-- 
2.39.5