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
|
function x = GaussPivot (A, b)
[m n] = size (A);
m ~= n && error ('Matrix A must be square');
nb = n + 1;
Aug = [A b];
% forward elimination
for k = 1 : n - 1
% partial pivoting
[big i] = max (abs (Aug(k:n, k)));
ipr = i + k - 1;
if ipr ~= k
Aug([k ipr], :) = Aug([ipr k], :);
end
for i = k + 1 : n
coeff = Aug(i, k) / Aug(k, k);
Aug(i, k:nb) = Aug(i, k:nb) - coeff * Aug(k, k:nb);
end
end
% back substitution
x = zeros (n, 1);
x(n) = Aug(n, nb) / Aug(n, n);
for i = n - 1 : -1 : 1
x(i) = (Aug(i, nb) - Aug(i, i+1:n) * x(i+1:n)) / Aug(i, i);
end
end
|
