Commit 836978e9 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

removed all trailing whitespaces in all git files

parent 5f0ca3b5
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
// Written by Jan Bos
......@@ -31,12 +31,12 @@
namespace itl {
/// Bi-Conjugate Gradient Stabilized(ell)
template < typename LinearOperator, typename Vector,
typename LeftPreconditioner, typename RightPreconditioner,
/// Bi-Conjugate Gradient Stabilized(ell)
template < typename LinearOperator, typename Vector,
typename LeftPreconditioner, typename RightPreconditioner,
typename Iteration >
int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
const LeftPreconditioner &L, const RightPreconditioner &R,
const LeftPreconditioner &L, const RightPreconditioner &R,
Iteration& iter, size_t l)
{
mtl::vampir_trace<7006> tracer;
......@@ -50,7 +50,7 @@ int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
Vector x0(resource(x)), y(resource(x));
mtl::vector::dense_vector<Vector> r_hat(l+1,Vector(resource(x))), u_hat(l+1,Vector(resource(x)));
// shift problem
// shift problem
x0= zero;
r_hat[0]= b;
if (two_norm(x) != zero) {
......@@ -64,7 +64,7 @@ int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
r_hat[0]= y;
u_hat[0]= zero;
Scalar rho_0(one), rho_1(zero), alpha(zero), Gamma(zero), beta(zero), omega(one);
Scalar rho_0(one), rho_1(zero), alpha(zero), Gamma(zero), beta(zero), omega(one);
mtl::matrix::dense2D<Scalar> tau(l+1, l+1);
mtl::vector::dense_vector<Scalar> sigma(l+1), gamma(l+1), gamma_a(l+1), gamma_aa(l+1);
......@@ -73,20 +73,20 @@ int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
rho_0= -omega * rho_0;
for (Size j= 0; j < l; ++j) {
rho_1= dot(r0_tilde, r_hat[j]);
rho_1= dot(r0_tilde, r_hat[j]);
beta= alpha * rho_1/rho_0; rho_0= rho_1;
for (Size i= 0; i <= j; ++i)
u_hat[i]= r_hat[i] - beta * u_hat[i];
y= A * Vector(solve(R, u_hat[j]));
u_hat[j+1]= solve(L, y);
Gamma= dot(r0_tilde, u_hat[j+1]);
Gamma= dot(r0_tilde, u_hat[j+1]);
alpha= rho_0 / Gamma;
for (Size i= 0; i <= j; ++i)
r_hat[i]-= alpha * u_hat[i+1];
if (iter.finished(r_hat[j])) {
x= solve(R, x);
x+= x0;
......@@ -94,7 +94,7 @@ int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
}
r_hat[j+1]= solve(R, r_hat[j]);
y= A * r_hat[j+1];
y= A * r_hat[j+1];
r_hat[j+1]= solve(L, y);
x+= alpha * u_hat[0];
}
......@@ -103,13 +103,13 @@ int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
irange i1m(1, imax);
mtl::vector::dense_vector<Vector> r_hat_tail(r_hat[i1m]);
tau[i1m][i1m]= orthogonalize_factors(r_hat_tail);
for (Size j= 1; j <= l; ++j)
for (Size j= 1; j <= l; ++j)
gamma_a[j]= dot(r_hat[j], r_hat[0]) / tau[j][j];
gamma[l]= gamma_a[l]; omega= gamma[l];
if (omega == zero) return iter.fail(3, "bicg breakdown #2");
// is this something like a tri-solve?
// is this something like a tri-solve?
for (Size j= l-1; j > 0; --j) {
Scalar sum= zero;
for (Size i=j+1;i<=l;++i)
......@@ -134,7 +134,7 @@ int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
/// Solver class for BiCGStab(ell) method; right preconditioner ignored (prints warning if not identity)
template < typename LinearOperator, typename Preconditioner,
template < typename LinearOperator, typename Preconditioner,
typename RightPreconditioner>
class bicgstab_ell_solver
{
......@@ -146,7 +146,7 @@ class bicgstab_ell_solver
bicgstab_ell_solver(const LinearOperator& A, size_t l, const Preconditioner& L) : A(A), l(l), L(L), R(A) {}
/// Construct solver from a linear operator and left preconditioner
bicgstab_ell_solver(const LinearOperator& A, size_t l, const Preconditioner& L, const RightPreconditioner& R)
bicgstab_ell_solver(const LinearOperator& A, size_t l, const Preconditioner& L, const RightPreconditioner& R)
: A(A), l(l), L(L), R(R) {}
/// Solve linear system approximately as specified by \p iter
......@@ -163,7 +163,7 @@ class bicgstab_ell_solver
itl::basic_iteration<double> iter(x, 1, 0, 0);
return solve(b, x, iter);
}
private:
const LinearOperator& A;
size_t l;
......
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_CG_INCLUDE
......@@ -33,9 +33,9 @@
namespace itl {
/// Conjugate Gradients without preconditioning
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Iteration& iter)
{
mtl::vampir_trace<7001> tracer;
......@@ -46,20 +46,20 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Scalar rho(0), rho_1(0), alpha(0), alpha_1(0);
Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
r = b - A*x;
rho = dot(r, r);
while (! iter.finished(Real(sqrt(abs(rho))))) {
++iter;
if (iter.first())
p = r;
else
p = r + (rho / rho_1) * p;
else
p = r + (rho / rho_1) * p;
// q = A * p; alpha = rho / dot(p, q);
(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
alpha= rho / alpha_1;
x += alpha * p;
rho_1 = rho;
(lazy(r) -= alpha * q) || (lazy(rho) = lazy_unary_dot(r));
......@@ -69,9 +69,9 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
}
/// Conjugate Gradients
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Preconditioner, typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& L, Iteration& iter)
{
using pc::is_identity;
......@@ -86,7 +86,7 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Scalar rho(0), rho_1(0), rr, alpha(0), alpha_1;
Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
r = b - A*x;
rr = dot(r, r);
while (! iter.finished(Real(sqrt(abs(rr))))) {
......@@ -95,12 +95,12 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
if (iter.first())
p = z;
else
else
p = z + (rho / rho_1) * p;
(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
alpha= rho / alpha_1;
x += alpha * p;
rho_1 = rho;
(lazy(r) -= alpha * q) || (lazy(rr) = lazy_unary_dot(r));
......@@ -109,29 +109,29 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
}
/// Conjugate Gradients with ignored right preconditioner to unify interface
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Preconditioner, typename RightPreconditioner, typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& L, const RightPreconditioner&, Iteration& iter)
{
return cg(A, x, b, L, iter);
}
/// Solver class for CG method; right preconditioner ignored (prints warning if not identity)
template < typename LinearOperator, typename Preconditioner,
template < typename LinearOperator, typename Preconditioner,
typename RightPreconditioner>
class cg_solver
{
public:
/// Construct solver from a linear operator; generate (left) preconditioner from it
explicit cg_solver(const LinearOperator& A) : A(A), L(A)
explicit cg_solver(const LinearOperator& A) : A(A), L(A)
{
if (!pc::static_is_identity<RightPreconditioner>::value)
std::cerr << "Right Preconditioner ignored!" << std::endl;
}
/// Construct solver from a linear operator and (left) preconditioner
cg_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
cg_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
{
if (!pc::static_is_identity<RightPreconditioner>::value)
std::cerr << "Right Preconditioner ignored!" << std::endl;
......@@ -151,7 +151,7 @@ class cg_solver
itl::basic_iteration<double> iter(x, 1, 0, 0);
return solve(b, x, iter);
}
private:
const LinearOperator& A;
Preconditioner L;
......
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_CGS_INCLUDE
......@@ -21,7 +21,7 @@
namespace itl {
/// Conjugate Gradient Squared
template < typename LinearOperator, typename Vector,
template < typename LinearOperator, typename Vector,
typename Preconditioner, typename Iteration >
int cgs(const LinearOperator &A, Vector &x, const Vector &b,
const Preconditioner &M, Iteration& iter)
......@@ -52,7 +52,7 @@ int cgs(const LinearOperator &A, Vector &x, const Vector &b,
u+= q;
uhat= solve(M, u);
x+= alpha * uhat;
qhat= A * uhat;
r-= alpha * qhat;
......@@ -63,20 +63,20 @@ int cgs(const LinearOperator &A, Vector &x, const Vector &b,
}
/// Solver class for CGS method; right preconditioner ignored (prints warning if not identity)
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
typename RightPreconditioner= pc::identity<LinearOperator> >
class cgs_solver
{
public:
/// Construct solver from a linear operator; generate (left) preconditioner from it
explicit cgs_solver(const LinearOperator& A) : A(A), L(A)
explicit cgs_solver(const LinearOperator& A) : A(A), L(A)
{
if (!pc::static_is_identity<RightPreconditioner>::value)
std::cerr << "Right Preconditioner ignored!" << std::endl;
}
/// Construct solver from a linear operator and (left) preconditioner
cgs_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
cgs_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
{
if (!pc::static_is_identity<RightPreconditioner>::value)
std::cerr << "Right Preconditioner ignored!" << std::endl;
......@@ -96,7 +96,7 @@ class cgs_solver
itl::basic_iteration<double> iter(x, 1, 0, 0);
return solve(b, x, iter);
}
private:
const LinearOperator& A;
Preconditioner L;
......
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_FSM_INCLUDE
......@@ -19,7 +19,7 @@ namespace itl {
/// Folded spectrum method
/** Computed and named as in http://en.wikipedia.org/wiki/Folded_spectrum_method **/
template < typename LinearOperator, typename VectorSpace, typename EigenValue,
template < typename LinearOperator, typename VectorSpace, typename EigenValue,
typename Damping, typename Iteration >
int fsm(const LinearOperator& H, VectorSpace& phi, EigenValue eps, Damping alpha, Iteration& iter)
{
......
......@@ -29,7 +29,7 @@
namespace itl {
/// Generalized Minimal Residual method (without restart)
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
/** 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 RightPreconditioner, typename Iteration >
int gmres_full(const Matrix &A, Vector &x, const Vector &b,
......@@ -45,8 +45,8 @@ int gmres_full(const Matrix &A, Vector &x, const Vector &b,
Scalar rho, nu, hr;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector r0(b - A *x), r(solve(L,r0)), va(resource(x)), va0(resource(x)), va00(resource(x));
mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1);
mtl::vector::dense_vector<Scalar> s(kmax+1, zero), c(kmax+1, zero), g(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1);
mtl::vector::dense_vector<Scalar> s(kmax+1, zero), c(kmax+1, zero), g(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
mtl::matrix::dense2D<Scalar> H(kmax+1, kmax); // dito
H= 0;
......@@ -60,7 +60,7 @@ int gmres_full(const Matrix &A, Vector &x, const Vector &b,
for (k= 0; k < kmax ; ++k, ++iter) {
va0= A * Vector(solve(R, V.vector(k)));
V.vector(k+1)= va= solve(L,va0);
// orth(V, V[k+1], false);
// orth(V, V[k+1], false);
// modified Gram Schmidt method
for (Size j= 0; j < k+1; j++) {
H[j][k]= dot(V.vector(j), V.vector(k+1));
......@@ -74,13 +74,13 @@ int gmres_full(const Matrix &A, Vector &x, const Vector &b,
V.vector(k+1)-= hr * V.vector(j);
}
H[k+1][k]= two_norm(V.vector(k+1));
if (H[k+1][k] != zero) // watch for breakdown
if (H[k+1][k] != zero) // watch for breakdown
V.vector(k+1)*= 1. / H[k+1][k];
// k Given's rotations
for(Size i= 0; i < k; i++)
mtl::matrix::givens<mtl::matrix::dense2D<Scalar> >(H, H[i][k-1], H[i+1][k-1]).trafo(i);
nu= sqrt(H[k][k]*H[k][k]+H[k+1][k]*H[k+1][k]);
if(nu != zero){
c[k]= H[k][k]/nu;
......@@ -91,7 +91,7 @@ int gmres_full(const Matrix &A, Vector &x, const Vector &b,
}
rho= abs(g[k+1]);
}
//reduce k, to get regular matrix
while (k > 0 && abs(g[k-1]<= iter.atol())) k--;
......@@ -99,18 +99,18 @@ int gmres_full(const Matrix &A, Vector &x, const Vector &b,
irange range(k);
for (; !range.empty(); --range) {
try {
y[range]= lu_solve(H[range][range], g[range]);
y[range]= lu_solve(H[range][range], g[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 "
std::cerr << "GMRES orhogonalized with " << k << " vectors but matrix singular, can only use "
<< range.finish() << " vectors!\n";
if (range.empty())
return iter.fail(2, "GMRES did not find any direction to correct x");
x+= Vector(solve(R, Vector(V.vector(range)*y[range])));
r= b - A*x;
return iter.terminate(r);
}
......@@ -121,7 +121,7 @@ template < typename Matrix, typename Vector, typename LeftPreconditioner,
int gmres(const Matrix &A, Vector &x, const Vector &b,
LeftPreconditioner &L, 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()));
......@@ -134,21 +134,21 @@ int gmres(const Matrix &A, Vector &x, const Vector &b,
}
/// Solver class for GMRES; right preconditioner ignored (prints warning if not identity)
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
typename RightPreconditioner= pc::identity<LinearOperator> >
class gmres_solver
{
public:
/// Construct solver from a linear operator; generate (left) preconditioner from it
explicit gmres_solver(const LinearOperator& A, size_t restart= 8)
explicit gmres_solver(const LinearOperator& A, size_t restart= 8)
: A(A), restart(restart), L(A), R(A) {}
/// Construct solver from a linear operator and left preconditioner
gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L)
gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L)
: A(A), restart(restart), L(L), R(A) {}
/// Construct solver from a linear operator and left preconditioner
gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L, const RightPreconditioner& R)
gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L, const RightPreconditioner& R)
: A(A), restart(restart), L(L), R(R) {}
/// Solve linear system approximately as specified by \p iter
......@@ -165,7 +165,7 @@ class gmres_solver
itl::basic_iteration<double> iter(x, 1, 0, 0);
return solve(b, x, iter);
}
private:
const LinearOperator& A;
size_t restart;
......
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
// Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems.
// Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems.
// SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
#ifndef ITL_IDR_S_INCLUDE
......@@ -29,12 +29,12 @@
namespace itl {
/// Induced Dimension Reduction on s dimensions (IDR(s))
template < typename LinearOperator, typename Vector,
typename LeftPreconditioner, typename RightPreconditioner,
/// Induced Dimension Reduction on s dimensions (IDR(s))
template < typename LinearOperator, typename Vector,
typename LeftPreconditioner, typename RightPreconditioner,
typename Iteration >
int idr_s(const LinearOperator &A, Vector &x, const Vector &b,
const LeftPreconditioner &, const RightPreconditioner &,
const LeftPreconditioner &, const RightPreconditioner &,
Iteration& iter, size_t s)
{
mtl::vampir_trace<7010> tracer;
......@@ -49,10 +49,10 @@ int idr_s(const LinearOperator &A, Vector &x, const Vector &b,
Scalar omega(zero);
Vector x0(x), y(resource(x)), v(resource(x)), t(resource(x)), q(resource(x)), r(b - A * x);
mtl::matrix::multi_vector<Vector> dR(Vector(resource(x), zero), s), dX(Vector(resource(x), zero), s), P(Vector(resource(x), zero), s);
mtl::vector::dense_vector<Scalar> m(s), c(s), dm(s); // replicated in distributed solvers
mtl::vector::dense_vector<Scalar> m(s), c(s), dm(s); // replicated in distributed solvers
mtl::matrix::dense2D<Scalar> M(s, s); // dito
random(P);
random(P);
P.vector(0)= r;
orth(P);
......@@ -61,20 +61,20 @@ int idr_s(const LinearOperator &A, Vector &x, const Vector &b,
omega= dot(v, r) / dot(v, v);
dX.vector(k)= omega * r;
dR.vector(k)= -omega * v;
x+= dX.vector(k);
x+= dX.vector(k);
r+= dR.vector(k);
if ((++iter).finished(r)) return iter;
M[iall][k]= trans(P) * dR.vector(k);
M[iall][k]= trans(P) * dR.vector(k);
}
Size oldest= 0;
m= trans(P) * r;
while (! iter.finished(r)) {
for (size_t k= 0; k < s; k++) {
c= lu_solve(M, m);
q= dR * -c;
q= dR * -c;
v= r + q;
if (k == 0) {
t= A * v;
......@@ -101,7 +101,7 @@ int idr_s(const LinearOperator &A, Vector &x, const Vector &b,
}
/// Solver class for IDR(s) method; right preconditioner ignored (prints warning if not identity)
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
typename RightPreconditioner= pc::identity<LinearOperator> >
class idr_s_solver
{
......@@ -113,7 +113,7 @@ class idr_s_solver
idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L) : A(A), s(s), L(L), R(A) {}
/// Construct solver from a linear operator and left preconditioner
idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L, const RightPreconditioner& R)
idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L, const RightPreconditioner& R)
: A(A), s(s), L(L), R(R) {}
/// Solve linear system approximately as specified by \p iter
......@@ -130,7 +130,7 @@ class idr_s_solver
itl::basic_iteration<double> iter(x, 1, 0, 0);
return solve(b, x, iter);
}
private:
const LinearOperator& A;
size_t s;
......
......@@ -26,7 +26,7 @@ namespace itl {
/// Quasi-Minimal Residual method
template < typename Matrix, typename Vector,typename LeftPreconditioner,
typename RightPreconditioner, typename Iteration >
int qmr(const Matrix& A, Vector& x, const Vector& b, LeftPreconditioner& L,
int qmr(const Matrix& A, Vector& x, const Vector& b, LeftPreconditioner& L,
const RightPreconditioner& R, Iteration& iter)
{
mtl::vampir_trace<7008> tracer;
......@@ -36,7 +36,7 @@ int qmr(const Matrix& A, Vector& x, const Vector& b, LeftPreconditioner& L,
const Scalar zero= math::zero(Scalar()), one= math::one(Scalar());
Scalar beta, gamma(one), gamma_1, delta, eta(-one), ep(one), rho_1, theta(zero), theta_1;
Vector r(b - A * x), v_tld(r), y(solve(L, v_tld)), w_tld(r), z(adjoint_solve(R,w_tld)), v(resource(x)), w(resource(x)),
Vector r(b - A * x), v_tld(r), y(solve(L, v_tld)), w_tld(r), z(adjoint_solve(R,w_tld)), v(resource(x)), w(resource(x)),
y_tld(resource(x)), z_tld(resource(x)), p(resource(x)), q(resource(x)), p_tld(resource(x)), d(resource(x)), s(resource(x));