// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) XXXX-2008 - INRIA
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt


// [x, flag, resNorm, iter, resVec] = gmres( A, b, x, M, restrt, max_it, tol )
//
// GMRES solves the linear system Ax=b
// using the Generalized Minimal RESidual ( GMRES ) method with restarts .
//
// input   A        REAL nonsymmetric positive definite matrix or function
//         x        REAL initial guess vector
//         b        REAL right hand side vector
//         M        REAL preconditioner matrix or function
//         restrt   INTEGER number of iterations between restarts
//         max_it   INTEGER maximum number of iterations
//         tol      REAL error tolerance
//
// output  x        REAL solution vector
//         flag     INTEGER: 0 = solution found to tolerance
//                           1 = no convergence given max_it
//         resNorm      REAL final residual norm
//         iter     INTEGER number of iterations performed
//         resVec      REAL residual vector

//     Details of this algorithm are described in
//
//     "Templates for the Solution of Linear Systems: Building Blocks
//     for Iterative Methods",
//     Barrett, Berry, Chan, Demmel, Donato, Dongarra, Eijkhout,
//     Pozo, Romine, and Van der Vorst, SIAM Publications, 1993
//     (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
//
//     "Iterative Methods for Sparse Linear Systems, Second Edition"
//     Saad, SIAM Publications, 2003
//     (ftp ftp.cs.umn.edu; cd dept/users/saad/PS; get all_ps.zip).

function [x, flag, resNorm, iter, resVec] = gmres(A, varargin)

// -----------------------
// Parsing input arguments
// -----------------------

[lhs,rhs]=argn(0);
if ( rhs < 2 ),
  error(msprintf(gettext("%s: Wrong number of input argument: At least %d expected.\n"),"gmres",2));
end

// Parsing the matrix A et the right hand side vector b
select type(A)
case 1 then
  matrixType = 1;
case 5 then
  matrixType = 1;
case 13 then
  matrixType = 0;
end
// If A is a matrix (full or sparse)
if (matrixType == 1),
  if (size(A,1) ~= size(A,2)),
    error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"gmres",1));
  end
end
b=varargin(1);
if (size(b,2) ~= 1),
  error(msprintf(gettext("%s: Wrong size for input argument #%d: Column vector expected.\n"),"gmres",2));
end
if (matrixType==1),
  if (size(b,1) ~= size(A,1)),
  error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"gmres",2,1));
  end
end

// Number of iterations between restarts
if (rhs >= 3),
  restrt=varargin(2);
  if (size(restrt) ~= [1 1]),
    error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"gmres",3));
  end
else
  restrt=20;
end

// Error tolerance tol
if (rhs >= 4),
  tol=varargin(3);
  if (size(tol) ~= [1 1]);
    error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"gmres",4));
  end
else
  tol = 1e-6;
end

// Maximum number of iterations max_it
if (rhs >= 5),
  max_it=varargin(4);
  if (size(max_it) ~= [1 1]),
    error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"),"gmres",5));
  end
else
  max_it=size(b,1);
end

// Parsing of the preconditioner matrix M
if (rhs >= 6),
  M = varargin(5);
  select type(M)
  case 1 then
    precondType = 1;
  case 5 then
    precondType = 1;
  case 13 then
    precondType = 0;
  end
  if (precondType == 1),
    if (size(M,1) ~= size(M,2)),
      error(msprintf(gettext("%s: Wrong size for input argument #%d: Square matrix expected.\n"),"gmres",4));
    end
    if (size(M,1) == 0),
      precondType = 2; // no preconditionning
    elseif ( size(M,1) ~= size(b,1) ),
      error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"gmres",4,2));
    end
  end
  if (precondType == 0),
    M = varargin(3);
  end
else
  precondType = 2; // no preconditionning
end

// Parsing of the initial vector x
if (rhs >= 7),
  x=varargin(6);
  if (size(x,2) ~= 1),
    error(msprintf(gettext("%s: Wrong size for input argument #%d: Column vector expected.\n"),"gmres",3));
  end
  if ( size(x,1) ~= size(b,1) ),
    error(msprintf(gettext("%s: Wrong size for input argument #%d: Same size as input argument #%d expected.\n"),"gmres",3,2));
  end
else
  x=zeros(b);
end

if (rhs > 7),
  error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"gmres",2,7));
end

// ------------
// Computations
// ------------

   j = 0;
   flag = 0;
   it2 = 0;

   bnrm2 = norm(b);
   if (bnrm2 == 0.0),
     x = zeros(b);
     resNorm = 0;
     iter = 0;
     resVec = resNorm;
     flag = 0;
     return
   end

   // r = M \ ( b-A*x );
   if (matrixType == 1),
     r = b - A*x;
   else
     r = b - A(x);
   end
   if (precondType == 1),
     r = M \ r;
   elseif (precondType == 0),
     r = M(r);
   end
   resNorm = norm(r)/bnrm2;
   resVec = resNorm;
   if (resNorm < tol),
     iter=0;
     return;
   end

   n = size(b,1);
   m = restrt;
   V(1:n,1:m+1) = zeros(n,m+1);
   H(1:m+1,1:m) = zeros(m+1,m);
   cs(1:m) = zeros(m,1);
   sn(1:m) = zeros(m,1);
   e1    = zeros(n,1);
   e1(1) = 1.0;

   for j = 1:max_it
     // r = M \ ( b-A*x );
     if (matrixType == 1),
       r = b - A*x;
     else
       r = b - A(x);
     end
     if (precondType == 1),
       r = M \ r;
     elseif (precondType == 0),
       r = M(r);
     end

     V(:,1) = r / norm( r );
     s = norm( r )*e1;
     for i = 1:m      // construct orthonormal
       it2 = it2 + 1; // basis using Gram-Schmidt
       // w = M \ (A*V(:,i));
       if (matrixType == 1),
         w = A*V(:,i);
       else
         w = A(V(:,i));
       end
       if (precondType == 1),
         w = M \ w;
       elseif (precondType == 0),
         w = M(w);
       end

       for k = 1:i
         H(k,i)= w'*V(:,k);
         w = w - H(k,i)*V(:,k);
       end
       H(i+1,i) = norm( w );
       V(:,i+1) = w / H(i+1,i);
       for k = 1:i-1 // apply Givens rotation
         temp     =  cs(k)*H(k,i) + sn(k)*H(k+1,i);
         H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
         H(k,i)   = temp;
       end
        // form i-th rotation matrix
       [tp1,tp2] = rotmat( H(i,i), H(i+1,i) );
       cs(i)  = tp1;
       sn(i)  = tp2;
       temp   = cs(i)*s(i);
       s(i+1) = -sn(i)*s(i);
       s(i)   = temp;
       H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
       H(i+1,i) = 0.0;
       resNorm  = abs(s(i+1)) / bnrm2;
       resVec = [resVec;resNorm];
       if ( resNorm <= tol ),
         y = H(1:i,1:i) \ s(1:i);
         x = x + V(:,1:i)*y;
         break;
       end
     end
     if (resNorm <= tol),
       iter = j-1+it2;
       break;
     end
     y = H(1:m,1:m) \ s(1:m);
     // update approximation
     x = x + V(:,1:m)*y;
     // r = M \ ( b-A*x )
     if (matrixType == 1),
       r = b - A*x;
     else
       r = b - A(x);
     end
     if (precondType == 1),
       r = M \ r;
     elseif (precondType == 0),
       r = M(r);
     end
     s(j+1) = norm(r);
     resNorm = s(j+1) / bnrm2;
     resVec = [resVec; resNorm];

     if ( resNorm <= tol ),
       iter = j+it2;
       break;
     end
     if ( j== max_it ),
       iter=j+it2;
     end
   end
   if ( resNorm > tol ),
     flag = 1;
     if (lhs < 2),
       warning(msprintf(gettext("%s: Did not converge.\n"),"gmres"));
     end
   end
endfunction //GMRES
 

//
// Compute the Givens rotation matrix parameters for a and b.
//
function [ c, s ] = rotmat( a, b )
if ( b == 0.0 ),
  c = 1.0;
  s = 0.0;
elseif ( abs(b) > abs(a) ),
  temp = a / b;
  s = 1.0 / sqrt( 1.0 + temp^2 );
  c = temp * s;
else
  temp = b / a;
  c = 1.0 / sqrt( 1.0 + temp^2 );
  s = temp * c;
end
endfunction //rotmat