From: Wolfgang Bangerth Date: Sat, 1 Jan 2022 11:56:55 +0000 (-0700) Subject: Parallelize more operations in the UMFPACK interfaces. X-Git-Tag: v9.4.0-rc1~671^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6bba4be23eec77066cfc3b2ca985b4f7a79e47fc;p=dealii.git Parallelize more operations in the UMFPACK interfaces. --- diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index a265776409..6807ce6433 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -189,34 +189,42 @@ SparseDirectUMFPACK::sort_arrays(const BlockSparseMatrix &matrix) // first. however, that means that there may be as many entries per row // in the wrong place as there are block columns. we can do the same // thing as above, but we have to do it multiple times - for (size_type row = 0; row < matrix.m(); ++row) - { - long int cursor = Ap[row]; - for (size_type block = 0; block < matrix.n_block_cols(); ++block) + parallel::apply_to_subranges( + 0, + matrix.m(), + [this, &matrix](const size_type row_begin, const size_type row_end) { + for (size_type row = row_begin; row < row_end; ++row) { - // find the next out-of-order element - while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] < Ai[cursor + 1])) - ++cursor; - - // if there is none, then just go on - if (cursor == Ap[row + 1] - 1) - break; - - // otherwise swap this entry with successive ones as long as - // necessary - long int element = cursor; - while ((element < Ap[row + 1] - 1) && (Ai[element] > Ai[element + 1])) + long int cursor = Ap[row]; + for (size_type block = 0; block < matrix.n_block_cols(); ++block) { - std::swap(Ai[element], Ai[element + 1]); - - std::swap(Ax[element], Ax[element + 1]); - if (numbers::NumberTraits::is_complex == true) - std::swap(Az[element], Az[element + 1]); - - ++element; + // find the next out-of-order element + while ((cursor < Ap[row + 1] - 1) && + (Ai[cursor] < Ai[cursor + 1])) + ++cursor; + + // if there is none, then just go on + if (cursor == Ap[row + 1] - 1) + break; + + // otherwise swap this entry with successive ones as long as + // necessary + long int element = cursor; + while ((element < Ap[row + 1] - 1) && + (Ai[element] > Ai[element + 1])) + { + std::swap(Ai[element], Ai[element + 1]); + + std::swap(Ax[element], Ax[element + 1]); + if (numbers::NumberTraits::is_complex == true) + std::swap(Az[element], Az[element + 1]); + + ++element; + } } } - } + }, + /* grain size = */ 50); } @@ -268,30 +276,31 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix) // then copy over matrix elements. note that for sparse matrices, // iterators are sorted so that they traverse each row from start to end - // before moving on to the next row. however, this isn't true for block - // matrices, so we have to do a bit of book keeping - - // loop over the elements of the matrix row by row, as suggested in the - // documentation of the sparse matrix iterator class - for (size_type row = 0; row < matrix.m(); ++row) - { - long int index = Ap[row]; - for (typename Matrix::const_iterator p = matrix.begin(row); - p != matrix.end(row); - ++p) + // before moving on to the next row. + parallel::apply_to_subranges( + 0, + matrix.m(), + [this, &matrix](const size_type row_begin, const size_type row_end) { + for (size_type row = row_begin; row < row_end; ++row) { - // write entry into the first free one for this row - Ai[index] = p->column(); - Ax[index] = std::real(p->value()); - if (numbers::NumberTraits::is_complex == true) - Az[index] = std::imag(p->value()); - - // then move pointer ahead - ++index; - } + long int index = Ap[row]; + for (typename Matrix::const_iterator p = matrix.begin(row); + p != matrix.end(row); + ++p) + { + // write entry into the first free one for this row + Ai[index] = p->column(); + Ax[index] = std::real(p->value()); + if (numbers::NumberTraits::is_complex == true) + Az[index] = std::imag(p->value()); - Assert(index == Ap[row + 1], ExcInternalError()); - } + // then move pointer ahead + ++index; + } + Assert(index == Ap[row + 1], ExcInternalError()); + } + }, + /* grain size = */ 50); // make sure that the elements in each row are sorted. we have to be more // careful for block sparse matrices, so ship this task out to a