function x = cg_squared ( A, b, x )

%*****************************************************************************80
%
%% cg_squared() applies the CG squared algorithm.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 March 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(n,n): the matrix
%
%    real b(n): the right hand side.
%
%    real x(n): an initial solution estimate.
%
%  Output:
%
%    real x(n): the solution estimate after CG squared was applied.
%
  tol = 0.0000000001;

  normb = norm ( b );
  if ( normb == 0.0 )
    normb = 1.0;
  end

  r = b - A * x;
  resid = norm ( r ) / normb;
  if ( resid <= tol )
    return;
  end

  rtilde = r;
  rho = 0.0;

  n = size ( A, 1 );
  max_iter = n;

  for i = 1 : max_iter

    rho_old = rho;
    rho = rtilde' * r;

    if ( rho == 0.0 )
      return;
    end

    if ( i == 1 )
      u = r;
      p = u;
    else
      beta = rho / rho_old;
      u = r + beta * q;
      p = u + beta * ( q + beta * p );
    end

    phat = p;
    vhat = A * phat;
    alpha = rho / ( rtilde' * vhat );
    q = u - alpha * vhat;
    uhat = u + q;
    x = x + alpha * uhat;
    qhat = A * uhat;

    r = r - alpha * qhat;
    resid = norm ( r ) / normb;
    if ( resid < tol )
      return;
    end
  
  end

  return
end

