const VectorType &b,
const PreconditionerType &preconditioner)
{
- deallog.push("Bicgstab");
+ LogStream::Prefix prefix("Bicgstab");
Vr = typename VectorMemory<VectorType>::Pointer(this->memory);
Vrbar = typename VectorMemory<VectorType>::Pointer(this->memory);
Vp = typename VectorMemory<VectorType>::Pointer(this->memory);
}
while (state.breakdown == true);
- deallog.pop();
-
// in case of failure: throw exception
AssertThrow(state.state == SolverControl::success,
SolverControl::NoConvergence (state.last_step,
{
SolverControl::State conv=SolverControl::iterate;
- deallog.push("cg");
+ LogStream::Prefix prefix("cg");
// Memory allocation
typename VectorMemory<VectorType>::Pointer g_pointer(this->memory);
int it=0;
double res = -std::numeric_limits<double>::max();
- try
+ double eigen_beta_alpha = 0;
+
+ // resize the vectors, but do not set
+ // the values since they'd be overwritten
+ // soon anyway.
+ g.reinit(x, true);
+ d.reinit(x, true);
+ h.reinit(x, true);
+
+ double gh,beta;
+
+ // compute residual. if vector is
+ // zero, then short-circuit the
+ // full computation
+ if (!x.all_zero())
{
- double eigen_beta_alpha = 0;
+ A.vmult(g,x);
+ g.add(-1.,b);
+ }
+ else
+ g.equ(-1.,b);
+ res = g.l2_norm();
- // resize the vectors, but do not set
- // the values since they'd be overwritten
- // soon anyway.
- g.reinit(x, true);
- d.reinit(x, true);
- h.reinit(x, true);
+ conv = this->iteration_status(0, res, x);
+ if (conv != SolverControl::iterate)
+ return;
- double gh,beta;
+ if (std::is_same<PreconditionerType,PreconditionIdentity>::value == false)
+ {
+ preconditioner.vmult(h,g);
- // compute residual. if vector is
- // zero, then short-circuit the
- // full computation
- if (!x.all_zero())
- {
- A.vmult(g,x);
- g.add(-1.,b);
- }
- else
- g.equ(-1.,b);
- res = g.l2_norm();
+ d.equ(-1.,h);
+
+ gh = g*h;
+ }
+ else
+ {
+ d.equ(-1.,g);
+ gh = res*res;
+ }
+
+ while (conv == SolverControl::iterate)
+ {
+ it++;
+ A.vmult(h,d);
+
+ double alpha = d*h;
+ Assert(alpha != 0., ExcDivideByZero());
+ alpha = gh/alpha;
+
+ x.add(alpha,d);
+ res = std::sqrt(g.add_and_dot(alpha, h, g));
- conv = this->iteration_status(0, res, x);
+ print_vectors(it, x, g, d);
+
+ conv = this->iteration_status(it, res, x);
if (conv != SolverControl::iterate)
- {
- deallog.pop();
- return;
- }
+ break;
- if (std::is_same<PreconditionerType,PreconditionIdentity>::value == false)
+ if (std::is_same<PreconditionerType,PreconditionIdentity>::value
+ == false)
{
preconditioner.vmult(h,g);
- d.equ(-1.,h);
-
- gh = g*h;
+ beta = gh;
+ Assert(beta != 0., ExcDivideByZero());
+ gh = g*h;
+ beta = gh/beta;
+ d.sadd(beta,-1.,h);
}
else
{
- d.equ(-1.,g);
+ beta = gh;
gh = res*res;
+ beta = gh/beta;
+ d.sadd(beta,-1.,g);
}
- while (conv == SolverControl::iterate)
+ this->coefficients_signal(alpha,beta);
+ // set up the vectors
+ // containing the diagonal
+ // and the off diagonal of
+ // the projected matrix.
+ if (do_eigenvalues)
{
- it++;
- A.vmult(h,d);
-
- double alpha = d*h;
- Assert(alpha != 0., ExcDivideByZero());
- alpha = gh/alpha;
-
- x.add(alpha,d);
- res = std::sqrt(g.add_and_dot(alpha, h, g));
-
- print_vectors(it, x, g, d);
-
- conv = this->iteration_status(it, res, x);
- if (conv != SolverControl::iterate)
- break;
-
- if (std::is_same<PreconditionerType,PreconditionIdentity>::value
- == false)
- {
- preconditioner.vmult(h,g);
-
- beta = gh;
- Assert(beta != 0., ExcDivideByZero());
- gh = g*h;
- beta = gh/beta;
- d.sadd(beta,-1.,h);
- }
- else
- {
- beta = gh;
- gh = res*res;
- beta = gh/beta;
- d.sadd(beta,-1.,g);
- }
-
- this->coefficients_signal(alpha,beta);
- // set up the vectors
- // containing the diagonal
- // and the off diagonal of
- // the projected matrix.
- if (do_eigenvalues)
- {
- diagonal.push_back(1./alpha + eigen_beta_alpha);
- eigen_beta_alpha = beta/alpha;
- offdiagonal.push_back(std::sqrt(beta)/alpha);
- }
- compute_eigs_and_cond(diagonal,offdiagonal,all_eigenvalues_signal,
- all_condition_numbers_signal);
+ diagonal.push_back(1./alpha + eigen_beta_alpha);
+ eigen_beta_alpha = beta/alpha;
+ offdiagonal.push_back(std::sqrt(beta)/alpha);
}
+ compute_eigs_and_cond(diagonal,offdiagonal,all_eigenvalues_signal,
+ all_condition_numbers_signal);
}
- catch (...)
- {
- deallog.pop();
- throw;
- }
+
compute_eigs_and_cond(diagonal,offdiagonal,eigenvalues_signal,
condition_number_signal);
- deallog.pop();
-
// in case of failure: throw exception
if (conv != SolverControl::success)
AssertThrow(false, SolverControl::NoConvergence (it, res));
VectorType &x,
const PreconditionerType &inverse_mass_matrix)
{
- deallog.push("FIRE");
+ LogStream::Prefix prefix("FIRE");
// FIRE algorithm constants
const double DELAYSTEP = 5;
} // While we need to iterate.
- deallog.pop();
-
// In the case of failure: throw exception.
if (conv != SolverControl::success)
AssertThrow (false,
//TODO:[?] Check, why there are two different start residuals.
//TODO:[GK] Make sure the parameter in the constructor means maximum basis size
- deallog.push("GMRES");
+ LogStream::Prefix prefix("GMRES");
const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
// Generate an object where basis vectors are stored.
if (!krylov_space_signal.empty())
krylov_space_signal(tmp_vectors);
- deallog.pop();
-
// in case of failure: throw exception
AssertThrow(iteration_state == SolverControl::success,
SolverControl::NoConvergence (accumulated_iterations,
const VectorType &b,
const PreconditionerType &preconditioner)
{
- deallog.push("FGMRES");
+ LogStream::Prefix prefix("FGMRES");
SolverControl::State iteration_state = SolverControl::iterate;
}
while (iteration_state == SolverControl::iterate);
- deallog.pop();
// in case of failure: throw exception
if (iteration_state != SolverControl::success)
AssertThrow(false, SolverControl::NoConvergence (accumulated_iterations,
const VectorType &b,
const PreconditionerType &preconditioner)
{
- deallog.push("minres");
+ LogStream::Prefix prefix("minres");
// Memory allocation
typename VectorMemory<VectorType>::Pointer Vu0 (this->memory);
delta[1] = delta[2];
}
- // Output
- deallog.pop ();
-
// in case of failure: throw exception
AssertThrow(conv == SolverControl::success,
SolverControl::NoConvergence (j, r_l2));
const VectorType &b,
const PreconditionerType &preconditioner)
{
- deallog.push("QMRS");
+ LogStream::Prefix prefix("QMRS");
// temporary vectors, allocated through the @p VectorMemory object at the
// start of the actual solution process and deallocated at the end.
}
while (state.state == SolverControl::iterate);
- // Output
- deallog.pop();
-
// in case of failure: throw exception
AssertThrow(state.state == SolverControl::success,
SolverControl::NoConvergence (step,
VectorType &d = *Vd;
d.reinit(x);
- deallog.push("Relaxation");
+ LogStream::Prefix prefix("Relaxation");
int iter=0;
- try
+ // Main loop
+ for (; conv==SolverControl::iterate; iter++)
{
- // Main loop
- for (; conv==SolverControl::iterate; iter++)
- {
- // Compute residual
- A.vmult(r,x);
- r.sadd(-1.,1.,b);
-
- // The required norm of the
- // (preconditioned)
- // residual is computed in
- // criterion() and stored
- // in res.
- conv = this->iteration_status (iter, r.l2_norm(), x);
- if (conv != SolverControl::iterate)
- break;
- R.step(x,b);
- }
+ // Compute residual
+ A.vmult(r,x);
+ r.sadd(-1.,1.,b);
+
+ // The required norm of the
+ // (preconditioned)
+ // residual is computed in
+ // criterion() and stored
+ // in res.
+ conv = this->iteration_status (iter, r.l2_norm(), x);
+ if (conv != SolverControl::iterate)
+ break;
+ R.step(x,b);
}
- catch (...)
- {
- deallog.pop();
- throw;
- }
- deallog.pop();
// in case of failure: throw exception
AssertThrow(conv == SolverControl::success,
VectorType &d = *Vd;
d.reinit(x);
- deallog.push("Richardson");
+ LogStream::Prefix prefix("Richardson");
- try
+ // Main loop
+ while (conv==SolverControl::iterate)
{
- // Main loop
- while (conv==SolverControl::iterate)
- {
- // Do not use residual,
- // but do it in 2 steps
- A.vmult(r,x);
- r.sadd(-1.,1.,b);
- preconditioner.vmult(d,r);
-
- // get the required norm of the (possibly preconditioned)
- // residual
- last_criterion = criterion(r, d);
- conv = this->iteration_status (iter, last_criterion, x);
- if (conv != SolverControl::iterate)
- break;
-
- x.add(additional_data.omega,d);
- print_vectors(iter,x,r,d);
-
- ++iter;
- }
- }
- catch (...)
- {
- deallog.pop();
- throw;
+ // Do not use residual,
+ // but do it in 2 steps
+ A.vmult(r,x);
+ r.sadd(-1.,1.,b);
+ preconditioner.vmult(d,r);
+
+ // get the required norm of the (possibly preconditioned)
+ // residual
+ last_criterion = criterion(r, d);
+ conv = this->iteration_status (iter, last_criterion, x);
+ if (conv != SolverControl::iterate)
+ break;
+
+ x.add(additional_data.omega,d);
+ print_vectors(iter,x,r,d);
+
+ ++iter;
}
- deallog.pop();
-
// in case of failure: throw exception
if (conv != SolverControl::success)
AssertThrow(false, SolverControl::NoConvergence (iter,
VectorType &d = *Vd;
d.reinit(x);
- deallog.push("RichardsonT");
+ LogStream::Prefix prefix("RichardsonT");
- try
- {
- // Main loop
- while (conv==SolverControl::iterate)
- {
- // Do not use Tresidual,
- // but do it in 2 steps
- A.Tvmult(r,x);
- r.sadd(-1.,1.,b);
- preconditioner.Tvmult(d,r);
-
- last_criterion = criterion(r, d);
- conv = this->iteration_status (iter, last_criterion, x);
- if (conv != SolverControl::iterate)
- break;
-
- x.add(additional_data.omega,d);
- print_vectors(iter,x,r,d);
-
- ++iter;
- }
- }
- catch (...)
+ // Main loop
+ while (conv==SolverControl::iterate)
{
- deallog.pop();
- throw;
+ // Do not use Tresidual,
+ // but do it in 2 steps
+ A.Tvmult(r,x);
+ r.sadd(-1.,1.,b);
+ preconditioner.Tvmult(d,r);
+
+ last_criterion = criterion(r, d);
+ conv = this->iteration_status (iter, last_criterion, x);
+ if (conv != SolverControl::iterate)
+ break;
+
+ x.add(additional_data.omega,d);
+ print_vectors(iter,x,r,d);
+
+ ++iter;
}
- deallog.pop();
// in case of failure: throw exception
if (conv != SolverControl::success)
AssertThrow(false, SolverControl::NoConvergence (iter, last_criterion));