From a4f43a69caf8681f3c6325735b3c1d19775ed84a Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 25 Oct 2006 00:01:44 +0000 Subject: [PATCH] Revamp the functionality of sorting each row for block matrices git-svn-id: https://svn.dealii.org/trunk@14075 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_direct.h | 18 +++ deal.II/lac/source/sparse_direct.cc | 156 +++++++++++++++++------- 2 files changed, 128 insertions(+), 46 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h index ddac248a23..685bdcfe8e 100644 --- a/deal.II/lac/include/lac/sparse_direct.h +++ b/deal.II/lac/include/lac/sparse_direct.h @@ -21,6 +21,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -1212,6 +1213,23 @@ class SparseDirectUMFPACK : public Subscriptor */ void clear (); + + /** + * Make sure that the arrays Ai + * and Ap are sorted in each + * row. UMFPACK wants it this + * way. We need to have two + * versions of this function, one + * for the usual SparseMatrix and + * one for the BlockSparseMatrix + * classes + */ + template + void sort_arrays (const SparseMatrix &); + + template + void sort_arrays (const BlockSparseMatrix &); + /** * The arrays in which we store the data * for the solver. diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc index 8f2bd38e3d..8c7a2c7bcf 100644 --- a/deal.II/lac/source/sparse_direct.cc +++ b/deal.II/lac/source/sparse_direct.cc @@ -1682,6 +1682,106 @@ SparseDirectUMFPACK::clear () +template +void +SparseDirectUMFPACK:: +sort_arrays (const SparseMatrix &) +{ + // do the copying around of entries + // so that the diagonal entry is in the + // right place. note that this is easy to + // detect: since all entries apart from the + // diagonal entry are sorted, we know that + // the diagonal entry is in the wrong place + // if and only if its column index is + // larger than the column index of the + // second entry in a row + // + // ignore rows with only one or no entry + for (unsigned int row=0; row Ai[cursor+1])) + { + std::swap (Ai[cursor], Ai[cursor+1]); + std::swap (Ax[cursor], Ax[cursor+1]); + ++cursor; + } + } +} + + + +template +void +SparseDirectUMFPACK:: +sort_arrays (const BlockSparseMatrix &matrix) +{ + // the case for block matrices is a + // bit more difficult, since all we + // know is that *within each + // block*, the diagonal of that + // block may come 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 (unsigned int row=0; row Ai[element+1])) + { + std::swap (Ai[element], Ai[element+1]); + std::swap (Ax[element], Ax[element+1]); + ++element; + } + } + } +} + + + template void SparseDirectUMFPACK:: @@ -1766,63 +1866,27 @@ factorize (const Matrix &matrix) Assert (row_pointers[i] == Ap[i+1], ExcInternalError()); } - // finally do the copying around of entries - // so that the diagonal entry is in the - // right place. note that this is easy to - // detect: since all entries apart from the - // diagonal entry are sorted, we know that - // the diagonal entry is in the wrong place - // if and only if its column index is - // larger than the column index of the - // second entry in a row - // - // ignore rows with only one or no entry - { - for (unsigned int row=0; row Ai[cursor+1])) - { - std::swap (Ai[cursor], Ai[cursor+1]); - std::swap (Ax[cursor], Ax[cursor+1]); - ++cursor; - } - } - } + // 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 different function + sort_array (); int status; - status = umfpack_di_symbolic (N, N, &Ap[0], &Ai[0], &Ax[0], &symbolic_decomposition, &control[0], 0); - AssertThrow (status == UMFPACK_OK, ExcUMFPACKError("umfpack_di_symbolic", status)); + AssertThrow (status == UMFPACK_OK, + ExcUMFPACKError("umfpack_di_symbolic", status)); status = umfpack_di_numeric (&Ap[0], &Ai[0], &Ax[0], symbolic_decomposition, &numeric_decomposition, &control[0], 0); - AssertThrow (status == UMFPACK_OK, ExcUMFPACKError("umfpack_di_numeric", status)); + AssertThrow (status == UMFPACK_OK, + ExcUMFPACKError("umfpack_di_numeric", status)); umfpack_di_free_symbolic (&symbolic_decomposition) ; } -- 2.39.5