Skip to content
Permalink
Newer
Older
100644 36 lines (36 sloc) 888 Bytes
Feb 23, 2017
1
function [L,U,P] = LU_pivot(A)
2
% GaussPivot: Gauss elimination pivoting
3
% x = GaussPivot(A,b): Gauss elimination with pivoting.
4
% input:
5
% A = coefficient matrix
6
% b = right hand side vector
7
% output:
8
% x = solution vector
9
[m,n]=size(A);
10
if m~=n, error('Matrix A must be square'); end
11
nb = n;
12
L=diag(ones(n,1));
13
U=A;
14
P=diag(ones(n,1));
15
% forward elimination
16
for k = 1:n-1
17
% partial pivoting
18
[big,i]=max(abs(U(k:n,k)));
19
ipr=i+k-1;
20
if ipr~=k
21
P([k,ipr],:)=P([ipr,k],:); % if the max is not the current index ipr, pivot count
22
%L([k,ipr],:)=L([ipr,k],:);
23
%U([k,ipr],:)=U([ipr,k],:);
24
end
25
for i = k+1:n
26
fik=U(i,k)/U(k,k);
27
L(i,k)=fik;
28
U(i,k:nb)=U(i,k:nb)-fik*U(k,k:nb);
29
end
30
end
31
% back substitution
32
%x=zeros(n,1);
33
%x(n)=Aug(n,nb)/Aug(n,n);
34
%for i = n-1:-1:1
35
% x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
36
%end
You can’t perform that action at this time.