* vectors automatically,
* avoiding restart.
*/
-
+
template <class VECTOR>
class TmpVectors
{
* occurs.
*/
VECTOR& operator[] (const unsigned int i) const;
-
+
/**
* Get vector number
* @p i. Allocate it if
*/
VECTOR& operator() (const unsigned int i,
const VECTOR &temp);
-
+
private:
/**
* Pool were vectors are
* obtained from.
*/
VectorMemory<VECTOR> &mem;
-
+
/**
* Field for storing the
* vectors.
*/
std::vector<VECTOR*> data;
-
+
/**
* Offset of the first
* vector. This is for later
* be implemented.
*/
unsigned int offset;
- };
+ };
}
}
* preconditioning. The default is left preconditioning. Finally it
* includes a flag indicating whether or not the default residual is
* used as stopping criterion.
-
+
* <h3>Left versus right preconditioning</h3>
*
* @p AdditionalData allows you to choose between left and right
* pipe additional data to the
* solver.
*/
- struct AdditionalData
+ struct AdditionalData
{
/**
* Constructor. By default, set the
AdditionalData (const unsigned int max_n_tmp_vectors = 30,
const bool right_preconditioning = false,
const bool use_default_residual = true);
-
+
/**
* Maximum number of
* temporary vectors. This
*/
bool use_default_residual;
};
-
+
/**
* Constructor.
*/
<< "The number of temporary vectors you gave ("
<< arg1 << ") is too small. It should be at least 10 for "
<< "any results, and much more for reasonable ones.");
-
+
protected:
/**
* Includes the maximum number of
* tmp vectors.
*/
AdditionalData additional_data;
-
+
/**
* Implementation of the computation of
* the norm of the residual.
*/
virtual double criterion();
-
+
/**
* Transformation of an upper
* Hessenberg matrix into
* Auxiliary matrix for inverting @p H
*/
FullMatrix<double> H1;
-
+
private:
/**
* No copy constructor.
* pipe additional data to the
* solver.
*/
- struct AdditionalData
+ struct AdditionalData
{
/**
* Constructor. By default,
:
max_basis_size(max_basis_size)
{}
-
+
/**
* Maximum number of
* tmp vectors.
*/
unsigned int max_basis_size;
};
-
+
/**
* Constructor.
*/
{
Assert (i+offset<data.size(),
ExcIndexRange(i, -offset, data.size()-offset));
-
+
Assert (data[i-offset] != 0, ExcNotInitialized());
return *data[i-offset];
}
SolverGMRES<VECTOR>::givens_rotation (Vector<double> &h,
Vector<double> &b,
Vector<double> &ci,
- Vector<double> &si,
+ Vector<double> &si,
int col) const
{
for (int i=0 ; i<col ; i++)
h(i) = c*dummy + s*h(i+1);
h(i+1) = -s*dummy + c*h(i+1);
};
-
+
const double r = 1./std::sqrt(h(col)*h(col) + h(col+1)*h(col+1));
si(col) = h(col+1) *r;
ci(col) = h(col) *r;
// don't make any guarantees to its
// optimality or that it even works
// as expected...
-
+
//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");
const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
// Generate an object where basis
// vectors are stored.
internal::SolverGMRES::TmpVectors<VECTOR> tmp_vectors (n_tmp_vectors, this->memory);
-
+
// number of the present iteration; this
// number is not reset to zero upon a
// restart
unsigned int accumulated_iterations = 0;
-
+
// matrix used for the orthogonalization
// process later
H.reinit(n_tmp_vectors, n_tmp_vectors-1);
ci (n_tmp_vectors-1),
si (n_tmp_vectors-1),
h (n_tmp_vectors-1);
-
+
unsigned int dim = 0;
SolverControl::State iteration_state = SolverControl::iterate;
-
+
// switch to determine whether we want a
// left or a right preconditioner. at
// present, left is default, but both
{
// reset this vector to the
// right size
- h.reinit (n_tmp_vectors-1);
-
+ h.reinit (n_tmp_vectors-1);
+
if (left_precondition)
{
A.vmult(p,x);
{
iteration_state = this->control().check (
accumulated_iterations, rho);
-
+
if (iteration_state != SolverControl::iterate)
break;
}
}
else
precondition.vmult(*r,v);
-
+
double res = r->l2_norm();
iteration_state = this->control().check (
accumulated_iterations, res);
break;
}
}
-
+
gamma(0) = rho;
-
+
v *= 1./rho;
// inner iteration doing at
++accumulated_iterations;
// yet another alias
VECTOR& vv = tmp_vectors(inner_iteration+1, x);
-
+
if (left_precondition)
{
A.vmult(p, tmp_vectors[inner_iteration]);
precondition.vmult(p, tmp_vectors[inner_iteration]);
A.vmult(vv,p);
};
-
+
dim = inner_iteration+1;
-
+
/* Orthogonalization */
for (unsigned int i=0 ; i<dim ; ++i)
{
h(i) = vv * tmp_vectors[i];
vv.add(-h(i), tmp_vectors[i]);
};
-
- double s = vv.l2_norm();
- h(inner_iteration+1) = s;
-
- /* Re-orthogonalization */
- for (unsigned i=0; i<dim; ++i)
+
+ /* Re-orthogonalization */
+ for (unsigned int i=0 ; i<dim ; ++i)
{
double htmp = vv * tmp_vectors[i];
h(i) += htmp;
vv.add(-htmp, tmp_vectors[i]);
- };
-
- s = vv.l2_norm();
+ }
+
+ const double s = vv.l2_norm();
h(inner_iteration+1) = s;
-
+
vv *= 1./s;
-
+
/* Transformation into
triagonal structure */
givens_rotation(h,gamma,ci,si,inner_iteration);
-
+
/* append vector on matrix */
for (unsigned int i=0; i<dim; ++i)
H(i,inner_iteration) = h(i);
-
+
/* default residual */
rho = std::fabs(gamma(dim));
else
{
deallog << "default_res=" << rho << std::endl;
-
+
dealii::Vector<double> h_(dim);
*x_=x;
*gamma_=gamma;
H1.reinit(dim+1,dim);
-
+
for (unsigned int i=0; i<dim+1; ++i)
for (unsigned int j=0; j<dim; ++j)
H1(i,j) = H(i,j);
-
+
H1.backward(h_,*gamma_);
-
+
if (left_precondition)
for (unsigned int i=0 ; i<dim; ++i)
x_->add(h_(i), tmp_vectors[i]);
if (left_precondition)
{
const double res=r->l2_norm();
-
+
iteration_state = this->control().check (
accumulated_iterations, res);
}
{
precondition.vmult(*x_, *r);
const double preconditioned_res=x_->l2_norm();
-
+
iteration_state = this->control().check (
accumulated_iterations, preconditioned_res);
}
delete gamma_;
}
-
+
deallog.pop();
// in case of failure: throw
// exception
template<class VECTOR>
double
-SolverGMRES<VECTOR>::criterion ()
+SolverGMRES<VECTOR>::criterion ()
{
// dummy implementation. this function is
// not needed for the present implementation
const PRECONDITIONER& precondition)
{
deallog.push("FGMRES");
-
- SolverControl::State iteration_state = SolverControl::iterate;
+
+ SolverControl::State iteration_state = SolverControl::iterate;
const unsigned int basis_size = additional_data.max_basis_size;
// vectors are stored.
typename internal::SolverGMRES::TmpVectors<VECTOR> v (basis_size, this->memory);
typename internal::SolverGMRES::TmpVectors<VECTOR> z (basis_size, this->memory);
-
+
// number of the present iteration; this
// number is not reset to zero upon a
// restart
unsigned int accumulated_iterations = 0;
-
+
// matrix used for the orthogonalization
// process later
H.reinit(basis_size+1, basis_size);
aux->reinit(x);
A.vmult(*aux, x);
aux->sadd(-1., 1., b);
-
+
double beta = aux->l2_norm();
if (this->control().check(accumulated_iterations,beta)
== SolverControl::success)
break;
-
+
H.reinit(basis_size+1, basis_size);
double a = beta;
-
+
for (unsigned int j=0;j<basis_size;++j)
{
v(j,x).equ(1./a, *aux);
-
+
precondition.vmult(z(j,x), v[j]);
A.vmult(*aux, z[j]);
-
+
// Gram-Schmidt
for (unsigned int i=0;i<=j;++i)
{
aux->add(-H(i,j), v[i]);
}
H(j+1,j) = a = aux->l2_norm();
-
+
// Compute projected solution
-
+
if (j>0)
{
H1.reinit(j+1,j);
projected_rhs.reinit(j+1);
- y.reinit(j);
+ y.reinit(j);
projected_rhs(0) = beta;
H1.fill(H);
Householder<double> house(H1);
// Update solution vector
for (unsigned int j=0;j<y.size();++j)
x.add(y(j), z[j]);
-
+
this->memory.free(aux);
} while (iteration_state == SolverControl::iterate);
-
+
deallog.pop();
// in case of failure: throw
// exception