From: Wolfgang Bangerth Date: Mon, 30 Dec 2019 19:07:39 +0000 (-0700) Subject: Implement the solution of complex linear systems using UMFPACK. X-Git-Tag: v9.2.0-rc1~612^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7d11989c51153c4d7e3e2a445500b49624fcd58f;p=dealii.git Implement the solution of complex linear systems using UMFPACK. --- diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index a354636531..d4cd647a06 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -244,6 +244,21 @@ public: void solve(Vector &rhs_and_solution, const bool transpose = false) const; + /** + * Like the previous function, but for a complex-valued right hand side + * and solution vector. + * + * If the matrix that was previously factorized had complex-valued entries, + * then the `rhs_and_solution` vector will, upon return from this function, + * simply contain the solution of the linear system $Ax=b$. If the matrix + * was real-valued, then this is also true, but the solution will simply + * be computed by applying the factorized $A^{-1}$ to both the + * real and imaginary parts of the right hand side vector. + */ + void + solve(Vector> &rhs_and_solution, + const bool transpose = false) const; + /** * Same as before, but for block vectors. */ @@ -251,6 +266,13 @@ public: solve(BlockVector &rhs_and_solution, const bool transpose = false) const; + /** + * Same as before, but for complex-valued block vectors. + */ + void + solve(BlockVector> &rhs_and_solution, + const bool transpose = false) const; + /** * Call the two functions factorize() and solve() in that order, i.e. * perform the whole solution process for the given right hand side vector. @@ -263,6 +285,15 @@ public: Vector &rhs_and_solution, const bool transpose = false); + /** + * Same as before, but for complex-valued solution vectors. + */ + template + void + solve(const Matrix & matrix, + Vector> &rhs_and_solution, + const bool transpose = false); + /** * Same as before, but for block vectors. */ @@ -272,6 +303,15 @@ public: BlockVector &rhs_and_solution, const bool transpose = false); + /** + * Same as before, but for complex-valued block vectors. + */ + template + void + solve(const Matrix & matrix, + BlockVector> &rhs_and_solution, + const bool transpose = false); + /** * @} */ @@ -362,11 +402,20 @@ private: sort_arrays(const BlockSparseMatrix &); /** - * The arrays in which we store the data for the solver. + * The arrays in which we store the data for the solver. These are documented + * in the descriptions of the umfpack_*_symbolic() and umfpack_*_numeric() + * functions, but in short: + * - `Ap` is the array saying which row starts where in `Ai` + * - `Ai` is the array that stores the column indices of nonzero entries + * - `Ax` is the array that stores the values of nonzero entries; if the + * matrix is complex-valued, then it stores the real parts + * - `Az` is the array that stores the imaginary parts of nonzero entries, + * and is used only if the matrix is complex-valued. */ std::vector Ap; std::vector Ai; std::vector Ax; + std::vector Az; /** * Control and work arrays for the solver routines. diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index 9a225192ab..b158a65616 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -14,6 +14,7 @@ // --------------------------------------------------------------------- #include +#include #include #include @@ -21,6 +22,7 @@ #include #include +#include #include #include #include @@ -93,6 +95,11 @@ SparseDirectUMFPACK::clear() tmp.swap(Ax); } + { + std::vector tmp; + tmp.swap(Az); + } + umfpack_dl_defaults(control.data()); } @@ -126,7 +133,11 @@ SparseDirectUMFPACK::sort_arrays(const SparseMatrix &matrix) while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1])) { std::swap(Ai[cursor], Ai[cursor + 1]); + std::swap(Ax[cursor], Ax[cursor + 1]); + if (numbers::NumberTraits::is_complex == true) + std::swap(Az[cursor], Az[cursor + 1]); + ++cursor; } } @@ -145,7 +156,11 @@ SparseDirectUMFPACK::sort_arrays(const SparseMatrixEZ &matrix) while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1])) { std::swap(Ai[cursor], Ai[cursor + 1]); + std::swap(Ax[cursor], Ax[cursor + 1]); + if (numbers::NumberTraits::is_complex == true) + std::swap(Az[cursor], Az[cursor + 1]); + ++cursor; } } @@ -181,7 +196,11 @@ SparseDirectUMFPACK::sort_arrays(const BlockSparseMatrix &matrix) 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[cursor], Az[cursor + 1]); + ++element; } } @@ -194,9 +213,11 @@ template void SparseDirectUMFPACK::factorize(const Matrix &matrix) { - Assert(matrix.m() == matrix.n(), ExcNotQuadratic()) + Assert(matrix.m() == matrix.n(), ExcNotQuadratic()); - clear(); + clear(); + + using number = typename Matrix::value_type; n_rows = matrix.m(); n_cols = matrix.n(); @@ -224,6 +245,8 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix) Ap.resize(N + 1); Ai.resize(matrix.n_nonzero_elements()); Ax.resize(matrix.n_nonzero_elements()); + if (numbers::NumberTraits::is_complex == true) + Az.resize(matrix.n_nonzero_elements()); // first fill row lengths array Ap[0] = 0; @@ -250,7 +273,9 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix) { // write entry into the first free one for this row Ai[row_pointers[row]] = p->column(); - Ax[row_pointers[row]] = p->value(); + Ax[row_pointers[row]] = std::real(p->value()); + if (numbers::NumberTraits::is_complex == true) + Az[row_pointers[row]] = std::imag(p->value()); // then move pointer ahead ++row_pointers[row]; @@ -268,24 +293,45 @@ SparseDirectUMFPACK::factorize(const Matrix &matrix) sort_arrays(matrix); int status; - status = umfpack_dl_symbolic(N, - N, - Ap.data(), - Ai.data(), - Ax.data(), - &symbolic_decomposition, - control.data(), - nullptr); + if (numbers::NumberTraits::is_complex == false) + status = umfpack_dl_symbolic(N, + N, + Ap.data(), + Ai.data(), + Ax.data(), + &symbolic_decomposition, + control.data(), + nullptr); + else + status = umfpack_zl_symbolic(N, + N, + Ap.data(), + Ai.data(), + Ax.data(), + Az.data(), + &symbolic_decomposition, + control.data(), + nullptr); AssertThrow(status == UMFPACK_OK, ExcUMFPACKError("umfpack_dl_symbolic", status)); - status = umfpack_dl_numeric(Ap.data(), - Ai.data(), - Ax.data(), - symbolic_decomposition, - &numeric_decomposition, - control.data(), - nullptr); + if (numbers::NumberTraits::is_complex == false) + status = umfpack_dl_numeric(Ap.data(), + Ai.data(), + Ax.data(), + symbolic_decomposition, + &numeric_decomposition, + control.data(), + nullptr); + else + status = umfpack_zl_numeric(Ap.data(), + Ai.data(), + Ax.data(), + Az.data(), + symbolic_decomposition, + &numeric_decomposition, + control.data(), + nullptr); AssertThrow(status == UMFPACK_OK, ExcUMFPACKError("umfpack_dl_numeric", status)); @@ -303,6 +349,13 @@ SparseDirectUMFPACK::solve(Vector &rhs_and_solution, Assert(Ai.size() != 0, ExcNotInitialized()); Assert(Ai.size() == Ax.size(), ExcNotInitialized()); + Assert(Az.size() == 0, + ExcMessage("You have previously factored a matrix using this class " + "that had complex-valued entries. This then requires " + "applying the factored matrix to a complex-valued " + "vector, but you are only providing a real-valued vector " + "here.")); + Vector rhs(rhs_and_solution.size()); rhs = rhs_and_solution; @@ -326,6 +379,111 @@ SparseDirectUMFPACK::solve(Vector &rhs_and_solution, } + +void +SparseDirectUMFPACK::solve(Vector> &rhs_and_solution, + const bool transpose /*=false*/) const +{ + // make sure that some kind of factorize() call has happened before + Assert(Ap.size() != 0, ExcNotInitialized()); + Assert(Ai.size() != 0, ExcNotInitialized()); + Assert(Ai.size() == Ax.size(), ExcNotInitialized()); + + // First see whether the matrix that was factorized was complex-valued. + // If so, just apply the complex factorization to the vector. + if (Az.size() != 0) + { + Assert(Ax.size() == Az.size(), ExcInternalError()); + + // It would be nice if we could just present a pointer to the + // first element of the complex-valued solution vector and let + // UMFPACK fill both the real and imaginary parts of the solution + // vector at that address. UMFPACK calls this the 'packed' format, + // and in those cases it only takes one pointer to the entire + // vector, rather than a pointer to the real and one pointer to + // the imaginary parts of the vector. The problem is that if we + // want to pack, then we also need to pack the matrix, and the + // functions above have already decided that we don't want to pack + // the matrix but instead deal with split format for the matrix, + // and then UMFPACK decides that it can't deal with a split matrix + // and a packed vector. We have to choose one or the other, not + // mix. + // + // So create four vectors, one each for the real and imaginary parts + // of the right hand side and solution. + Vector rhs_re(rhs_and_solution.size()); + Vector rhs_im(rhs_and_solution.size()); + for (unsigned int i = 0; i < rhs_and_solution.size(); ++i) + { + rhs_re(i) = std::real(rhs_and_solution(i)); + rhs_im(i) = std::imag(rhs_and_solution(i)); + } + + Vector solution_re(rhs_and_solution.size()); + Vector solution_im(rhs_and_solution.size()); + + // Solve the system. note that since UMFPACK wants compressed column + // storage instead of the compressed row storage format we use in + // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead + // + // Conversely, if we solve for the transpose, we have to use UMFPACK_A + // instead. + // + // Note that for the complex case here, the transpose is selected using + // UMFPACK_Aat, not UMFPACK_At. + const int status = umfpack_zl_solve(transpose ? UMFPACK_A : UMFPACK_Aat, + Ap.data(), + Ai.data(), + Ax.data(), + Az.data(), + solution_re.data(), + solution_im.data(), + rhs_re.data(), + rhs_im.data(), + numeric_decomposition, + control.data(), + nullptr); + AssertThrow(status == UMFPACK_OK, + ExcUMFPACKError("umfpack_dl_solve", status)); + + // Now put things back together into the output vector + for (unsigned int i = 0; i < rhs_and_solution.size(); ++i) + rhs_and_solution(i) = {solution_re(i), solution_im(i)}; + } + else + { + // We have factorized a real-valued matrix, but the rhs and solution + // vectors are complex-valued. UMFPACK does not natively support this + // case, but we can just apply the factorization to real and imaginary + // parts of the right hand side separately + const Vector> rhs = rhs_and_solution; + + // Get the real part of the right hand side, solve with it, and copy the + // results into the result vector by just copying the real output + // into the complex-valued result vector (which implies setting the + // imaginary parts to zero): + Vector rhs_real_or_imag(rhs_and_solution.size()); + for (unsigned int i = 0; i < rhs.size(); ++i) + rhs_real_or_imag(i) = std::real(rhs(i)); + + solve(rhs_real_or_imag, transpose); + + rhs_and_solution = rhs_real_or_imag; + + // Then repeat the whole thing with the imaginary part. The copying step + // is more complicated now because we can only touch the imaginary + // component of the output vector: + for (unsigned int i = 0; i < rhs.size(); ++i) + rhs_real_or_imag(i) = std::imag(rhs(i)); + + solve(rhs_real_or_imag, transpose); + + for (unsigned int i = 0; i < rhs.size(); ++i) + rhs_and_solution(i).imag(rhs_real_or_imag(i)); + } +} + + void SparseDirectUMFPACK::solve(BlockVector &rhs_and_solution, const bool transpose /*=false*/) const @@ -341,6 +499,21 @@ SparseDirectUMFPACK::solve(BlockVector &rhs_and_solution, +void +SparseDirectUMFPACK::solve(BlockVector> &rhs_and_solution, + const bool transpose /*=false*/) const +{ + // the UMFPACK functions want a contiguous array of elements, so + // there is no way around copying data around. thus, just copy the + // data into a regular vector and back + Vector> tmp(rhs_and_solution.size()); + tmp = rhs_and_solution; + solve(tmp, transpose); + rhs_and_solution = tmp; +} + + + template void SparseDirectUMFPACK::solve(const Matrix & matrix, @@ -352,6 +525,19 @@ SparseDirectUMFPACK::solve(const Matrix & matrix, } + +template +void +SparseDirectUMFPACK::solve(const Matrix & matrix, + Vector> &rhs_and_solution, + const bool transpose /*=false*/) +{ + factorize(matrix); + solve(rhs_and_solution, transpose); +} + + + template void SparseDirectUMFPACK::solve(const Matrix & matrix, @@ -363,6 +549,18 @@ SparseDirectUMFPACK::solve(const Matrix & matrix, } + +template +void +SparseDirectUMFPACK::solve(const Matrix & matrix, + BlockVector> &rhs_and_solution, + const bool transpose /*=false*/) +{ + factorize(matrix); + solve(rhs_and_solution, transpose); +} + + #else @@ -402,6 +600,17 @@ SparseDirectUMFPACK::solve(Vector &, const bool) const +void +SparseDirectUMFPACK::solve(Vector> &, const bool) const +{ + AssertThrow( + false, + ExcMessage( + "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); +} + + + void SparseDirectUMFPACK::solve(BlockVector &, const bool) const { @@ -412,6 +621,19 @@ SparseDirectUMFPACK::solve(BlockVector &, const bool) const } + +void +SparseDirectUMFPACK::solve(BlockVector> &, + const bool) const +{ + AssertThrow( + false, + ExcMessage( + "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); +} + + + template void SparseDirectUMFPACK::solve(const Matrix &, Vector &, const bool) @@ -424,6 +646,20 @@ SparseDirectUMFPACK::solve(const Matrix &, Vector &, const bool) +template +void +SparseDirectUMFPACK::solve(const Matrix &, + Vector> &, + const bool) +{ + AssertThrow( + false, + ExcMessage( + "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); +} + + + template void SparseDirectUMFPACK::solve(const Matrix &, BlockVector &, const bool) @@ -434,6 +670,20 @@ SparseDirectUMFPACK::solve(const Matrix &, BlockVector &, const bool) "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); } + + +template +void +SparseDirectUMFPACK::solve(const Matrix &, + BlockVector> &, + const bool) +{ + AssertThrow( + false, + ExcMessage( + "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); +} + #endif @@ -497,17 +747,23 @@ SparseDirectUMFPACK::n() const // explicit instantiations for SparseMatrixUMFPACK -#define InstantiateUMFPACK(MatrixType) \ - template void SparseDirectUMFPACK::factorize(const MatrixType &); \ - template void SparseDirectUMFPACK::solve(const MatrixType &, \ - Vector &, \ - bool); \ - template void SparseDirectUMFPACK::solve(const MatrixType &, \ - BlockVector &, \ - bool); \ - template void SparseDirectUMFPACK::initialize(const MatrixType &, \ +#define InstantiateUMFPACK(MatrixType) \ + template void SparseDirectUMFPACK::factorize(const MatrixType &); \ + template void SparseDirectUMFPACK::solve(const MatrixType &, \ + Vector &, \ + const bool); \ + template void SparseDirectUMFPACK::solve(const MatrixType &, \ + Vector> &, \ + const bool); \ + template void SparseDirectUMFPACK::solve(const MatrixType &, \ + BlockVector &, \ + const bool); \ + template void SparseDirectUMFPACK::solve( \ + const MatrixType &, BlockVector> &, const bool); \ + template void SparseDirectUMFPACK::initialize(const MatrixType &, \ const AdditionalData) +// Instantiate everything for real-valued matrices InstantiateUMFPACK(SparseMatrix); InstantiateUMFPACK(SparseMatrix); InstantiateUMFPACK(SparseMatrixEZ); @@ -515,5 +771,11 @@ InstantiateUMFPACK(SparseMatrixEZ); InstantiateUMFPACK(BlockSparseMatrix); InstantiateUMFPACK(BlockSparseMatrix); +// Now also for complex-valued matrices +InstantiateUMFPACK(SparseMatrix>); +InstantiateUMFPACK(SparseMatrix>); +InstantiateUMFPACK(BlockSparseMatrix>); +InstantiateUMFPACK(BlockSparseMatrix>); + DEAL_II_NAMESPACE_CLOSE