}
+%-------------------------------------------------------------------------------
+% Step 93
+%-------------------------------------------------------------------------------
+
+@inbook{Battermann_1998,
+ title={Preconditioners for Karush-Kuhn-Tucker Matrices Arising in the Optimal Control of Distributed Systems},
+ url={http://dx.doi.org/10.1007/978-3-0348-8849-3_2},
+ DOI={10.1007/978-3-0348-8849-3_2},
+ booktitle={Control and Estimation of Distributed Parameter Systems},
+ publisher={Birkh{\"a}user Basel},
+ author={Battermann, A. and Heinkenschloss, M.},
+ year={1998},
+ pages={15-32}
+}
+
+
%-------------------------------------------------------------------------------
% References used elsewhere
%-------------------------------------------------------------------------------
0 & \mathcal{F} & 0
\end{array}\right)
\left(\begin{array}{c}
- U\\\Lambda\\ C
+ U\\ \Lambda \\ C
\end{array}\right) =
\left(\begin{array}{c}
\overline{\mathcal{U}}\\0\\0
diagonal, preconditioners like PreconditionJacobi will not work
because they divide by diagonal entries. Instead, block
preconditioners such as those discussed in step-20 or step-22 (among
-many others) will likely be useful.
+many others) will likely be useful. For block preconditioners, the
+key realizations is that the blocks of the system matrix
+@f{align*}{
+ \left(\begin{array}{c c c}
+ \mathcal{M} & -\mathcal{N}^T & 0\\
+ -\mathcal{N} & 0 & \mathcal{F}^T\\
+ 0 & \mathcal{F} & 0
+ \end{array}\right)
+@f}
+can often individually be solved with quite efficiently; for example,
+$\mathcal{M}$ is a mass matrix that is easily solved with using a CG
+iteration, and $\mathcal N$ is a Laplace matrix for which CG with
+a geometric or algebraic multigrid preconditioner is very effective.
+Using the ideas of step-20 and step-22, we should then create a
+$3\times 3$ block preconditioner in which some blocks correspond to the
+inverses of $\mathcal{M}$ or $\mathcal{N}$, or some kind of Schur
+complement. A starting point for this kind of consideration is
+@cite Battermann_1998 .
2. To validate the optimization problem is working correctly, we could
try to match a target function which is itself a solution to the
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_minres.h>
-#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
// solver. For smaller problems, one can also use a direct solver
// (see step-29) for which you would just replace the main part of
// this function by the following three lines of code:
- //
- // `SparseDirectUMFPACK direct_solver;`
- //
- // `direct_solver.initialize(system_matrix);`
- //
- // `direct_solver.vmult(solution, system_rhs);`
+ // @code
+ // SparseDirectUMFPACK direct_solver;
+ // direct_solver.initialize(system_matrix);
+ // direct_solver.vmult(solution, system_rhs);
+ // @endcode
template <int dim>
void Step93<dim>::solve()
{
std::cout << "Beginning solve..." << std::endl;
- Timer timer;
+ {
+ Timer timer;
- SolverControl solver_control(1'000'000, 1e-6 * system_rhs.l2_norm());
- SolverMinRes<Vector<double>> solver(solver_control);
+ SolverControl solver_control(1'000'000, 1e-6 * system_rhs.l2_norm());
+ SolverMinRes<Vector<double>> solver(solver_control);
- solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+ solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
- timer.stop();
+ timer.stop();
+ }
std::cout << "Wall time: " << timer.wall_time() << "s" << std::endl;
std::cout << "Solved in " << solver_control.last_step()
<< " MINRES iterations." << std::endl;