% %============================================================================%
% % Duke University                                                            %
% % K. P. Trofatter                                                            %
% % kpt2@duke.edu                                                              %
% %============================================================================%
% zgd_quad() - complex gradient descent on quadratic form
% 
% USAGE:
%   [x, rr] = zgd_quad(A, b, x0=zeros(size(b)), tol=1.0e-3, imax=numel(b))
%
% INPUT:
%   [n,n] complex  | A    | conjugate-symmetric positive-definite matrix
%   [1,1] function | A    | conjugate-symmetric positive-definite function
%   [n,1] complex  | b    | data vector
%   [n,1] complex  | x0   | initial guess vector
%   [1,1] double   | tol  | tolerance of residual magnitude squared
%   [1,1] double   | imax | maximum number of iterations
%
% OUTPUT:
%   [n,1] complex  | x    | approximate solution vector
%   [1,1] double   | rr   | residual magnitude squared


function [x, rr] = zgd_quad(A, b, x0, tol, imax)
    
    % default arguments
    if ~exist('x0', 'var') || isempty(x0)
        x0 = zeros(size(b));
    end
    if ~exist('tol', 'var') || isempty(tol)
        tol = 1.0e-3;
    end
    if ~exist('imax', 'var') || isempty(imax)
        imax = numel(b);
    end
    
    % build matrix multiplication function
    if isa(A, 'function_handle')
        Af = A;
    else
        Af = @(x) A * x;
    end
    
    % compute initial residual
    x = x0;
    r = b - Af(x);
    
    % iterate
    rr = zeros(1, imax);
    for i = 1 : imax
        
        % tolerance test
        rr(i) = r' * r;
        if rr(i) <= tol ^ 2
            break
        end
        
        % update
        Ar = Af(r);
        alpha = rr(i) / (r' * Ar);
        x = x + alpha * r;
        r = r - alpha * Ar;
        
    end
    
    % clip residual magnitude squared history
    rr = rr(1 : i);
    rr = sqrt(rr);
    
end


%==============================================================================%
%                                                                              %
%                                                                              %
%                                                                              %
%==============================================================================%