solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
protected:
/**
template <typename MatrixType, typename PreconditionerType>
IterationResult
iterate(const MatrixType &A,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
};
/*@}*/
template <typename MatrixType, typename PreconditionerType>
typename SolverBicgstab<VectorType>::IterationResult
SolverBicgstab<VectorType>::iterate(const MatrixType &A,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
//TODO:[GK] Implement "use the length of the computed orthogonal residual" in the BiCGStab method.
SolverControl::State state = SolverControl::iterate;
p.add(-beta*omega, v);
}
- precondition.vmult(y,p);
+ preconditioner.vmult(y,p);
A.vmult(v,y);
rhobar = rbar * v;
return IterationResult(false, SolverControl::success, step, res);
}
- precondition.vmult(z,r);
+ preconditioner.vmult(z,r);
A.vmult(t,z);
rhobar = t*r;
omega = rhobar/(t*t);
SolverBicgstab<VectorType>::solve(const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
deallog.push("Bicgstab");
Vr = this->memory.alloc();
state.state = SolverControl::success;
break;
}
- state = iterate(A, precondition);
+ state = iterate(A, preconditioner);
++step;
}
while (state.breakdown == true);
solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* Connect a slot to retrieve the CG coefficients. The slot will be called
SolverCG<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
SolverControl::State conv=SolverControl::iterate;
if (std::is_same<PreconditionerType,PreconditionIdentity>::value == false)
{
- precondition.vmult(h,g);
+ preconditioner.vmult(h,g);
d.equ(-1.,h);
if (std::is_same<PreconditionerType,PreconditionIdentity>::value
== false)
{
- precondition.vmult(h,g);
+ preconditioner.vmult(h,g);
beta = gh;
Assert(beta != 0., ExcDivideByZero());
void solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
protected:
void SolverFIRE<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
std::function<double(VectorType &, const VectorType &)> compute_func =
return 0.5*A.matrix_norm_square(x) - x*b;
};
- this->solve (compute_func, x, precondition);
+ this->solve (compute_func, x, preconditioner);
}
solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* Connect a slot to retrieve the estimated condition number. Called on each
solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
private:
SolverGMRES<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
// this code was written a very long time ago by people not associated with
// deal.II. we don't make any guarantees to its optimality or that it even
{
A.vmult(p,x);
p.sadd(-1.,1.,b);
- precondition.vmult(v,p);
+ preconditioner.vmult(v,p);
}
else
{
r->sadd(-1.,1.,b);
}
else
- precondition.vmult(*r,v);
+ preconditioner.vmult(*r,v);
double res = r->l2_norm();
last_res = res;
if (left_precondition)
{
A.vmult(p, tmp_vectors[inner_iteration]);
- precondition.vmult(vv,p);
+ preconditioner.vmult(vv,p);
}
else
{
- precondition.vmult(p, tmp_vectors[inner_iteration]);
+ preconditioner.vmult(p, tmp_vectors[inner_iteration]);
A.vmult(vv,p);
}
p = 0.;
for (unsigned int i=0; i<dim; ++i)
p.add(h_(i), tmp_vectors[i]);
- precondition.vmult(*r,p);
+ preconditioner.vmult(*r,p);
x_->add(1.,*r);
};
A.vmult(*r,*x_);
}
else
{
- precondition.vmult(*x_, *r);
+ preconditioner.vmult(*x_, *r);
const double preconditioned_res=x_->l2_norm();
last_res = preconditioned_res;
p = 0.;
for (unsigned int i=0; i<dim; ++i)
p.add(h(i), tmp_vectors[i]);
- precondition.vmult(v,p);
+ preconditioner.vmult(v,p);
x.add(1.,v);
};
// end of outer iteration. restart if no convergence and the number of
SolverFGMRES<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
deallog.push("FGMRES");
v(j,x) = 0.;
- precondition.vmult(z(j,x), v[j]);
+ preconditioner.vmult(z(j,x), v[j]);
A.vmult(*aux, z[j]);
// Gram-Schmidt
solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* @addtogroup Exceptions
SolverMinRes<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
deallog.push("minres");
// positive definite and symmetric
// M v = u[1]
- precondition.vmult (v,*u[1]);
+ preconditioner.vmult (v,*u[1]);
delta[1] = v * (*u[1]);
// Preconditioner positive
// precondition: solve M v = u[2]
// Preconditioner has to be positive
// definite and symmetric.
- precondition.vmult(v,*u[2]);
+ preconditioner.vmult(v,*u[2]);
delta[2] = v * (*u[2]);
* 1e-16.
*/
explicit
- AdditionalData(bool exact_residual = false,
- double breakdown=1.e-16) :
+ AdditionalData(const bool exact_residual = false,
+ const double breakdown = 1.e-16) :
exact_residual(exact_residual),
breakdown(breakdown)
{}
solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* Interface for derived class. This function gets the current iteration
template <typename MatrixType, typename PreconditionerType>
IterationResult
iterate (const MatrixType &A,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* Number of the current iteration (accumulated over restarts)
SolverQMRS<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
deallog.push("QMRS");
{
if (step > 0)
deallog << "Restart step " << step << std::endl;
- state = iterate(A, precondition);
+ state = iterate(A, preconditioner);
}
while (state.state == SolverControl::iterate);
template <typename MatrixType, typename PreconditionerType>
typename SolverQMRS<VectorType>::IterationResult
SolverQMRS<VectorType>::iterate(const MatrixType &A,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
/* Remark: the matrix A in the article is the preconditioned matrix.
* Therefore, we have to precondition x before we compute the first residual.
d.reinit(x);
// Apply right preconditioning to x
- precondition.vmult(q,x);
+ preconditioner.vmult(q,x);
// Preconditioned residual
A.vmult(v,q);
v.sadd(-1.,1.,b);
p = v;
- precondition.vmult(q,p);
+ preconditioner.vmult(q,p);
tau = v.norm_sqr();
rho = q*v;
return IterationResult(SolverControl::iterate, std::fabs(rho));
// Step 7
const double rho_old = rho;
- precondition.vmult(q,v);
+ preconditioner.vmult(q,v);
rho = q*v;
const double beta = rho/rho_old;
p.sadd(beta,v);
- precondition.vmult(q,p);
+ preconditioner.vmult(q,p);
}
return IterationResult(SolverControl::success, std::fabs(rho));
}
solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* Solve $A^Tx=b$ for $x$.
Tsolve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition);
+ const PreconditionerType &preconditioner);
/**
* Set the damping-coefficient. Default is 1., i.e. no damping.
SolverRichardson<VectorType>::solve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
SolverControl::State conv=SolverControl::iterate;
// but do it in 2 steps
A.vmult(r,x);
r.sadd(-1.,1.,b);
- precondition.vmult(d,r);
+ preconditioner.vmult(d,r);
// The required norm of the
// (preconditioned)
SolverRichardson<VectorType>::Tsolve (const MatrixType &A,
VectorType &x,
const VectorType &b,
- const PreconditionerType &precondition)
+ const PreconditionerType &preconditioner)
{
SolverControl::State conv=SolverControl::iterate;
double last_criterion = -std::numeric_limits<double>::max();
// but do it in 2 steps
A.Tvmult(r,x);
r.sadd(-1.,1.,b);
- precondition.Tvmult(d,r);
+ preconditioner.Tvmult(d,r);
last_criterion = criterion();
conv = this->iteration_status (iter, last_criterion, x);