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

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
```