Introduction to Computational Linear Algebra

72 Introduction to Computational Linear Algebra

Algorithm 3.3 Optimal Storage Algorithm for Scaled Partial Pivoting

function [L,U,P]=mylu(A,tol)

% Input: A a square matrix (invertible or not)

% Output: either:

% L unit lower triangular

% U upper triangular

% P: Permutation matrix

% P*A=L*U and in case

% A singular display matrix rank

n=length(A);r=n;

P=eye(n);

% Compute the scales

s=max(abs(A’));s=s(:);

IV=1:n;

for k=1:n

[p,ip]=max(abs(A(IV(k:n),k))./s(IV(k:n)));

ip=ip+k-1;

if abs(A(IV(ip),k))>tol

IV([ip k])=IV([k ip]);

piv=sign(A(IV(k),k))*p*s(IV(k));

piv=1/piv;% 1 flop

A(IV(k+1:n),k)=piv*A(IV(k+1:n),k);%(n-k) flops

B=A(IV(k+1:n),k)*A(IV(k),k+1:n);%(n-k)^2 flops ...