function [x,tnorm,flag]=CMRH(A,b,tol,maxiter,x0,iprint)
%
% function [x,tnorm,flag] = (A,b,tol,maxiter,x0,iprint)
%
% This function solves the linear system Ax=b using the CMRH method
%
% Input :   
%           A       : the n times n matrix
%           b       : the right hand side vector
%           tol     : the stopping criterion 
%           maxiter : the maximum number of iterations
%           x0      : the initial approximation
%           iprint  : iprint = 1 to print some components of the 
%                     approximated solution and the norm of the residual
%                     in the Matlab Command Window.
%                   : iprint = 0 if we do not want to print the components
%                     and the norm of the residuals in the matlab Command 
%                     Window. 
%
% Output :  
%           x       : the solution 
%           tnorm   : the pseudo residual norm array
%           flag    : the stopping criterion uses
%                   if flag = 0  : we use tol
%                   if flag = 1  : we use maxiter 
%

n=length(A(:,1));

if (nargin < 5)
    x = x0;
else
    x = zeros(n,1);
end
if (nargin < 4)
    maxiter = n;
end
if (nargin < 3)
    tol = 1.d-7;
end

rs = zeros(n,1); s = zeros(n,1); c = zeros(n,1);
p = 1:n; b = b-A*x;
[~,indice] = max(abs(b(1:n)));
beta = b(indice); D_inv = 1/beta; rs(1) = beta;
dnorm = abs(beta); tnorm = dnorm;
if ((indice < 1)||(indice > 1))
    temp = p(1); p(1) = p(indice); p(indice) = temp;
    temp = b(1); b(1) = b(indice); b(indice) = temp;
    Ai = A(:,indice);
    A(:,indice) = A(:,1); A(:,1) = Ai; 
    Aip = A(indice,:); A(indice,:)=A(1,:); A(1,:) = Aip;
end
b = D_inv.*b;

k = 0;
while((dnorm > tol)&&(k < maxiter))
    k = k+1;
    n_k = n-k;
    k1 = k+1;
    if (iprint==1)
       fprintf(' %3.0f %6.4e \n',k-1,dnorm);
    end
    work = A(:,k);
    work = A(:,k1:n)*b(k1:n)+work;
    A(k1:n,k) = b(k1:n);
    for j = 1:k
        A(j,k) = work(j);
        work(j) = 0.d0;
        j1 = j+1;
        work(j1:n) = -A(j,k).*A(j1:n,j)+work(j1:n);
    end
    if (k < n)
        [~,indice] = max(abs(work(k1:n)));
        indice = k+indice;
        H_k1 = work(indice);
        if ((indice < k1)||(indice > k1))
            temp = p(k1); p(k1) = p(indice); p(indice) = temp;
            temp = work(k1); work(k1) = work(indice); work(indice) = temp;
            Ai = A(:,indice); A(:,indice) = A(:,k1);
            A(:,k1) = Ai; Aip = A(indice,:); A(indice,:) = A(k1,:);
            A(k1,:) = Aip;
        end
        if ((H_k1 < 0)||(H_k1 > 0))
            D_inv = 1/H_k1; b = work; b = D_inv.*b;
        end
    else
        H_k1 = 0;
        b = zeros(n,1);
    end
    
    % solving the least squares problem    
    if (k > 1)
        for i = 2:k
            i_1 = i-1;
            dtemp = A(i_1,k);
            A(i_1,k) = c(i_1)*dtemp+s(i_1)*A(i,k);
            A(i,k) = -s(i_1)*dtemp+c(i_1)*A(i,k);
        end
    end
    if (k <= n)
        dgam = sqrt(A(k,k)^2+H_k1^2);
        c(k) = abs(A(k,k))/dgam;
        s(k) = (H_k1*abs(A(k,k)))/(A(k,k)*dgam);
        rs(k1) = -s(k)*rs(k);
        rs(k) = c(k)*rs(k);
        A(k,k) = c(k)*A(k,k)+s(k)*H_k1;
        dnorm = abs(rs(k1));
    end
    tnorm = [tnorm dnorm];
end
if (k==maxiter)
    flag=1;
else
    flag=0;
end
if (iprint==1)
    fprintf(' %3.0f %6.4e \n',k,dnorm);
end

% computing d the solution of the least square problem
for j = k:-1:1
    rs(j) = rs(j)/A(j,j);
    rs(1:j-1) = rs(1:j-1)-rs(j).*A(1:j-1,j);
end
work(1:k) = rs(1:k);
for j = k:-1:1
    rs(j+1:k) = rs(j+1:k)+rs(j).*A(j+1:k,j);
end
c = A(k1:n,1:k)*work(1:k);
x(1:k) = rs(1:k)+x(1:k);
x(k1:n) = c(1:n_k)+x(k1:n);

% Reorder the components of x to get the solution
xx(p(1:n)) = x(1:n);
x = xx';