What is wrong with this algorithm

Question

Crout factorization:

n=10 A = full(gallery('tridiag',n,-1,2,-1)) i = 2:(n-1); bmid = i.^2 / ((n+1).^4) b = [1+1/(n+1)^4, bmid, 6+(n^2)/(n+1)^4]' for i = 1:n L(i,1) = A(i,1) end

for j = 1:n U(1,j) = A(1,j)/L(1,1) end

for j = 2:n for i = j:n sum = 0.0 for k = 1:(j-1) sum = sum + L(i,k) * U(k,j) end L(i,j) = A(i,j) - sum end

U(j,j) = 1;
for i = (j+1):n
    sum = 0.0
    for k = 1:(j-1)
        sum = sum + L(j,k) * U(k,i)
    end
    U(j,i) = (A(j,i) - sum)/L(j,j)
end

end x=A\b error=max(abs(b-A*x))

This is taking me more than an hour when I try to perform n = 1000. And, my instructor told me that it shouldn't take that long to do Crout factorization. Please I need help with this. Thanks.

Answer 1

It's so slow because your inner $k$ loop compute the sum one element at a time. Instead you should use array slicing. Here is modified code which runs in 14 sec for $n=1000$.

tic
n=1000;
A = full(gallery('tridiag',n,-1,2,-1));
L = zeros(n); U = zeros(n);  % CHANGE: pre-allocate L and U
i = 2:(n-1);
bmid = i.^2 / ((n+1).^4);
b = [1+1/(n+1)^4, bmid, 6+(n^2)/(n+1)^4]';
for i = 1:n
    L(i,1) = A(i,1);
end

for j = 1:n
    U(1,j) = A(1,j)/L(1,1);
end

for j = 2:n
    for i = j:n
        sum = L(i,1:(j-1)) * U(1:(j-1),j); % CHANGE: use array slicing
        %sum = 0.0;
        %for k = 1:(j-1)
        %   sum = sum + L(i,k) * U(k,j);
        %end
        L(i,j) = A(i,j) - sum;
    end
    U(j,j) = 1;
    for i = (j+1):n
        sum = L(j,1:(j-1)) * U(1:(j-1),i); % CHANGE: use array slicing
        %sum = 0.0;
        %for k = 1:(j-1)
        %   sum = sum + L(j,k) * U(k,i);
        %end
        U(j,i) = (A(j,i) - sum)/L(j,j);
    end
end
x=A\b;
error=max(abs(b-A*x))
toc