* //...initialize and build A
*
* std::vector<unsigned int> permutation(x.size());
+ * std::vector<unsigned int> inverse_permutation(x.size());
*
- * //...fill permutation with reasonable values
+ * //...fill permutation and its inverse with reasonable values
*
* // Define and initialize preconditioner
*
* PreconditionPSOR<SparseMatrix<double> > precondition;
- * precondition.initialize (A, permutation, .6);
+ * precondition.initialize (A, permutation, inverse_permutation, .6);
*
* solver.solve (A, x, b, precondition);
* @end{itemize}
*/
void initialize (const MATRIX& A,
const std::vector<unsigned int>& permutation,
+ const std::vector<unsigned int>& inverse_permutation,
typename PreconditionRelaxation<MATRIX>::AdditionalData
parameters = typename PreconditionRelaxation<MATRIX>::AdditionalData());
* Storage for the permutation vector.
*/
const std::vector<unsigned int>* permutation;
+ /**
+ * Storage for the inverse
+ * permutation vector.
+ */
+ const std::vector<unsigned int>* inverse_permutation;
};
PreconditionPSOR<MATRIX>::initialize (
const MATRIX &rA,
const std::vector<unsigned int>& p,
+ const std::vector<unsigned int>& ip,
typename PreconditionRelaxation<MATRIX>::AdditionalData parameters)
{
permutation = &p;
+ inverse_permutation = &ip;
PreconditionRelaxation<MATRIX>::initialize(rA, parameters);
}
{
Assert (this->A!=0, ExcNotInitialized());
dst = src;
- this->A->PSOR (dst, *permutation, this->relaxation);
+ this->A->PSOR (dst, *permutation, *inverse_permutation, this->relaxation);
}
{
Assert (this->A!=0, ExcNotInitialized());
dst = src;
- this->A->TPSOR (dst, *permutation, this->relaxation);
+ this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation);
}
const TYPEMAT) const;
template void SparseMatrix<TYPEMAT>::SSOR<TYPE2> (Vector<TYPE2> &,
const TYPEMAT) const;
+template void SparseMatrix<TYPEMAT>::PSOR<TYPE2> (
+ Vector<TYPE2> &,
+ const std::vector<unsigned int>&,
+ const std::vector<unsigned int>&,
+ const TYPEMAT) const;
+template void SparseMatrix<TYPEMAT>::TPSOR<TYPE2> (
+ Vector<TYPE2> &,
+ const std::vector<unsigned int>&,
+ const std::vector<unsigned int>&,
+ const TYPEMAT) const;
template void SparseMatrix<TYPEMAT>::SOR_step<TYPE2> (Vector<TYPE2> &, const Vector<TYPE2> &,
const TYPEMAT) const;
template void SparseMatrix<TYPEMAT>::TSOR_step<TYPE2> (Vector<TYPE2> &,
* 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 <typename somenumber>
void PSOR (Vector<somenumber> &v,
- const std::vector<unsigned int> permutation,
+ const std::vector<unsigned int>& permutation,
+ const std::vector<unsigned int>& inverse_permutation,
const number om = 1.) const;
/**
* 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 <typename somenumber>
void TPSOR (Vector<somenumber> &v,
- const std::vector<unsigned int> permutation,
+ const std::vector<unsigned int>& permutation,
+ const std::vector<unsigned int>& inverse_permutation,
const number om = 1.) const;
/**
{
somenumber s = dst(row);
for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[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]];
}
}
void
SparseMatrix<number>::PSOR (
Vector<somenumber>& dst,
- const std::vector<unsigned int> permutation,
+ const std::vector<unsigned int>& permutation,
+ const std::vector<unsigned int>& inverse_permutation,
const number om) const
{
Assert (cols != 0, ExcMatrixNotInitialized());
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<m(); ++urow)
{
somenumber s = dst(row);
for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[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);
}
void
SparseMatrix<number>::TPSOR (
Vector<somenumber>& dst,
- const std::vector<unsigned int> permutation,
+ const std::vector<unsigned int>& permutation,
+ const std::vector<unsigned int>& inverse_permutation,
const number om) const
{
Assert (cols != 0, ExcMatrixNotInitialized());
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;)
{
somenumber s = dst(row);
for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[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);
}
template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
void
-check_method( SOLVER& solver, const MATRIX& A,
+check_solve( SOLVER& solver, const MATRIX& A,
VECTOR& u, VECTOR& f, const PRECONDITION& P)
{
u = 0.;
}
}
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+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");
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);
prec_ssor.initialize(A, 1.2);
std::vector<unsigned int> permutation(dim);
+ std::vector<unsigned int> inverse_permutation(dim);
for (unsigned int i=0;i<dim;++i)
- permutation(i) = dim-i-1;
+ permutation[i] = dim-i-1;
+ for (unsigned int i=0;i<dim;++i)
+ inverse_permutation[permutation[i]] = i;
+
+ for (unsigned int i=0;i<dim;++i)
+ permutation[i] = dim-i-1;
PreconditionPSOR<> prec_psor;
- prec_psor.initialize(A, permutation, 1.2);
+ prec_psor.initialize(A, permutation, inverse_permutation, 1.2);
Vector<double> f(dim);
Vector<double> u(dim);
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();
}
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
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
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
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
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
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
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
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