computation/solver/iterative/Gmres.cxx

// 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
{


  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