From: Martin Kronbichler Date: Fri, 5 Feb 2010 11:59:37 +0000 (+0000) Subject: When initializing a sparse decomposition (ILU etc) with an additional sparsity patter... X-Git-Tag: v8.0.0~6543 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=54a3e87f4bfa37a2cac6a5396962e69780c34ac5;p=dealii.git 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 --- 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, * causing this sparsity to * be used. * - * Note that the sparsity - * structures of - * *use_this_sparsity and + * Note that the sparsity structures + * of *use_this_sparsity 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 +#include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -185,8 +188,11 @@ SparseLUDecomposition::copy_from (const SparseMatrix& 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::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::copy_from (const SparseMatrix& 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; rowm(); ++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++]; + } + } }