Hi all,

I discussed the preconditioned CG method today in class (with SSOR and
optimal relaxation parameter applied to the beloved Laplace problem.)

I was quite surprised that I was able to solve the problem for 1 million
degrees of freedom in about 9 seconds and 80 CG iterations with Julia
(that I use in class). Granted the tolerance for the convergence check
wasn't set very high...

    julia> run_test()
    8.372644 seconds (1.97 k allocations: 7.065 GiB, 40.45% gc time)
    Number of iterations: 80
    Error L2: 0.009811802487601433

But bottom line is, I checked how fast I can solve the very same problem
with deal.II's builtin linear algebra:

    $ time ./step
    DEAL::Number of degrees of freedom: 1048576
    DEAL:cg::Starting value 1024.00
    DEAL:cg::Convergence step 81 value 0.00899413
    ./step  8.15s user 0.14s system 210% cpu 3.936 total

(Startup and initalization and assembly - without the solve - takes
roughly 0.9s - So we are more at 3.0s for the solve step. I was simply
too lazy to include a timer.)

That's about as optimal as it gets I'd say.

Code is attached if someone feels the urge to improve above results (or
maybe quickly benchmark the other two linear algebra interfaces we have


The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Attachment: step.jl
Description: Binary data

#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>

using namespace dealii;

int main()
  constexpr unsigned int m = 1024;
  constexpr unsigned int n = m * m;

  deallog << "Number of degrees of freedom: " << n << std::endl;

  DynamicSparsityPattern dsp(n);
  for (unsigned int i = 0; i < n; ++i) {
    dsp.add(i, i);
    if (i < n - 1) {
      dsp.add(i, i + 1);
      dsp.add(i + 1, i);
    if (i < n - m - 1) {
      dsp.add(i, i + m);
      dsp.add(i + m, i);
  SparsityPattern sparsity_pattern;

  SparseMatrix<double> system_matrix(sparsity_pattern);
  for (unsigned int i = 0; i < n; ++i) {
    system_matrix.set(i, i, 4.);
    if (i < n - 1) {
      system_matrix.set(i, i + 1, -1.);
      system_matrix.set(i + 1, i, -1.);
    if (i < n - m - 1) {
      system_matrix.set(i, i + m, -1.);
      system_matrix.set(i + m, i, -1.);

  Vector<double> system_rhs(n);
  for (auto &it : system_rhs)
    it = 1.0;

  Vector<double> solution(n);

  SolverControl solver_control(10000, 1e-2);
  SolverCG<> solver(solver_control);

  constexpr double optimal_omega = 1.9938888033081086;
  PreconditionSSOR<> preconditioner;
  preconditioner.initialize(system_matrix, optimal_omega);

  solver.solve(system_matrix, solution, system_rhs, preconditioner);

Reply via email to