]> https://gitweb.dealii.org/ - dealii.git/commitdiff
When initializing a sparse decomposition (ILU etc) with an additional sparsity patter...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 5 Feb 2010 11:59:37 +0000 (11:59 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 5 Feb 2010 11:59:37 +0000 (11:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@20504 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_decomposition.h
deal.II/lac/include/lac/sparse_decomposition.templates.h

index 32ec0de4161e1255763a4741417d28f3e339eb0c..a946d17481d668c8237a238491cd83454c1a259c 100644 (file)
@@ -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;
     };
index 9ed89e7ecd3329d596064dfd17229cac676298e0..f9937f111595570b7f2d7d52d5319d19fc87389d 100644 (file)
 #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++];
+       }
+    }
 }
 
 

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.