www.pudn.com > starter.zip > conjGrad.m, change:2011-01-04,size:1845b


function [x,k,res,negCurv] = cg(A,b,optTol,maxIter,verbose,precFunc,precArgs,matrixVectFunc,matrixVectArgs) 
% [x,k,res,negCurv] = 
% cg(A,b,optTol,maxIter,verbose,precFunc,precArgs,matrixVectFunc,matrixVect 
% Args) 
% Linear Conjugate Gradient, where optionally we use 
% - preconditioner on vector v with precFunc(v,precArgs{:}) 
% - matrix multipled by vector with matrixVectFunc(v,matrixVectArgs{:}) 
 
x = zeros(size(b)); 
r = -b; 
 
% Apply preconditioner (if supplied) 
if nargin >= 7 && ~isempty(precFunc) 
    y = precFunc(r,precArgs{:}); 
else 
    y = r; 
end 
 
ry = r'*y; 
p = -y; 
k = 0; 
 
res = norm(r); 
done = 0; 
negCurv = []; 
while res > optTol & k < maxIter & ~done 
    % Compute Matrix-vector product 
    if nargin >= 9 
        Ap = matrixVectFunc(p,matrixVectArgs{:}); 
    else 
        Ap = A*p; 
    end 
    pAp = p'*Ap; 
 
    % Check for negative Curvature 
    if pAp <= 1e-16 
        if verbose 
            fprintf('Negative Curvature Detected!\n'); 
        end 
         
        if nargout == 4 
           if pAp < 0 
              negCurv = p; 
              return 
           end 
        end 
         
        if k == 0 
            if verbose 
                fprintf('First-Iter, Proceeding...\n'); 
            end 
            done = 1; 
        else 
            if verbose 
                fprintf('Stopping\n'); 
            end 
            break; 
        end 
    end 
 
    % Conjugate Gradient 
    alpha = ry/(pAp); 
    x = x + alpha*p; 
    r = r + alpha*Ap; 
     
    % If supplied, apply preconditioner 
    if nargin >= 7 && ~isempty(precFunc) 
        y = precFunc(r,precArgs{:}); 
    else 
        y = r; 
    end 
     
    ry_new = r'*y; 
    beta = ry_new/ry; 
    p = -y + beta*p; 
    k = k + 1; 
 
    % Update variables 
    ry = ry_new; 
    res = norm(r); 
end 
end