* //...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);
}
{
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();
}