// Copyright (C) 2003-2009 Marc Duruflé
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 of the License, or (at your option)
// any later version.
//
// Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
// more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Seldon. If not, see http://www.gnu.org/licenses/.
#ifndef SELDON_FILE_ITERATIVE_GMRES_CXX
namespace Seldon
{
//! Solves a linear system by using Generalized Minimum Residual (GMRES)
/*!
Solves the unsymmetric linear system Ax = b using restarted GMRES.
return value of 0 indicates convergence within the
maximum number of iterations (determined by the iter object).
return value of 1 indicates a failure to converge.
See: Y. Saad and M. Schulter. GMRES: A generalized minimum residual
algorithm for solving nonsysmmetric linear systems, SIAM
J. Sci. Statist. Comp. 7(1986), pp, 856-869
\param[in] A Complex General Matrix
\param[in,out] x Vector on input it is the initial guess
on output it is the solution
\param[in] b Vector right hand side of the linear system
\param[in] M Right preconditioner
\param[in] outer Iteration parameters
*/
template <class Titer, class MatrixSparse, class Vector1, class Preconditioner>
int Gmres(MatrixSparse& A, Vector1& x, const Vector1& b,
Preconditioner& M, Iteration<Titer> & outer)
{
const int N = A.GetM();
if (N <= 0)
return 0;
typedef typename Vector1::value_type Complexe;
Complexe zero(0);
int m = outer.GetRestart();
// V is the array of orthogonal basis contructed
// from the Krylov subspace (v0,A*v0,A^2*v0,...,A^m*v0)
std::vector<Vector1> V(m+1, b);
// Upper triangular hessenberg matrix
// we don't store the sub-diagonal
// we apply rotations to eliminate this sub-diagonal
Matrix<Complexe, General, ColUpTriang> H(m+1,m+1); H.Fill(zero);
// s is the vector of residual norm for each inner iteration
// w is used in the Arnoldi algorithm
// u is a temporary vector which contains the product A*v(i)
// r is the residual
Vector1 w(b), r(b), u(b);
Vector<Complexe> s(m+1);
s.Fill(zero); w.Fill(zero); r.Fill(zero); u.Fill(zero);
for (int i = 0; i < m+1; i++)
V[i].Fill(zero);
Vector<Complexe> rotations_sin(m+1);
rotations_sin.Fill(zero);
Vector<Titer> rotations_cos(m+1);
rotations_cos.Fill(Titer(0));
// we compute residual
Copy(b, w);
if (!outer.IsInitGuess_Null())
MltAdd(Complexe(-1), A, x, Complexe(1), w);
else
x.Fill(zero);
// preconditioning
M.Solve(A, w, r);
Titer beta = Norm2(r);
// we initialize outer
int success_init = outer.Init(r);
if (success_init != 0)
return outer.ErrorCode();
// the coefficient H(m+1,m)
Complexe hi_ip1;
outer.SetNumberIteration(0);
// Loop until the stopping criteria are reached
while (! outer.Finished(beta))
{
// we normalize V(0) and we init s
Copy(r, V[0]);
Mlt(Complexe(Complexe(1)/beta), V[0]);
s.Fill(zero);
s(0) = beta;
int i = 0, k;
// we initialize the iter iteration
// m is the maximum number of inner iterations
Iteration<Titer> inner(outer);
inner.SetNumberIteration(outer.GetNumberIteration());
inner.SetMaxNumberIteration(outer.GetNumberIteration()+m);
do
{
// product matrix vector u=A*V(i)
Mlt(A, V[i], u);
// preconditioning
M.Solve(A, u, w);
// Arnoldi algorithm
for (k = 0; k <= i; k++)
{
// h_{k,i} = \bar{v(k)} w
H.Val(k, i) = DotProdConj(V[k], w);
Add(-H(k,i), V[k], w);
}
// we compute h(i+1,i)
hi_ip1 = Norm2(w);
Copy(w, V[i+1]);
// we normalize V(i+1)
if (hi_ip1 != zero)
Mlt(Complexe(1)/hi_ip1, V[i+1]);
// we apply precedent generated rotations
// to the last column we computed.
for (k = 0; k < i; k++)
ApplyRot(H.Val(k,i), H.Val(k+1,i),
rotations_cos(k), rotations_sin(k));
// we generate a new rotation Omega=[c s;-conj(s) c] in order to
// cancel h(i+1,i) and we store this rotation
if (hi_ip1 != zero)
{
GenRot(H.Val(i,i), hi_ip1,
rotations_cos(i), rotations_sin(i));
// After this call we must have hi_ip1=0
// GenRot must modify the entries H(i,i) hi_ip1
// we apply the rotation to the right hand side s
ApplyRot(s(i), s(i+1), rotations_cos(i), rotations_sin(i));
}
++inner, ++outer, ++i;
} while (! inner.Finished(abs(s(i))));
// Now we solve the triangular system H
H.Resize(i, i); s.Resize(i);
Solve(H, s); s.Resize(m+1);
// new iterate x = x + sum_0^{i-1} s(k)*V(k)
for (k = 0; k < i; k++)
Add(s(k), V[k], x);
// we compute the new residual
Copy(b, w);
MltAdd(Complexe(-1), A, x, Complexe(1), w);
M.Solve(A, w, r);
// residual norm
beta = Norm2(r);
}
return outer.ErrorCode();
}
} // end namespace
#define SELDON_FILE_ITERATIVE_GMRES_CXX
#endif