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

## Advertising

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 :-)) Best, Matthias -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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.

**
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.depth_console(2); 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; sparsity_pattern.copy_from(dsp); 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); }