Commit a6b85d35 authored by Praetorius, Simon's avatar Praetorius, Simon

gmres solvers improved: householder orthogonalization added; corrections in...

gmres solvers improved: householder orthogonalization added; corrections in sequential PETSc solvers/matrixstreams
parent 806b9b01
project(AMDIS)
cmake_minimum_required(VERSION 2.6)
cmake_policy(SET CMP0017 OLD)
SET(LIB_DIR ${AMDIS_SOURCE_DIR}/lib)
SET(SOURCE_DIR ${AMDIS_SOURCE_DIR}/src)
......@@ -416,7 +417,9 @@ if(ENABLE_EXTENSIONS)
${EXTENSIONS_DIR}/POperators.cc
${EXTENSIONS_DIR}/SingularDirichletBC2.cc
${EXTENSIONS_DIR}/time/ExtendedRosenbrockStationary.cc
${EXTENSIONS_DIR}/pugixml/src/pugixml.cpp)
${EXTENSIONS_DIR}/pugixml/src/pugixml.cpp
${EXTENSIONS_DIR}/preconditioner/PhaseFieldCrystal_.cc
${EXTENSIONS_DIR}/preconditioner/CahnHilliard_.cc)
if(ENABLE_PETSC)
list(APPEND EXTENSIONS_SRC
......@@ -427,9 +430,7 @@ if(ENABLE_EXTENSIONS)
if(ENABLE_PARALLEL_DOMAIN)
list(APPEND EXTENSIONS_SRC
# ${EXTENSIONS_DIR}/preconditioner/PetscSolverNavierStokes2.cc
${EXTENSIONS_DIR}/preconditioner/PetscSolverPfc.cc
${EXTENSIONS_DIR}/preconditioner/PhaseFieldCrystal_.cc
${EXTENSIONS_DIR}/preconditioner/CahnHilliard_.cc)
${EXTENSIONS_DIR}/preconditioner/PetscSolverPfc.cc)
endif(ENABLE_PARALLEL_DOMAIN)
list(APPEND COMPILEFLAGS "-DHAVE_EXTENSIONS=1")
......@@ -708,8 +709,6 @@ list(APPEND deb_add_dirs "share/amdis/doc")
list(REMOVE_DUPLICATES deb_add_dirs)
# TESTING
ENABLE_TESTING()
INCLUDE(CTest)
INCLUDE(Dart)
......
......@@ -56,6 +56,7 @@ namespace AMDiS {
const ElementVector& uh_loc,
WorldMatrix<double>* D2_uh) const
{
// TODO: REMOVE STATIC
static WorldMatrix<double> D2(DEFAULT_VALUE, 0.0);
DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0);
DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0);
......
......@@ -296,7 +296,7 @@ namespace AMDiS {
{
FUNCNAME_DBG("CoarseningManager3d::getCoarsenPatch()");
static unsigned char next_el[6][2] = {{3,2},
static const unsigned char next_el[6][2] = {{3,2},
{1,3},
{1,2},
{0,3},
......
......@@ -531,7 +531,7 @@ namespace AMDiS {
DOFAdmin *admin = feSpace->getAdmin();
// define result vector
static WorldVector<DOFVector<double>*> *result = nullptr;
static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
if (grad) {
result = grad;
......@@ -630,7 +630,7 @@ namespace AMDiS {
TEST_EXIT_DBG(vec)("no vector\n");
int dow = Global::getGeo(WORLD);
static WorldVector<DOFVector<double>*> *result = nullptr;
static WorldVector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
if (!res && !result) {
result = new WorldVector<DOFVector<double>*>;
......
......@@ -1247,15 +1247,15 @@ namespace AMDiS {
template<typename T>
const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
{
static DOFVector<T> result;
return mult(d, v, result);
static DOFVector<T> result; // TODO: REMOVE STATIC
return mult(d, v, result);
}
template<typename T>
const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
{
static DOFVector<T> result;
static DOFVector<T> result; // TODO: REMOVE STATIC
return mult(d, v, result);
}
......@@ -1263,7 +1263,7 @@ namespace AMDiS {
template<typename T>
const DOFVector<T>& operator+(const DOFVector<T> &v1 , const DOFVector<T> &v2)
{
static DOFVector<T> result;
static DOFVector<T> result; // TODO: REMOVE STATIC
return add(v1, v2, result);
}
......@@ -1765,7 +1765,7 @@ namespace AMDiS {
const FiniteElemSpace *feSpace = DOFVector<T>::feSpace;
// define result vector
static DOFVector<typename GradientType<T>::type> *result = nullptr;
static DOFVector<typename GradientType<T>::type> *result = nullptr; // TODO: REMOVE STATIC
if (grad) {
result = grad;
......@@ -1850,7 +1850,7 @@ namespace AMDiS {
int dim = DOFVector<T>::dim;
// define result vector
static DOFVector<typename GradientType<T>::type> *vec = nullptr;
static DOFVector<typename GradientType<T>::type> *vec = nullptr; // TODO: REMOVE STATIC
DOFVector<typename GradientType<T>::type> *result = grad;
......@@ -1922,7 +1922,7 @@ namespace AMDiS {
TEST_EXIT_DBG(vec)("no vector\n");
static std::vector<DOFVector<double>*> *result = nullptr;
static std::vector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
int len = num_rows(GradientType<T>::getValues((*vec)[0]));
......
......@@ -795,7 +795,7 @@ Recovery::recoveryUh(DOFVector<double> *uh, const FiniteElemSpace *fe_space)
}
// define result vector
static DOFVector<double> *vec = nullptr;
static DOFVector<double> *vec = nullptr;// TODO: REMOVE STATIC
DOFVector<double> *result = nullptr;
// Allocate memory for result vector
......@@ -868,7 +868,7 @@ Recovery::recovery(DOFVector<double> *uh, const FiniteElemSpace *fe_space,
}
// define result vector
static DOFVector<WorldVector<double> > *vec = nullptr;
static DOFVector<WorldVector<double> > *vec = nullptr;// TODO: REMOVE STATIC
DOFVector<WorldVector<double> > *result = nullptr;
// Allocate memory for result vector
......@@ -920,7 +920,7 @@ Recovery::recovery(DOFVector<double> *uh,
const FiniteElemSpace *fe_space = uh->getFeSpace();
// define result vector
static DOFVector<WorldVector<double> > *vec = nullptr;
static DOFVector<WorldVector<double> > *vec = nullptr;// TODO: REMOVE STATIC
DOFVector<WorldVector<double> > *result = nullptr;
// Allocate memory for result vector
......
......@@ -131,6 +131,7 @@ namespace AMDiS {
void Tetrahedron::sortFaceIndices(int face, FixVec<int,WORLD> &vec) const
{
// TODO: REMOVE STATIC
static MatrixOfFixVecs<FixVec<int,WORLD> > *sorted_3d = nullptr;
if (sorted_3d == nullptr) {
......
......@@ -66,7 +66,7 @@ namespace AMDiS {
virtual int getPositionOfVertex(int side, int vertex) const
{
static int positionOfVertex[4][4] = {{-1, 0, 1, 2},
static const int positionOfVertex[4][4] = {{-1, 0, 1, 2},
{0, -1, 1, 2},
{0, 1, -1, 2},
{0, 1, 2, -1}};
......
......@@ -536,6 +536,7 @@ namespace AMDiS {
int sav_neighbour = neighbour;
// father.neigh[coarse_nb[i][j]] == child[i - 1].neigh[j]
// TODO: REMOVE STATIC
static int coarse_nb[3][3][4] = {{{-2, -2, -2, -2}, {-1, 2, 3, 1}, {-1, 3, 2, 0}},
{{-2, -2, -2, -2}, {-1, 2, 3, 1}, {-1, 2, 3, 0}},
{{-2, -2, -2, -2}, {-1, 2, 3, 1}, {-1, 2, 3, 0}}};
......@@ -799,6 +800,7 @@ namespace AMDiS {
int sav_neighbour = neighbour;
// father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j]
// TODO: REMOVE STATIC
static int coarse_nb[3][3] = {{-2, -2, -2}, {2, -1, 1}, {-1, 2, 0}};
TEST_EXIT_DBG(stack_used > 0)("no current element");
......
......@@ -62,6 +62,7 @@ namespace AMDiS {
void Triangle::sortFaceIndices(int face, FixVec<int, WORLD> &vec) const
{
// TODO: REMOVE STATIC
static MatrixOfFixVecs<FixVec<int, WORLD> > *sorted_2d = nullptr;
if (sorted_2d == nullptr) {
......
......@@ -160,7 +160,7 @@ namespace AMDiS {
TEST_EXIT_DBG(side >= 0 && side <= 2)("Wrong side number %d!\n", side);
TEST_EXIT_DBG(vertex >= 0 && vertex <= 2)("Wrong vertex number %d!\n", vertex);
static int positionOfVertex[3][3] = {{-1, 0, 1}, {1, -1, 0}, {0, 1, -1}};
static const int positionOfVertex[3][3] = {{-1, 0, 1}, {1, -1, 0}, {0, 1, -1}};
return positionOfVertex[side][vertex];
}
......
......@@ -52,7 +52,8 @@ namespace AMDiS {
ofstream outfile;
outfile.open(filename.c_str());
outfile.precision(10);
outfile.setf( std::ios::scientific, std:: ios::floatfield );
outfile.precision(16);
DOFIterator<WorldVector<double> > it(&coordDof, USED_DOFS);
for (it.reset(); !it.end(); ++it) {
......
......@@ -144,7 +144,7 @@ namespace AMDiS { namespace Parallel {
int dim = componentSpaces[pressureComponent]->getMesh()->getDim();
vector<int> velocityComponents;
for (size_t i = 0; i < dim; i++)
for (int i = 0; i < dim; i++)
velocityComponents.push_back(i);
PCSetType(pc, PCFIELDSPLIT);
......
......@@ -41,7 +41,9 @@
#include "solver/itl/minres.hpp"
#include "solver/itl/gcr.hpp"
#include "solver/itl/fgmres.hpp"
#include "solver/itl/fgmres_householder.hpp"
#include "solver/itl/gmres2.hpp"
#include "solver/itl/gmres_householder.hpp"
#include "solver/itl/preonly.hpp"
......@@ -264,20 +266,37 @@ namespace AMDiS {
* The method is not preconditioned
*/
enum ORTHOGONALIZATION {
GRAM_SCHMIDT = 1,
HOUSEHOLDER = 2
};
class gmres_type
{
int restart;
int restart;
ORTHOGONALIZATION orthogonalization;
public:
gmres_type(std::string name) : restart(30)
gmres_type(std::string name) : restart(30), orthogonalization(GRAM_SCHMIDT)
{
Parameters::get(name + "->restart", restart);
Parameters::get(name + "->orthogonalization", orthogonalization);
}
template < class LinOp, class X, class B, class L, class R, class I >
int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter)
{ return itl::gmres2(A, x, b, l, r, iter, restart); }
{
switch (orthogonalization) {
default:
case GRAM_SCHMIDT:
return itl::gmres2(A, x, b, l, r, iter, restart); break;
case HOUSEHOLDER:
return itl::gmres_householder(A, x, b, l, iter, restart); break;
}
}
};
typedef ITL_Solver< gmres_type > GMResSolver;
// ===================================================================================
/**
......@@ -343,6 +362,7 @@ namespace AMDiS {
class gcr_type
{
int restart;
public:
gcr_type(std::string name) : restart(30)
{
......@@ -369,14 +389,25 @@ namespace AMDiS {
class fgmres_type
{
int restart;
ORTHOGONALIZATION orthogonalization;
public:
fgmres_type(std::string name) : restart(30)
fgmres_type(std::string name) : restart(30), orthogonalization(GRAM_SCHMIDT)
{
Parameters::get(name + "->restart", restart);
Parameters::get(name + "->orthogonalization", orthogonalization);
}
template < class LinOp, class X, class B, class L, class R, class I >
int operator()(const LinOp& A, X& x, const B& b, const L& l, const R& r, I& iter)
{ return itl::fgmres(A, x, b, l, r, iter, restart); }
{
switch (orthogonalization) {
default:
case GRAM_SCHMIDT:
return itl::fgmres(A, x, b, l, r, iter, restart); break;
case HOUSEHOLDER:
return itl::fgmres_householder(A, x, b, r, iter, restart); break;
}
}
};
typedef ITL_Solver< fgmres_type > FGMResSolver;
......
......@@ -67,15 +67,15 @@ namespace AMDiS {
KrylovPreconditioner(std::string name) : fullMatrix(nullptr), solver(nullptr), runner(nullptr)
{
#if defined HAVE_PARALLEL_PETSC
#if defined HAVE_PARALLEL_PETSC
string backend("p_petsc");
#elif defined HAVE_PARALLEL_MTL
#elif defined HAVE_PARALLEL_MTL
string backend("p_mtl");
#elif defined HAVE_PETSC
#elif defined HAVE_PETSC
string backend("petsc");
#else
#else
string backend("mtl");
#endif
#endif
// === read backend-name ===
string initFileStr = name + "->solver";
......
......@@ -240,7 +240,11 @@ namespace AMDiS {
/// fill PETSc vector from DOFVector
inline void operator<<(Vec& petscVec, /* const */SystemVector& vec)
{
const VecType vecType; // TODO: type changed from petsc 3.2 to petsc 3.3
#ifdef VecType // PETSc uses MACROS instead of typedefs in Versions 3.3x
const VecType vecType;
#else
VecType vecType;
#endif
VecGetType(petscVec, &vecType);
if (strcmp(vecType, VECNEST) == 0) {
......@@ -265,7 +269,11 @@ namespace AMDiS {
/// fill SystemVector from PETSc vector
inline void operator>>(const Vec& petscVec, SystemVector& vec)
{
#ifdef VecType // PETSc uses MACROS instead of typedefs in Versions 3.3x
const VecType vecType;
#else
VecType vecType;
#endif
VecGetType(petscVec, &vecType);
if (strcmp(vecType, VECNEST) == 0) {
......
/******************************************************************************
*
* AMDiS - Adaptive multidimensional simulations
*
* Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
* Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
*
* Authors:
* Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
*
* This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
* WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
*
* This file is part of AMDiS
*
* See also license.opensource.txt in the distribution.
*
******************************************************************************/
// Written by Simon Praetorius (adopted from previous implementation)
#ifndef ITL_FGMRES_HOUSEHOLDER_INCLUDE
#define ITL_FGMRES_HOUSEHOLDER_INCLUDE
#include <algorithm>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/matrix/multi_vector.hpp>
#include <boost/numeric/mtl/operation/givens.hpp>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/irange.hpp>
#include "solver/itl/details.hpp"
namespace itl {
/// Flexible Generalized Minimal Residual method (without restart) using householder othogonalization
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
regardless on whether the termination criterion is reached or not. **/
template < typename Matrix, typename Vector, typename RightPreconditioner, typename Iteration >
int fgmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
RightPreconditioner &R, Iteration& iter)
{
using mtl::irange; using std::abs; using math::reciprocal; using mtl::iall; using mtl::imax;
using mtl::signum; using mtl::vector::dot; using mtl::conj;
typedef typename mtl::Collection<Vector>::value_type Scalar;
typedef typename mtl::Collection<Vector>::size_type Size;
if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16;
Scalar rho, bnrm2, temp, alpha, beta;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector w(solve(R, b)), r(b - A*x);
mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1);
mtl::matrix::multi_vector<Vector> Z(Vector(resource(x), zero), kmax+1);
mtl::vector::dense_vector<Scalar> sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
mtl::matrix::dense2D<Scalar> H(kmax, kmax);
bnrm2 = two_norm(w);
if (bnrm2 < dbl_tol)
bnrm2 = 1.0;
// TODO: if rhs=0 => solution=0 ??
temp = two_norm(r); // norm of preconditioned residual
rho = temp * reciprocal(bnrm2);
if (iter.finished(rho)) // initial guess is good enough solution
return iter;
// u = r + sing(r(0))*||r||*e0
beta = signum(r[0])*temp;
w = r;
w[0] += beta;
w *= reciprocal(two_norm(w));
V.vector(0) = w;
H = zero;
s[0] = -beta;
// GMRES iteration
for (k= 0; k < kmax && !iter.finished(rho); ++k, ++iter) {
w = (-2.0 * V.vector(k)[k])*V.vector(k);
w[k] += 1.0;
// v := P_0*...*P_{k-2}*(P_{k-1} * e_k)
for (Size i= k; i > 0; i--) {
temp = 2.0 * dot(V.vector(i-1), w);
w -= temp * V.vector(i-1);
}
temp = two_norm(w);
if (temp == zero)
return iter.fail(2, "GMRES: breakdown");
// Explicitly normalize v to reduce the effects of round-off.
Z.vector(k) = solve(R, w);
w = A * Z.vector(k);
// P_{k-1}*...*P_0*Av
for (Size i = 0; i <= k; i++) {
temp = 2.0 * dot(V.vector(i), w);
w -= temp * V.vector(i);
}
temp = two_norm(w);
if (temp == zero)
return iter.fail(3, "GMRES: breakdown");
irange range_to_end(k+1,imax);
set_to_zero(V.vector(k+1));
V.vector(k+1)[range_to_end] = w[range_to_end];
alpha = two_norm(V.vector(k+1));
if (alpha != 0.0) {
alpha *= signum(w[k+1]);
V.vector(k+1)[k+1] += alpha;
V.vector(k+1) *= reciprocal(two_norm(V.vector(k+1)));
w[k+1] = -alpha;
}
for (Size i= 0; i < k; i++) {
temp = conj(cs[i])*w[i] + conj(sn[i])*w[i+1];
w[i+1] =- sn[i]*w[i] + cs[i]*w[i+1];
w[i] = temp;
}
details::rotmat(w[k], w[k+1], cs[k], sn[k]);
s[k+1] = -sn[k]*s[k];
s[k] = conj(cs[k])*s[k];
w[k] = cs[k]*w[k] + sn[k]*w[k+1];
w[k+1] = 0.0;
irange range(num_rows(H));
H[iall][k] = w[range];
rho = std::abs(s[k+1]) / bnrm2;
}
// reduce k, to get regular matrix
// while (k > 0 && std::abs(s[k-1]) <= iter.atol()) k--;
// iteration is finished -> compute x: solve H*y=s as far as rank of H allows
irange range(k);
for (; !range.empty(); --range) {
try {
y[range] = upper_trisolve(H[range][range], s[range]);
} catch (mtl::matrix_singular) { continue; } // if singular then try with sub-matrix
break;
}
if (range.finish() < k)
std::cerr << "GMRES orhogonalized with " << k << " vectors but matrix singular, can only use "
<< range.finish() << " vectors!\n";
if (range.empty())
return iter.fail(3, "GMRES did not find any direction to correct x");
x += Z.vector(range) * y[range];
r = b - A*x;
return iter.terminate(r);
}
/// Felxible Generalized Minimal Residual method with restart
template < typename Matrix, typename Vector,
typename RightPreconditioner, typename Iteration >
int fgmres_householder(const Matrix &A, Vector &x, const Vector &b,
RightPreconditioner &R,
Iteration& iter, typename mtl::Collection<Vector>::size_type restart)
{
do {
Iteration inner(iter);
inner.set_max_iterations(std::min(int(iter.iterations()+restart), iter.max_iterations()));
inner.suppress_resume(true);
fgmres_householder_full(A, x, b, R, inner);
iter.update_progress(inner);
} while (!iter.finished());
return iter;
}
} // namespace itl
#endif // ITL_FGMRES_HOUSEHOLDER_INCLUDE
......@@ -62,10 +62,11 @@ int gmres2_full(const Matrix &A, Vector &x, const Vector &b,
bnrm2 = two_norm(b);
if (bnrm2 < dbl_tol)
bnrm2 = 1.0;
// TODO: if rhs=0 => solution=0 ??
temp = two_norm(r);
temp = two_norm(r); // norm of preconditioned residual
rho = temp * reciprocal(bnrm2);
if (iter.finished(rho))
if (iter.finished(rho)) // initial guess is good enough solution
return iter;
V.vector(0) = r * reciprocal(temp);
......@@ -78,7 +79,7 @@ int gmres2_full(const Matrix &A, Vector &x, const Vector &b,
r = A * Vector(solve(R, V.vector(k)));
w = solve(L, r);
temp = two_norm(w);
// orth(V, w, false);
// modified Gram Schmidt method
for (Size i= 0; i < k+1; i++) {
......
/******************************************************************************
*
* AMDiS - Adaptive multidimensional simulations
*
* Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
* Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
*
* Authors:
* Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
*
* This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
* WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
*
* This file is part of AMDiS
*
* See also license.opensource.txt in the distribution.
*
******************************************************************************/
// Written by Simon Praetorius (adopted from previous implementation)
#ifndef ITL_GMRES_HOUSEHOLDER_INCLUDE
#define ITL_GMRES_HOUSEHOLDER_INCLUDE
#include <algorithm>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/matrix/multi_vector.hpp>
#include <boost/numeric/mtl/operation/givens.hpp>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/irange.hpp>
#include "solver/itl/details.hpp"
namespace itl {
/// Generalized Minimal Residual method (without restart) using householder othogonalization.
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
regardless on whether the termination criterion is reached or not. **/
template < typename Matrix, typename Vector, typename LeftPreconditioner, typename Iteration >
int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
LeftPreconditioner &L, Iteration& iter)
{
using mtl::irange; using std::abs; using math::reciprocal; using mtl::iall; using mtl::imax;
using mtl::signum; using mtl::vector::dot; using mtl::conj;
typedef typename mtl::Collection<Vector>::value_type Scalar;
typedef typename mtl::Collection<Vector>::size_type Size;
if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16; // kappa in [1.25, 100]
Scalar rho, bnrm2, temp, alpha, beta;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector w(b - A *x), r(solve(L,w));
mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1);
mtl::vector::dense_vector<Scalar> sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
mtl::matrix::dense2D<Scalar> H(kmax, kmax);
bnrm2 = two_norm(b);
if (bnrm2 < dbl_tol)
bnrm2 = 1.0;
// TODO: if rhs=0 => solution=0 ??
temp = two_norm(r); // norm of preconditioned residual
rho = temp * reciprocal(bnrm2);
if (iter.finished(rho)) // initial guess is good enough solution
return iter;
// u = r + sing(r(0))*||r||*e0
beta = signum(r[0])*temp;
w = r;
w[0] += beta;
w *= reciprocal(two_norm(w));
V.vector(0) = w;
H = zero;
s[0] = -beta;
// GMRES iteration
for (k= 0; k < kmax && !iter.finished(rho); ++k, ++iter) {
w = (-2.0 * V.vector(k)[k])*V.vector(k);
w[k] += 1.0;
// v := P_0*...*P_{k-2}*(P_{k-1} * e_k)
for (Size i= k; i > 0; i--) {
temp = 2.0 * dot(V.vector(i-1), w);
w -= temp * V.vector(i-1);
}
temp = two_norm(w);
if (temp == zero)
return iter.fail(2, "GMRES: breakdown");
// Explicitly normalize v to reduce the effects of round-off.
w *= reciprocal(temp);
w = solve(L, Vector(A*w));
// P_{k-1}*...*P_0*Av
for (Size i = 0; i <= k; i++) {
temp = 2.0 * dot(V.vector(i), w);
w -= temp * V.vector(i);
}
temp = two_norm(w);
if (temp == zero)
return iter.fail(3, "GMRES: breakdown");
irange range_to_end(k+1,imax);
set_to_zero(V.vector(k+1));
V.vector(k+1)[range_to_end] = w[range_to_end];
alpha = two_norm(V.vector(k+1));
if (alpha != 0.0) {
alpha *= signum(w[k+1]);
V.vector(k+1)[k+1] += alpha;
V.vector(k+1) *= reciprocal(two_norm(V.vector(k+1)));
w[k+1]