implement orthomin(1) in matlab

main.m
function [x, flag, relres, iter, resvec] = orthomin1(A, b, rtol, maxit)
    % Initialize vectors and variables
    n = length(b);
    x = zeros(n, 1);
    r = b;
    p = r;
    rsold = r' * r;
    for iter = 1:maxit
        % Orthomin(1) iteration
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        resvec(iter) = sqrt(rsnew);
        relres = sqrt(rsnew) / norm(b);
        if relres < rtol
            break;
        end
        beta = rsnew / rsold;
        rsold = rsnew;
        p = r + beta * p;
    end
    
    % Check convergence
    if relres <= rtol
        flag = 0; % Converged
    else
        flag = 1; % Not converged
    end
end
737 chars
32 lines

related categories

gistlibby LogSnag