1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
% cgsolve.m
%
% Solve a symmetric positive definite system Ax = b via conjugate gradients.
%
% Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
%
% A - Either an NxN matrix, or a function handle.
%
% b - N vector
%
% tol - Desired precision. Algorithm terminates when
% norm(Ax-b)/norm(b) < tol .
%
% maxiter - Maximum number of iterations.
%
% verbose - If 0, do not print out progress messages.
% If and integer greater than 0, print out progress every 'verbose' iters.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%
function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
if (nargin < 5), verbose = 1; end
implicit = isa(A,'function_handle');
x = zeros(length(b),1);
r = b;
d = r;
delta = r'*r;
delta0 = b'*b;
numiter = 0;
bestx = x;
bestres = sqrt(delta/delta0);
while ((numiter < maxiter) && (delta > tol^2*delta0))
% q = A*d
if (implicit), q = A(d); else q = A*d; end
alpha = delta/(d'*q);
x = x + alpha*d;
if (mod(numiter+1,50) == 0)
% r = b - Aux*x
if (implicit), r = b - A(x); else r = b - A*x; end
else
r = r - alpha*q;
end
deltaold = delta;
delta = r'*r;
beta = delta/deltaold;
d = r + beta*d;
numiter = numiter + 1;
if (sqrt(delta/delta0) < bestres)
bestx = x;
bestres = sqrt(delta/delta0);
end
if ((verbose) && (mod(numiter,verbose)==0))
disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ...
numiter, bestres, sqrt(delta/delta0)));
end
end
if (verbose)
disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres));
end
x = bestx;
res = bestres;
iter = numiter;
|
