From: Guido Kanschat Date: Fri, 25 Apr 2003 12:52:33 +0000 (+0000) Subject: PSOR X-Git-Tag: v8.0.0~16658 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1b74fbc9550dc4c97dcd5a578a3cd21fe145fa2e;p=dealii.git PSOR git-svn-id: https://svn.dealii.org/trunk@7478 0785d39b-7218-0410-832d-ea1e28bc413d --- 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(); }