// As already mentioned in the introduction, the system matrix is neither
// symmetric nor definite, and so it is not quite obvious how to come up
// with an iterative solver and a preconditioner that do a good job on this
- // matrix. We chose instead to go a different way and solve the linear
+ // matrix. (For more on this topic, see also the
+ // <a href="#extensions">Possibilities for extensions</a> section below.)
+ // We chose instead to go a different way and solve the linear
// system with the sparse LU decomposition provided by UMFPACK. This is
// often a good first choice for 2d problems and works reasonably well even
- // for a large number of DoFs. The deal.II interface to UMFPACK is given by
- // the SparseDirectUMFPACK class, which is very easy to use and allows us to
- // solve our linear system with just 3 lines of code.
-
+ // for moderately large numbers of DoFs. The deal.II interface to UMFPACK
+ // is implemented in the SparseDirectUMFPACK class, which is very easy to
+ // use and allows us to solve our linear system with just 3 lines of code.
+ //
// Note again that for compiling this example program, you need to have the
// deal.II library built with UMFPACK support.
template <int dim>
Timer timer;
// The code to solve the linear system is short: First, we allocate an
- // object of the right type. The following <code>initialize</code> call
- // provides the matrix that we would like to invert to the
- // SparseDirectUMFPACK object, and at the same time kicks off the
- // LU-decomposition. Hence, this is also the point where most of the
- // computational work in this program happens.
+ // object of the right type. The following call to
+ // SparseDirectUMFPACK::solve() takes as argument the matrix to decompose,
+ // and a vector that upon input equals the right hand side of the linear
+ // system to be solved, and upon output contains the solution of the linear
+ // system. To satisfy this input/output requirement, we first assign the
+ // right hand side vector to the `solution` variable.
SparseDirectUMFPACK A_direct;
- A_direct.initialize(system_matrix);
- // After the decomposition, we can use <code>A_direct</code> like a matrix
- // representing the inverse of our system matrix, so to compute the
- // solution we just have to multiply with the right hand side vector:
- A_direct.vmult(solution, system_rhs);
+ solution = system_rhs;
+ A_direct.solve(system_matrix, solution);
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;