文档介绍：Introduction to puting -- A Matrix Vector Approach Using Matlab Written by Charles Loan
陈文斌
复旦大学
Chapter 6 Linear systems
Triangular Problems
Banded Problems
Full Problems
Analysis
Triangular Problems
Forward substitution
for i=1:n
x(i)=b(i);
for j=1:i-1
x(i)=x(i)-L(i,j)*x(j);
end
x(i)=x(i)/L(i,i);
end
x(1)=b(1)/L(1,1);
for i=2:n
x(i)=(b(i)-L(i,1:i-1)*x(1:i-1))/L(i,i);
end
Row oriented
floops
Column-oriented
saxpy
function x = LTriSol(L,b)
n = length(b);
x = zeros(n,1);
for j=1:n-1
x(j) = b(j)/L(j,j);
b(j+1:n) = b(j+1:n) - L(j+1:n,j)*x(j);
end
x(n) = b(n)/L(n,n);
The version involves n2 flops, just like the row-oriented, dot product version developed eariler.
Backward Substitution
x(n)=b(n)/U(n,n);
for i=n-1:-1:1
x(i)=(b(i)-U(i,i+1:n)*x(i+1:n))/U(i,i)
end
function x = UTriSol(U,b)
n = length(b);
x = zeros(n,1);
for j=n:-1:2
x(j) = b(j)/U(j,j);
b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = b(1)/U(1,1);
Column-oriented saxpy
Backward Substitution
Multiple Right-Hand Sides
X=zeros(n,r);
for k=1:r
X(:,k)=LTriSol(L,B(:,k));
end
x = zeros(n,1);
For k=1:r
for j=1:n-1
X(j,k) = B(j,k)/L(j,j);
B(j+1:n,k) = B(j+1:n,k) - L(j+1:n,j)*X(j,k);
end
X(n,k) =B(n,k)/L(n,n);
End