From def0748424f34e208eef7501b290a302fe60d523 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 25 Apr 2003 12:52:33 +0000 Subject: [PATCH] PSOR git-svn-id: https://svn.dealii.org/trunk@7478 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/precondition.h | 17 +++- .../lac/include/lac/sparse_matrix.2.templates | 10 +++ deal.II/lac/include/lac/sparse_matrix.h | 14 +++- .../lac/include/lac/sparse_matrix.templates.h | 27 ++++-- tests/lac/solver.cc | 83 +++++++++++++------ .../lac/solver.output | 53 +++++++++++- 6 files changed, 162 insertions(+), 42 deletions(-) diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h index b30f2d119f..7259f27fce 100644 --- a/deal.II/lac/include/lac/precondition.h +++ b/deal.II/lac/include/lac/precondition.h @@ -374,13 +374,14 @@ class PreconditionSSOR : public PreconditionRelaxation * //...initialize and build A * * std::vector permutation(x.size()); + * std::vector inverse_permutation(x.size()); * - * //...fill permutation with reasonable values + * //...fill permutation and its inverse with reasonable values * * // Define and initialize preconditioner * * PreconditionPSOR > precondition; - * precondition.initialize (A, permutation, .6); + * precondition.initialize (A, permutation, inverse_permutation, .6); * * solver.solve (A, x, b, precondition); * @end{itemize} @@ -412,6 +413,7 @@ class PreconditionPSOR : public PreconditionRelaxation */ void initialize (const MATRIX& A, const std::vector& permutation, + const std::vector& inverse_permutation, typename PreconditionRelaxation::AdditionalData parameters = typename PreconditionRelaxation::AdditionalData()); @@ -432,6 +434,11 @@ class PreconditionPSOR : public PreconditionRelaxation * Storage for the permutation vector. */ const std::vector* permutation; + /** + * Storage for the inverse + * permutation vector. + */ + const std::vector* inverse_permutation; }; @@ -712,9 +719,11 @@ inline void PreconditionPSOR::initialize ( const MATRIX &rA, const std::vector& p, + const std::vector& ip, typename PreconditionRelaxation::AdditionalData parameters) { permutation = &p; + inverse_permutation = &ip; PreconditionRelaxation::initialize(rA, parameters); } @@ -726,7 +735,7 @@ PreconditionPSOR::vmult (VECTOR& dst, const VECTOR& src) const { Assert (this->A!=0, ExcNotInitialized()); dst = src; - this->A->PSOR (dst, *permutation, this->relaxation); + this->A->PSOR (dst, *permutation, *inverse_permutation, this->relaxation); } @@ -738,7 +747,7 @@ PreconditionPSOR::Tvmult (VECTOR& dst, const VECTOR& src) const { Assert (this->A!=0, ExcNotInitialized()); dst = src; - this->A->TPSOR (dst, *permutation, this->relaxation); + this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation); } diff --git a/deal.II/lac/include/lac/sparse_matrix.2.templates b/deal.II/lac/include/lac/sparse_matrix.2.templates index caa15d8dcb..0f664b8a34 100644 --- a/deal.II/lac/include/lac/sparse_matrix.2.templates +++ b/deal.II/lac/include/lac/sparse_matrix.2.templates @@ -67,6 +67,16 @@ template void SparseMatrix::TSOR (Vector &, const TYPEMAT) const; template void SparseMatrix::SSOR (Vector &, const TYPEMAT) const; +template void SparseMatrix::PSOR ( + Vector &, + const std::vector&, + const std::vector&, + const TYPEMAT) const; +template void SparseMatrix::TPSOR ( + Vector &, + const std::vector&, + const std::vector&, + const TYPEMAT) const; template void SparseMatrix::SOR_step (Vector &, const Vector &, const TYPEMAT) const; template void SparseMatrix::TSOR_step (Vector &, diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index e69fc38545..b35345df64 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -814,14 +814,17 @@ class SparseMatrix : public Subscriptor * that is, first the row * @p{permutation[0]}, then * @p{permutation[1]} and so - * on. + * on. For efficiency reasons, + * the permutation as well as its + * inverse are required. * * @p{omega} is the relaxation * parameter. */ template void PSOR (Vector &v, - const std::vector permutation, + const std::vector& permutation, + const std::vector& inverse_permutation, const number om = 1.) const; /** @@ -834,14 +837,17 @@ class SparseMatrix : public Subscriptor * that is, first the row * @p{permutation[m()-1]}, then * @p{permutation[m()-2]} and so - * on. + * on. For efficiency reasons, + * the permutation as well as its + * inverse are required. * * @p{omega} is the relaxation * parameter. */ template void TPSOR (Vector &v, - const std::vector permutation, + const std::vector& permutation, + const std::vector& inverse_permutation, const number om = 1.) const; /** diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index f5cf298da6..3a9659c2db 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -962,9 +962,12 @@ SparseMatrix::SOR (Vector& dst, const number om) const { somenumber s = dst(row); for (unsigned int j=cols->rowstart[row]; jrowstart[row+1]; ++j) - if (cols->colnums[j] < row) - s -= val[j] * dst(cols->colnums[j]); - + { + const unsigned int col = cols->colnums[j]; + if (col < row) + s -= val[j] * dst(col); + } + dst(row) = s * om / val[cols->rowstart[row]]; } } @@ -998,7 +1001,8 @@ template void SparseMatrix::PSOR ( Vector& dst, - const std::vector permutation, + const std::vector& permutation, + const std::vector& inverse_permutation, const number om) const { Assert (cols != 0, ExcMatrixNotInitialized()); @@ -1007,6 +1011,8 @@ SparseMatrix::PSOR ( Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert (m() == permutation.size(), ExcDimensionMismatch(m(), permutation.size())); + Assert (m() == inverse_permutation.size(), + ExcDimensionMismatch(m(), inverse_permutation.size())); for (unsigned int urow=0; urow::PSOR ( somenumber s = dst(row); for (unsigned int j=cols->rowstart[row]; jrowstart[row+1]; ++j) { - const unsigned int col = permutation[cols->colnums[j]]; - if (col < row) + const unsigned int col = cols->colnums[j]; + if (inverse_permutation[col] < urow) s -= val[j] * dst(col); } @@ -1029,7 +1035,8 @@ template void SparseMatrix::TPSOR ( Vector& dst, - const std::vector permutation, + const std::vector& permutation, + const std::vector& inverse_permutation, const number om) const { Assert (cols != 0, ExcMatrixNotInitialized()); @@ -1038,6 +1045,8 @@ SparseMatrix::TPSOR ( Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert (m() == permutation.size(), ExcDimensionMismatch(m(), permutation.size())); + Assert (m() == inverse_permutation.size(), + ExcDimensionMismatch(m(), inverse_permutation.size())); for (unsigned int urow=m(); urow != 0;) { @@ -1046,8 +1055,8 @@ SparseMatrix::TPSOR ( somenumber s = dst(row); for (unsigned int j=cols->rowstart[row]; jrowstart[row+1]; ++j) { - const unsigned int col = permutation[cols->colnums[j]]; - if (col > row) + const unsigned int col = cols->colnums[j]; + if (inverse_permutation[col] > urow) s -= val[j] * dst(col); } diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index dc54ef4c80..538d9c71b6 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -31,7 +31,7 @@ template void -check_method( SOLVER& solver, const MATRIX& A, +check_solve( SOLVER& solver, const MATRIX& A, VECTOR& u, VECTOR& f, const PRECONDITION& P) { u = 0.; @@ -46,6 +46,23 @@ check_method( SOLVER& solver, const MATRIX& A, } } +template +void +check_Tsolve(SOLVER& solver, const MATRIX& A, + VECTOR& u, VECTOR& f, const PRECONDITION& P) +{ + u = 0.; + f = 1.; + try + { + solver.Tsolve(A,u,f,P); + } + catch (std::exception& e) + { + deallog << e.what() << std::endl; + } +} + int main() { std::ofstream logfile("solver.output"); @@ -59,6 +76,8 @@ int main() SolverControl verbose_control(100, 1.e-3, true); SolverCG<> cg(control, mem); SolverGMRES<> gmres(control, mem, 8); + SolverGMRES<>::AdditionalData(8, true); + SolverGMRES<> gmresright(control, mem, 8); SolverBicgstab<> bicgstab(control, mem); SolverRichardson<> rich(control, mem); SolverQMRS<> qmrs(control, mem); @@ -84,11 +103,17 @@ int main() prec_ssor.initialize(A, 1.2); std::vector permutation(dim); + std::vector inverse_permutation(dim); for (unsigned int i=0;i prec_psor; - prec_psor.initialize(A, permutation, 1.2); + prec_psor.initialize(A, permutation, inverse_permutation, 1.2); Vector f(dim); Vector u(dim); @@ -110,46 +135,56 @@ int main() deallog.push("no-fail"); control.set_max_steps(10); - check_method(cg,A,u,f,prec_no); - check_method(bicgstab,A,u,f,prec_no); - check_method(gmres,A,u,f,prec_no); - check_method(qmrs,A,u,f,prec_no); + check_solve(cg,A,u,f,prec_no); + check_solve(bicgstab,A,u,f,prec_no); + check_solve(gmres,A,u,f,prec_no); + check_solve(gmresright,A,u,f,prec_no); + check_solve(qmrs,A,u,f,prec_no); control.set_max_steps(100); deallog.pop(); deallog.push("no"); - check_method(cg,A,u,f,prec_no); - check_method(bicgstab,A,u,f,prec_no); - check_method(gmres,A,u,f,prec_no); - check_method(qmrs,A,u,f,prec_no); + check_solve(cg,A,u,f,prec_no); + check_solve(bicgstab,A,u,f,prec_no); + check_solve(gmres,A,u,f,prec_no); + check_solve(gmresright,A,u,f,prec_no); + check_solve(qmrs,A,u,f,prec_no); deallog.pop(); deallog.push("ssor"); - check_method(rich,A,u,f,prec_ssor); - check_method(cg,A,u,f,prec_ssor); - check_method(bicgstab,A,u,f,prec_ssor); - check_method(gmres,A,u,f,prec_ssor); - check_method(qmrs,A,u,f,prec_ssor); + check_Tsolve(rich,A,u,f,prec_ssor); + check_solve(rich,A,u,f,prec_ssor); + check_solve(cg,A,u,f,prec_ssor); + check_solve(bicgstab,A,u,f,prec_ssor); + check_solve(gmres,A,u,f,prec_ssor); + check_solve(gmresright,A,u,f,prec_ssor); + check_solve(qmrs,A,u,f,prec_ssor); deallog.pop(); deallog.push("sor"); - check_method(rich,A,u,f,prec_sor); - check_method(cg,A,u,f,prec_sor); - check_method(bicgstab,A,u,f,prec_sor); - check_method(gmres,A,u,f,prec_sor); + check_Tsolve(rich,A,u,f,prec_sor); + check_solve(rich,A,u,f,prec_sor); + check_solve(cg,A,u,f,prec_sor); + check_solve(bicgstab,A,u,f,prec_sor); + check_solve(gmres,A,u,f,prec_sor); + check_solve(gmresright,A,u,f,prec_sor); + + deallog.pop(); deallog.push("psor"); - check_method(rich,A,u,f,prec_psor); - check_method(cg,A,u,f,prec_psor); - check_method(bicgstab,A,u,f,prec_psor); - check_method(gmres,A,u,f,prec_psor); + check_Tsolve(rich,A,u,f,prec_psor); + check_solve(rich,A,u,f,prec_psor); + check_solve(cg,A,u,f,prec_psor); + check_solve(bicgstab,A,u,f,prec_psor); + check_solve(gmres,A,u,f,prec_psor); + check_solve(gmresright,A,u,f,prec_psor); deallog.pop(); } diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/lac/solver.output b/tests/results/i686-pc-linux-gnu+gcc3.2/lac/solver.output index 4513f30b2e..a964277c44 100644 --- a/tests/results/i686-pc-linux-gnu+gcc3.2/lac/solver.output +++ b/tests/results/i686-pc-linux-gnu+gcc3.2/lac/solver.output @@ -7,6 +7,8 @@ DEAL:no-fail:Bicgstab::Starting value 3.000 DEAL:no-fail:Bicgstab::Convergence step 3 value 7.850e-16 DEAL:no-fail:GMRES::Starting value 3.000 DEAL:no-fail:GMRES::Convergence step 3 value 2.266e-17 +DEAL:no-fail:GMRES::Starting value 3.000 +DEAL:no-fail:GMRES::Convergence step 3 value 2.266e-17 DEAL:no-fail:QMRS::Starting value 3.000 DEAL:no-fail:QMRS::Convergence step 3 value 1.183e-16 DEAL:no:cg::Starting value 3.000 @@ -15,8 +17,12 @@ DEAL:no:Bicgstab::Starting value 3.000 DEAL:no:Bicgstab::Convergence step 3 value 7.850e-16 DEAL:no:GMRES::Starting value 3.000 DEAL:no:GMRES::Convergence step 3 value 2.266e-17 +DEAL:no:GMRES::Starting value 3.000 +DEAL:no:GMRES::Convergence step 3 value 2.266e-17 DEAL:no:QMRS::Starting value 3.000 DEAL:no:QMRS::Convergence step 3 value 1.183e-16 +DEAL:ssor:RichardsonT::Starting value 3.000 +DEAL:ssor:RichardsonT::Convergence step 10 value 0.0004418 DEAL:ssor:Richardson::Starting value 3.000 DEAL:ssor:Richardson::Convergence step 10 value 0.0004418 DEAL:ssor:cg::Starting value 3.000 @@ -25,8 +31,12 @@ DEAL:ssor:Bicgstab::Starting value 3.000 DEAL:ssor:Bicgstab::Convergence step 2 value 0.0003670 DEAL:ssor:GMRES::Starting value 1.600 DEAL:ssor:GMRES::Convergence step 3 value 0.0002978 +DEAL:ssor:GMRES::Starting value 1.600 +DEAL:ssor:GMRES::Convergence step 3 value 0.0002978 DEAL:ssor:QMRS::Starting value 3.000 DEAL:ssor:QMRS::Convergence step 4 value 0.0001345 +DEAL:sor:RichardsonT::Starting value 3.000 +DEAL:sor:RichardsonT::Convergence step 7 value 0.0004339 DEAL:sor:Richardson::Starting value 3.000 DEAL:sor:Richardson::Convergence step 7 value 0.0004339 DEAL:sor:cg::Starting value 3.000 @@ -36,6 +46,21 @@ DEAL:sor:Bicgstab::Starting value 3.000 DEAL:sor:Bicgstab::Convergence step 5 value 1.632e-15 DEAL:sor:GMRES::Starting value 1.462 DEAL:sor:GMRES::Convergence step 5 value 1.278e-18 +DEAL:sor:GMRES::Starting value 1.462 +DEAL:sor:GMRES::Convergence step 5 value 1.278e-18 +DEAL:psor:RichardsonT::Starting value 3.000 +DEAL:psor:RichardsonT::Convergence step 7 value 0.0004339 +DEAL:psor:Richardson::Starting value 3.000 +DEAL:psor:Richardson::Convergence step 7 value 0.0004339 +DEAL:psor:cg::Starting value 3.000 +DEAL:psor:cg::Failure step 100 value 0.2585 +DEAL:psor::Iterative method reported convergence failure in step 100 with residual 0.258509 +DEAL:psor:Bicgstab::Starting value 3.000 +DEAL:psor:Bicgstab::Convergence step 5 value 1.797e-15 +DEAL:psor:GMRES::Starting value 1.462 +DEAL:psor:GMRES::Convergence step 5 value 1.208e-18 +DEAL:psor:GMRES::Starting value 1.462 +DEAL:psor:GMRES::Convergence step 5 value 1.208e-18 DEAL::Size 12 Unknowns 121 DEAL::SOR-diff:0.000 DEAL:no-fail:cg::Starting value 11.00 @@ -47,6 +72,9 @@ DEAL:no-fail::Iterative method reported convergence failure in step 10 with resi DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.841354 +DEAL:no-fail:GMRES::Starting value 11.00 +DEAL:no-fail:GMRES::Failure step 10 value 0.8414 +DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.841354 DEAL:no-fail:QMRS::Starting value 11.00 DEAL:no-fail:QMRS::Failure step 10 value 0.4215 DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.421504 @@ -56,8 +84,12 @@ DEAL:no:Bicgstab::Starting value 11.00 DEAL:no:Bicgstab::Convergence step 11 value 0.0003990 DEAL:no:GMRES::Starting value 11.00 DEAL:no:GMRES::Convergence step 43 value 0.0009788 +DEAL:no:GMRES::Starting value 11.00 +DEAL:no:GMRES::Convergence step 43 value 0.0009788 DEAL:no:QMRS::Starting value 11.00 DEAL:no:QMRS::Convergence step 16 value 0.0002583 +DEAL:ssor:RichardsonT::Starting value 11.00 +DEAL:ssor:RichardsonT::Convergence step 59 value 0.0008892 DEAL:ssor:Richardson::Starting value 11.00 DEAL:ssor:Richardson::Convergence step 59 value 0.0008892 DEAL:ssor:cg::Starting value 11.00 @@ -66,8 +98,12 @@ DEAL:ssor:Bicgstab::Starting value 11.00 DEAL:ssor:Bicgstab::Convergence step 5 value 9.315e-05 DEAL:ssor:GMRES::Starting value 11.05 DEAL:ssor:GMRES::Convergence step 7 value 0.0006838 +DEAL:ssor:GMRES::Starting value 11.05 +DEAL:ssor:GMRES::Convergence step 7 value 0.0006838 DEAL:ssor:QMRS::Starting value 11.00 DEAL:ssor:QMRS::Convergence step 9 value 0.0004140 +DEAL:sor:RichardsonT::Starting value 11.00 +DEAL:sor:RichardsonT::Convergence step 88 value 0.0009636 DEAL:sor:Richardson::Starting value 11.00 DEAL:sor:Richardson::Convergence step 88 value 0.0009636 DEAL:sor:cg::Starting value 11.00 @@ -77,5 +113,20 @@ DEAL:sor:Bicgstab::Starting value 11.00 DEAL:sor:Bicgstab::Convergence step 14 value 0.0009987 DEAL:sor:GMRES::Starting value 7.322 DEAL:sor:GMRES::Convergence step 23 value 0.0006981 -DEAL::GrowingVectorMemory:Overall allocated vectors: 188 +DEAL:sor:GMRES::Starting value 7.322 +DEAL:sor:GMRES::Convergence step 23 value 0.0006981 +DEAL:psor:RichardsonT::Starting value 11.00 +DEAL:psor:RichardsonT::Convergence step 88 value 0.0009636 +DEAL:psor:Richardson::Starting value 11.00 +DEAL:psor:Richardson::Convergence step 88 value 0.0009636 +DEAL:psor:cg::Starting value 11.00 +DEAL:psor:cg::Failure step 100 value 5.643 +DEAL:psor::Iterative method reported convergence failure in step 100 with residual 5.64329 +DEAL:psor:Bicgstab::Starting value 11.00 +DEAL:psor:Bicgstab::Convergence step 14 value 0.0009987 +DEAL:psor:GMRES::Starting value 7.322 +DEAL:psor:GMRES::Convergence step 23 value 0.0006981 +DEAL:psor:GMRES::Starting value 7.322 +DEAL:psor:GMRES::Convergence step 23 value 0.0006981 +DEAL::GrowingVectorMemory:Overall allocated vectors: 312 DEAL::GrowingVectorMemory:Maximum allocated vectors: 8 -- 2.39.5