void
solve(Vector<double> &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<std::complex<double>> &rhs_and_solution,
+ const bool transpose = false) const;
+
/**
* Same as before, but for block vectors.
*/
solve(BlockVector<double> &rhs_and_solution,
const bool transpose = false) const;
+ /**
+ * Same as before, but for complex-valued block vectors.
+ */
+ void
+ solve(BlockVector<std::complex<double>> &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.
Vector<double> &rhs_and_solution,
const bool transpose = false);
+ /**
+ * Same as before, but for complex-valued solution vectors.
+ */
+ template <class Matrix>
+ void
+ solve(const Matrix & matrix,
+ Vector<std::complex<double>> &rhs_and_solution,
+ const bool transpose = false);
+
/**
* Same as before, but for block vectors.
*/
BlockVector<double> &rhs_and_solution,
const bool transpose = false);
+ /**
+ * Same as before, but for complex-valued block vectors.
+ */
+ template <class Matrix>
+ void
+ solve(const Matrix & matrix,
+ BlockVector<std::complex<double>> &rhs_and_solution,
+ const bool transpose = false);
+
/**
* @}
*/
sort_arrays(const BlockSparseMatrix<number> &);
/**
- * 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<types::suitesparse_index> Ap;
std::vector<types::suitesparse_index> Ai;
std::vector<double> Ax;
+ std::vector<double> Az;
/**
* Control and work arrays for the solver routines.
// ---------------------------------------------------------------------
#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/numbers.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/vector.h>
#include <cerrno>
+#include <complex>
#include <iostream>
#include <list>
#include <typeinfo>
tmp.swap(Ax);
}
+ {
+ std::vector<double> tmp;
+ tmp.swap(Az);
+ }
+
umfpack_dl_defaults(control.data());
}
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<number>::is_complex == true)
+ std::swap(Az[cursor], Az[cursor + 1]);
+
++cursor;
}
}
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<number>::is_complex == true)
+ std::swap(Az[cursor], Az[cursor + 1]);
+
++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<number>::is_complex == true)
+ std::swap(Az[cursor], Az[cursor + 1]);
+
++element;
}
}
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();
Ap.resize(N + 1);
Ai.resize(matrix.n_nonzero_elements());
Ax.resize(matrix.n_nonzero_elements());
+ if (numbers::NumberTraits<number>::is_complex == true)
+ Az.resize(matrix.n_nonzero_elements());
// first fill row lengths array
Ap[0] = 0;
{
// 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<number>::is_complex == true)
+ Az[row_pointers[row]] = std::imag(p->value());
// then move pointer ahead
++row_pointers[row];
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<number>::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<number>::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));
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<double> rhs(rhs_and_solution.size());
rhs = rhs_and_solution;
}
+
+void
+SparseDirectUMFPACK::solve(Vector<std::complex<double>> &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<double> rhs_re(rhs_and_solution.size());
+ Vector<double> 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<double> solution_re(rhs_and_solution.size());
+ Vector<double> 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<std::complex<double>> 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<double> 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<double> &rhs_and_solution,
const bool transpose /*=false*/) const
+void
+SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &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<std::complex<double>> tmp(rhs_and_solution.size());
+ tmp = rhs_and_solution;
+ solve(tmp, transpose);
+ rhs_and_solution = tmp;
+}
+
+
+
template <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix & matrix,
}
+
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix & matrix,
+ Vector<std::complex<double>> &rhs_and_solution,
+ const bool transpose /*=false*/)
+{
+ factorize(matrix);
+ solve(rhs_and_solution, transpose);
+}
+
+
+
template <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix & matrix,
}
+
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix & matrix,
+ BlockVector<std::complex<double>> &rhs_and_solution,
+ const bool transpose /*=false*/)
+{
+ factorize(matrix);
+ solve(rhs_and_solution, transpose);
+}
+
+
#else
+void
+SparseDirectUMFPACK::solve(Vector<std::complex<double>> &, 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<double> &, const bool) const
{
}
+
+void
+SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &,
+ 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 <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, const bool)
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix &,
+ Vector<std::complex<double>> &,
+ 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 <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, 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 <class Matrix>
+void
+SparseDirectUMFPACK::solve(const Matrix &,
+ BlockVector<std::complex<double>> &,
+ 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
// explicit instantiations for SparseMatrixUMFPACK
-#define InstantiateUMFPACK(MatrixType) \
- template void SparseDirectUMFPACK::factorize(const MatrixType &); \
- template void SparseDirectUMFPACK::solve(const MatrixType &, \
- Vector<double> &, \
- bool); \
- template void SparseDirectUMFPACK::solve(const MatrixType &, \
- BlockVector<double> &, \
- bool); \
- template void SparseDirectUMFPACK::initialize(const MatrixType &, \
+#define InstantiateUMFPACK(MatrixType) \
+ template void SparseDirectUMFPACK::factorize(const MatrixType &); \
+ template void SparseDirectUMFPACK::solve(const MatrixType &, \
+ Vector<double> &, \
+ const bool); \
+ template void SparseDirectUMFPACK::solve(const MatrixType &, \
+ Vector<std::complex<double>> &, \
+ const bool); \
+ template void SparseDirectUMFPACK::solve(const MatrixType &, \
+ BlockVector<double> &, \
+ const bool); \
+ template void SparseDirectUMFPACK::solve( \
+ const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
+ template void SparseDirectUMFPACK::initialize(const MatrixType &, \
const AdditionalData)
+// Instantiate everything for real-valued matrices
InstantiateUMFPACK(SparseMatrix<double>);
InstantiateUMFPACK(SparseMatrix<float>);
InstantiateUMFPACK(SparseMatrixEZ<double>);
InstantiateUMFPACK(BlockSparseMatrix<double>);
InstantiateUMFPACK(BlockSparseMatrix<float>);
+// Now also for complex-valued matrices
+InstantiateUMFPACK(SparseMatrix<std::complex<double>>);
+InstantiateUMFPACK(SparseMatrix<std::complex<float>>);
+InstantiateUMFPACK(BlockSparseMatrix<std::complex<double>>);
+InstantiateUMFPACK(BlockSparseMatrix<std::complex<float>>);
+
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff