fixed_point_algorithm.m

function [x, defect, niterations] = fixed_point_algorithm(A, f, eval_nonlin, tol, maxiter, damping, verbose)
%function [x, defect, niterations] = fixed_point_algorithm(A, f, eval_nonlin, tol, maxiter, damping, verbose)
%
% used in stokes_detailed_simulation, stokes_rb_simulation,
% rbhdd_detailed_simulation
% A can be an LUPQ_Handle holding the l-u-p-q-factors if model is linear

% IM 28.05.2014

has_nonlinearity = ~isempty(eval_nonlin);

start_iteration = tic;
x = zeros(length(f), 1);
defect = f;
niterations = 0;

if nargin < 7
    verbose = 0;
end
if nargin < 6 || isempty(damping)
    damping = 1;
end
if nargin < 5 || isempty(maxiter)
    maxiter = 16;
end
if nargin < 4 || isempty(tol)
    tol = 1e-9;
end

if has_nonlinearity && verbose 
    disp('start fixed point iteration');
end

if has_nonlinearity
    A_lin = A;
end

while norm(defect) > tol
    
    x = x + damping * (A \ defect);
    
    if has_nonlinearity
        A = A_lin + eval_nonlin(x);
    end
    
    defect = f - A * x;
    
    niterations = niterations + 1;
    
    if has_nonlinearity
        if niterations >= maxiter
            warning(['iteration maximum in fixed point iteration, norm(defect): ', ...
                num2str(norm(defect))])
            break
        end
        
        if norm(defect) > 1e12
            warning('fixed point iteration diverged!')
            break
        end
        
        if verbose > 1
            disp([num2str(niterations), ': norm(defect) = ', num2str(norm(defect))]);
        elseif verbose
            fprintf('.');
        end
    end
end

if has_nonlinearity && verbose
    fprintf('\n');
    disp(['terminate fixed point iteration after ', num2str(toc(start_iteration)), ...
        ' seconds.']);
end

end