deallog.depth_console(0);
deallog.threshold_double(1.e-10);
- SolverControl control(100, 1.e-3);
+ SolverControl control(10, 1.e-3);
SolverRichardson<> rich(control);
SolverRelaxation<> relax(control);
DEAL::Size 33 Unknowns 1024
DEAL::Block size 2
DEAL:Jacobi:Richardson::Starting value 32.00
-DEAL:Jacobi:Richardson::Failure step 100 value 16.51
+DEAL:Jacobi:Richardson::Failure step 10 value 27.72
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 16.51
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.72
DEAL:Jacobi::diff 0
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 16.51
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.72
DEAL:Jacobi::diff 0
DEAL:SOR :Richardson::Starting value 32.00
-DEAL:SOR :Richardson::Failure step 100 value 12.27
+DEAL:SOR :Richardson::Failure step 10 value 26.93
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 12.27
+DEAL:SOR :Relaxation::Failure step 10 value 26.93
DEAL:SOR ::diff 0
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 12.27
+DEAL:SOR :Relaxation::Failure step 10 value 26.93
DEAL:SOR ::diff 0
DEAL:SSOR :Richardson::Starting value 32.00
-DEAL:SSOR :Richardson::Failure step 100 value 5.660
+DEAL:SSOR :Richardson::Failure step 10 value 23.78
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 5.660
+DEAL:SSOR :Relaxation::Failure step 10 value 23.78
DEAL:SSOR ::diff 0
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 5.660
+DEAL:SSOR :Relaxation::Failure step 10 value 23.78
DEAL:SSOR ::diff 0
DEAL::Block size 4
DEAL:Jacobi:Richardson::Starting value 32.00
-DEAL:Jacobi:Richardson::Failure step 100 value 15.27
+DEAL:Jacobi:Richardson::Failure step 10 value 27.72
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 15.27
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.72
DEAL:Jacobi::diff 0
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 15.27
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.72
DEAL:Jacobi::diff 0
DEAL:SOR :Richardson::Starting value 32.00
-DEAL:SOR :Richardson::Failure step 100 value 10.60
+DEAL:SOR :Richardson::Failure step 10 value 26.59
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 10.60
+DEAL:SOR :Relaxation::Failure step 10 value 26.59
DEAL:SOR ::diff 0
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 10.60
+DEAL:SOR :Relaxation::Failure step 10 value 26.59
DEAL:SOR ::diff 0
DEAL:SSOR :Richardson::Starting value 32.00
-DEAL:SSOR :Richardson::Failure step 100 value 4.196
+DEAL:SSOR :Richardson::Failure step 10 value 23.33
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 4.196
+DEAL:SSOR :Relaxation::Failure step 10 value 23.33
DEAL:SSOR ::diff 0
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 4.196
+DEAL:SSOR :Relaxation::Failure step 10 value 23.33
DEAL:SSOR ::diff 0
DEAL::Block size 8
DEAL:Jacobi:Richardson::Starting value 32.00
-DEAL:Jacobi:Richardson::Failure step 100 value 14.29
+DEAL:Jacobi:Richardson::Failure step 10 value 27.14
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 14.29
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.14
DEAL:Jacobi::diff 0
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 14.29
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.14
DEAL:Jacobi::diff 0
DEAL:SOR :Richardson::Starting value 32.00
-DEAL:SOR :Richardson::Failure step 100 value 9.455
+DEAL:SOR :Richardson::Failure step 10 value 25.81
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 9.455
+DEAL:SOR :Relaxation::Failure step 10 value 25.81
DEAL:SOR ::diff 0
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 9.455
+DEAL:SOR :Relaxation::Failure step 10 value 25.81
DEAL:SOR ::diff 0
DEAL:SSOR :Richardson::Starting value 32.00
-DEAL:SSOR :Richardson::Failure step 100 value 3.335
+DEAL:SSOR :Richardson::Failure step 10 value 22.60
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 3.335
+DEAL:SSOR :Relaxation::Failure step 10 value 22.60
DEAL:SSOR ::diff 0
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 3.335
+DEAL:SSOR :Relaxation::Failure step 10 value 22.60
DEAL:SSOR ::diff 0
DEAL::Block size 16
DEAL:Jacobi:Richardson::Starting value 32.00
-DEAL:Jacobi:Richardson::Failure step 100 value 13.60
+DEAL:Jacobi:Richardson::Failure step 10 value 26.75
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 13.60
+DEAL:Jacobi:Relaxation::Failure step 10 value 26.75
DEAL:Jacobi::diff 0
DEAL:Jacobi:Relaxation::Starting value 32.00
-DEAL:Jacobi:Relaxation::Failure step 100 value 13.60
+DEAL:Jacobi:Relaxation::Failure step 10 value 26.75
DEAL:Jacobi::diff 0
DEAL:SOR :Richardson::Starting value 32.00
-DEAL:SOR :Richardson::Failure step 100 value 8.752
+DEAL:SOR :Richardson::Failure step 10 value 25.13
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 8.752
+DEAL:SOR :Relaxation::Failure step 10 value 25.13
DEAL:SOR ::diff 0
DEAL:SOR :Relaxation::Starting value 32.00
-DEAL:SOR :Relaxation::Failure step 100 value 8.752
+DEAL:SOR :Relaxation::Failure step 10 value 25.13
DEAL:SOR ::diff 0
DEAL:SSOR :Richardson::Starting value 32.00
-DEAL:SSOR :Richardson::Failure step 100 value 2.874
+DEAL:SSOR :Richardson::Failure step 10 value 21.81
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 2.874
+DEAL:SSOR :Relaxation::Failure step 10 value 21.81
DEAL:SSOR ::diff 0
DEAL:SSOR :Relaxation::Starting value 32.00
-DEAL:SSOR :Relaxation::Failure step 100 value 2.874
+DEAL:SSOR :Relaxation::Failure step 10 value 21.81
DEAL:SSOR ::diff 0
--- /dev/null
+//----------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+
+// Compare preconditioned Richardson with block relaxation in
+// different permutation modes. All output diffs
+// should be zero.
+
+#include "../tests.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/solver_relaxation.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/relaxation_block.h>
+
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+double
+check_solve( SOLVER& solver, const MATRIX& A,
+ VECTOR& u, VECTOR& f, const PRECONDITION& P)
+{
+ double result = 0.;
+ u = 0.;
+ f = 1.;
+ try
+ {
+ solver.solve(A,u,f,P);
+ }
+ catch (SolverControl::NoConvergence& e)
+ {
+ result = e.last_residual;
+ }
+ return result;
+}
+
+int main()
+{
+ const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ std::ofstream logfile(logname.c_str());
+// logfile.setf(std::ios::fixed);
+ deallog << std::setprecision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ SolverControl control(10, 1.e-3);
+ SolverRichardson<> rich(control);
+ SolverRelaxation<> relax(control);
+
+ for (unsigned int size=33; size <= 33; size *= 3)
+ {
+ unsigned int dim = (size-1)*(size-1);
+
+ deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+ // Make matrix
+ FDMatrix testproblem(size, size);
+ SparsityPattern structure(dim, dim, 5);
+ testproblem.five_point_structure(structure);
+ structure.compress();
+ SparseMatrix<double> A(structure);
+ testproblem.five_point(A, true);
+
+ // The permutation vectors;
+ std::vector<unsigned int> perm(dim);
+ std::vector<unsigned int> iperm(dim);
+
+ for (unsigned int blocksize = 4; blocksize < 32; blocksize <<= 1)
+ {
+ deallog << "Block size " << blocksize << std::endl;
+
+ const unsigned int n_blocks = dim/blocksize;
+ BlockList blist;
+ blist.initialize(n_blocks);
+ std::vector<unsigned int> indices(blocksize);
+
+ // The permutation vectors;
+ std::vector<unsigned int> bperm(n_blocks);
+ std::vector<unsigned int> ibperm(n_blocks);
+
+ for (unsigned int block=0;block<n_blocks;++block)
+ {
+ bperm[block] = n_blocks-block-1;
+ for (unsigned int i=0;i<blocksize;++i)
+ indices[i] = i+block*blocksize;
+ blist.add(block, indices);
+ }
+ // Currently, bperm is just
+ // the inversion.
+ // Now swap the last two
+ // and cyclically exchange
+ // the first three
+ unsigned int swap=bperm[n_blocks-1];
+ bperm[n_blocks-1] = bperm[n_blocks-2];
+ bperm[n_blocks-2] = swap;
+
+ swap = bperm[0];
+ bperm[0] = bperm[1];
+ bperm[1] = bperm[2];
+ bperm[2] = swap;
+
+ for (unsigned int block=0;block<n_blocks;++block)
+ {
+ ibperm[bperm[block]] = block;
+ for (unsigned int i=0;i<blocksize;++i)
+ perm[i+block*blocksize] = i+bperm[block]*blocksize;
+ }
+
+ for (unsigned int i=0;i<dim;++i)
+ iperm[perm[i]] = i;
+
+ deallog << "Size " << bperm.size() << std::endl;
+
+ RelaxationBlock<SparseMatrix<double>,double>::AdditionalData relax_data(blist, 0.8);
+ PreconditionBlock<SparseMatrix<double>,double>::AdditionalData prec_data(blocksize, 0.8);
+
+ PreconditionBlockJacobi<SparseMatrix<double>,double> prec_jacobi;
+ prec_jacobi.initialize(A, perm, iperm, prec_data);
+ PreconditionBlockSOR<SparseMatrix<double>,double> prec_sor;
+ prec_sor.initialize(A, perm, iperm, prec_data);
+ PreconditionBlockSSOR<SparseMatrix<double>,double> prec_ssor;
+ prec_ssor.initialize(A, perm, iperm, prec_data);
+
+ PreconditionBlockJacobi<SparseMatrix<double>,double> prec_bjacobi;
+ prec_bjacobi.initialize(A, bperm, ibperm, prec_data);
+ PreconditionBlockSOR<SparseMatrix<double>,double> prec_bsor;
+ prec_bsor.initialize(A, bperm, ibperm, prec_data);
+ PreconditionBlockSSOR<SparseMatrix<double>,double> prec_bssor;
+ prec_bssor.initialize(A, bperm, ibperm, prec_data);
+
+// RelaxationBlockSOR<SparseMatrix<double>,double> relax_sor;
+// relax_sor.set_permutation(bperm,ibperm);
+// relax_sor.initialize(A, relax_data);
+
+// RelaxationBlockSSOR<SparseMatrix<double>,double> relax_ssor;
+// relax_sor.set_permutation(bperm,ibperm);
+// relax_ssor.initialize(A, relax_data);
+
+ Vector<double> f(dim);
+ Vector<double> u(dim);
+ Vector<double> res(dim);
+
+ f = 1.;
+ u = 1.;
+
+ try
+ {
+ double r1, r2;
+
+ deallog.push("Jacobi");
+ r1 = check_solve(rich,A,u,f,prec_jacobi);
+ r2 = check_solve(relax,A,u,f,prec_jacobi);
+ deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+ r2 = check_solve(relax,A,u,f,prec_bjacobi);
+ deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+ deallog.pop();
+
+ deallog.push("SOR ");
+ r1 = check_solve(rich,A,u,f,prec_sor);
+ r2 = check_solve(relax,A,u,f,prec_sor);
+ deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+ r2 = check_solve(relax,A,u,f,prec_bsor);
+ deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+ deallog.pop();
+
+ deallog.push("SSOR ");
+ r1 = check_solve(rich,A,u,f,prec_ssor);
+ r2 = check_solve(relax,A,u,f,prec_ssor);
+ deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+ r2 = check_solve(relax,A,u,f,prec_bssor);
+ deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+ deallog.pop();
+
+ }
+ catch (std::exception& e)
+ {
+ std::cerr << "Exception: " << e.what() << std::endl;
+ }
+ }
+ }
+}
+
--- /dev/null
+
+DEAL::Size 33 Unknowns 1024
+DEAL::Block size 4
+DEAL::Size 256
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 10 value 27.53
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.53
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 10 value 27.53
+DEAL:Jacobi::diff 0
+DEAL:SOR :Richardson::Starting value 32.00
+DEAL:SOR :Richardson::Failure step 10 value 24.60
+DEAL:SOR :Relaxation::Starting value 32.00
+DEAL:SOR :Relaxation::Failure step 10 value 24.60
+DEAL:SOR ::diff 0
+DEAL:SOR :Relaxation::Starting value 32.00
+DEAL:SOR :Relaxation::Failure step 10 value 24.60
+DEAL:SOR ::diff 0
+DEAL:SSOR :Richardson::Starting value 32.00
+DEAL:SSOR :Richardson::Failure step 10 value 21.41
+DEAL:SSOR :Relaxation::Starting value 32.00
+DEAL:SSOR :Relaxation::Failure step 10 value 21.41
+DEAL:SSOR ::diff 0
+DEAL:SSOR :Relaxation::Starting value 32.00
+DEAL:SSOR :Relaxation::Failure step 10 value 21.41
+DEAL:SSOR ::diff 0
+DEAL::Block size 8
+DEAL::Size 128
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 10 value 26.69
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 10 value 26.69
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 10 value 26.69
+DEAL:Jacobi::diff 0
+DEAL:SOR :Richardson::Starting value 32.00
+DEAL:SOR :Richardson::Failure step 10 value 23.66
+DEAL:SOR :Relaxation::Starting value 32.00
+DEAL:SOR :Relaxation::Failure step 10 value 23.66
+DEAL:SOR ::diff 0
+DEAL:SOR :Relaxation::Starting value 32.00
+DEAL:SOR :Relaxation::Failure step 10 value 23.66
+DEAL:SOR ::diff 0
+DEAL:SSOR :Richardson::Starting value 32.00
+DEAL:SSOR :Richardson::Failure step 10 value 19.22
+DEAL:SSOR :Relaxation::Starting value 32.00
+DEAL:SSOR :Relaxation::Failure step 10 value 19.22
+DEAL:SSOR ::diff 0
+DEAL:SSOR :Relaxation::Starting value 32.00
+DEAL:SSOR :Relaxation::Failure step 10 value 19.22
+DEAL:SSOR ::diff 0
+DEAL::Block size 16
+DEAL::Size 64
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 10 value 25.76
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 10 value 25.76
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 10 value 25.76
+DEAL:Jacobi::diff 0
+DEAL:SOR :Richardson::Starting value 32.00
+DEAL:SOR :Richardson::Failure step 10 value 22.95
+DEAL:SOR :Relaxation::Starting value 32.00
+DEAL:SOR :Relaxation::Failure step 10 value 22.95
+DEAL:SOR ::diff 0
+DEAL:SOR :Relaxation::Starting value 32.00
+DEAL:SOR :Relaxation::Failure step 10 value 22.95
+DEAL:SOR ::diff 0
+DEAL:SSOR :Richardson::Starting value 32.00
+DEAL:SSOR :Richardson::Failure step 10 value 17.41
+DEAL:SSOR :Relaxation::Starting value 32.00
+DEAL:SSOR :Relaxation::Failure step 10 value 17.41
+DEAL:SSOR ::diff 0
+DEAL:SSOR :Relaxation::Starting value 32.00
+DEAL:SSOR :Relaxation::Failure step 10 value 17.41
+DEAL:SSOR ::diff 0