From: kronbichler Date: Thu, 16 Apr 2009 11:30:47 +0000 (+0000) Subject: Make SparseILU::vmult faster. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8302e66d498193c89ca70b1e778fb2a3ab71ab6a;p=dealii-svn.git Make SparseILU::vmult faster. git-svn-id: https://svn.dealii.org/trunk@18622 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 851560f9a6..f25798f2b3 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -593,6 +593,16 @@ inconvenience this causes.

lac

    +
  1. +

    + Updated: The SparseILU::vmult kernel has been re-written using the same + data structures as in SparseMatrix::vmult, which reduces the count of + operations by one third and the execution times on typical problems by + ten to twenty precent. +
    + (Martin Kronbichler 2009/04/16) +

    +

  2. New: There is now a new class VectorView that allows views of @@ -610,7 +620,7 @@ inconvenience this causes.

  3. Updated: The local_to_global functions in ConstraintMatrix got smarter, - which imcreases the speed of sparsity pattern generation with + which increases the speed of sparsity pattern generation with ConstraintMatrix argument, and makes writing into sparse matrices using distribute_local_to_global faster, in case there are many constraints.
    diff --git a/deal.II/lac/include/lac/sparse_decomposition.templates.h b/deal.II/lac/include/lac/sparse_decomposition.templates.h index 9050076bad..42b0b536fa 100644 --- a/deal.II/lac/include/lac/sparse_decomposition.templates.h +++ b/deal.II/lac/include/lac/sparse_decomposition.templates.h @@ -178,6 +178,18 @@ template void SparseLUDecomposition::copy_from (const SparseMatrix& matrix) { + // check whether we use the same sparsity + // pattern as the input matrix + if (&this->get_sparsity_pattern() == &matrix.get_sparsity_pattern()) + { + 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; + return; + } + // preset the elements std::fill_n (&this->global_entry(0), this->n_nonzero_elements(), diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index c7cd70031b..c19807c119 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -54,7 +54,6 @@ template void SparseILU::decompose (const SparseMatrix &matrix, const double strengthen_diagonal) { - SparseLUDecomposition::decompose (matrix, strengthen_diagonal); Assert (matrix.m()==matrix.n(), ExcNotQuadratic ()); Assert (this->m()==this->n(), ExcNotQuadratic ()); Assert (matrix.m()==this->m(), ExcDimensionMismatch(matrix.m(), this->m())); @@ -62,7 +61,7 @@ void SparseILU::decompose (const SparseMatrix &matrix, Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal)); - this->copy_from (matrix); + SparseLUDecomposition::decompose (matrix, strengthen_diagonal); if (strengthen_diagonal>0) this->strengthen_diagonal_impl(); @@ -77,7 +76,7 @@ void SparseILU::decompose (const SparseMatrix &matrix, const std::size_t * const ia = sparsity.get_rowstart_indices(); const unsigned int * const ja = sparsity.get_column_numbers(); - number * luval = &this->global_entry (0); + number * luval = this->SparseMatrix::val; const unsigned int N = this->m(); @@ -178,6 +177,7 @@ void SparseILU::vmult (Vector &dst, = this->get_sparsity_pattern().get_rowstart_indices(); const unsigned int * const column_numbers = this->get_sparsity_pattern().get_column_numbers(); + // solve LUx=b in two steps: // first Ly = b, then // Ux = y @@ -199,9 +199,13 @@ void SparseILU::vmult (Vector &dst, // find the position where the part // right of the diagonal starts const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row]; - - for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col) - dst(row) -= this->global_entry (col-column_numbers) * dst(*col); + + somenumber dst_row = 0; + const number * luval = this->SparseMatrix::val + + (rowstart - column_numbers); + for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col, ++luval) + dst_row += *luval * dst(*col); + dst(row) -= dst_row; } // now the backward solve. same @@ -220,8 +224,13 @@ void SparseILU::vmult (Vector &dst, // right of the diagonal starts const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row]; - for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col) - dst(row) -= this->global_entry (col-column_numbers) * dst(*col); + somenumber dst_row = 0; + const number * luval = this->SparseMatrix::val + + (first_after_diagonal - column_numbers); + for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col, ++luval) + dst_row += *luval * dst(*col); + + dst(row) -= dst_row; // scale by the diagonal element. // note that the diagonal element @@ -244,6 +253,7 @@ void SparseILU::Tvmult (Vector &dst, = this->get_sparsity_pattern().get_rowstart_indices(); const unsigned int * const column_numbers = this->get_sparsity_pattern().get_column_numbers(); + // solve (LU)'x=b in two steps: // first U'y = b, then // L'x = y @@ -268,9 +278,12 @@ void SparseILU::Tvmult (Vector &dst, // find the position where the part // right of the diagonal starts const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row]; - - for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col) - tmp(*col) += this->global_entry (col-column_numbers) * dst(row); + + const somenumber dst_row = dst (row); + const number * luval = this->SparseMatrix::val + + (first_after_diagonal - column_numbers); + for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col, ++luval) + tmp(*col) += *luval * dst_row; } // now the backward solve. same @@ -292,9 +305,12 @@ void SparseILU::Tvmult (Vector &dst, // find the position where the part // right of the diagonal starts const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row]; - - for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col) - tmp(*col) += this->global_entry (col-column_numbers) * dst(row); + + const somenumber dst_row = dst (row); + const number * luval = this->SparseMatrix::val + + (rowstart - column_numbers); + for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col, ++luval) + tmp(*col) += *luval * dst_row; } }