Matlab's ltitr vs. Julia

ltitr_test.jl

# Test for ltitr implemented in Julia. I get about 45 ms.

function ltitr(A, B, U)
    N = size(U, 2)
    X = zeros(size(A, 1), N)
    for i = 2:N
        X[:, i] = A*X[:, i-1] + B*U[:, i-1]
    end
    return X
end

# Generated from Matlab's drss.
A = [
    0.0948    0.3140   -0.0264    0.3122    0.0508   -0.1067
    0.1942   -0.0640    0.3212   -0.3044   -0.1294    0.1517
    0.1020    0.0272   -0.2001    0.0469   -0.0988   -0.3590
    0.2762   -0.3470    0.1809    0.1321   -0.0954    0.2293
   -0.0388   -0.1235    0.1149   -0.0617    0.2265   -0.1043
   -0.2912    0.1957    0.0466    0.3083   -0.0804   -0.4048
    ]

B = [
   -0.9738    0.5583    0.7618
   -1.6347    0.7752    1.2038
         0   -0.7823    1.5441
    3.2967   -0.5119   -1.0941
   -0.4837    0.3074    1.4009
    0.3281   -0.7417   -0.7114
    ]

N = 10000
U = randn(3, N)

gc_disable()
tic()
for i=1:50 X = ltitr(A, B, U); end
timing = toc()
println("ltitr took ", (timing / 50) * 1000, " ms")

ltitr_test.m

function ltitr_test
% Function to test built-in and hand-coded versions of ltitr. I get about
% 3.5 ms for the built-in and 48 ms for the hand-coded.

A = [
    0.0948    0.3140   -0.0264    0.3122    0.0508   -0.1067
    0.1942   -0.0640    0.3212   -0.3044   -0.1294    0.1517
    0.1020    0.0272   -0.2001    0.0469   -0.0988   -0.3590
    0.2762   -0.3470    0.1809    0.1321   -0.0954    0.2293
   -0.0388   -0.1235    0.1149   -0.0617    0.2265   -0.1043
   -0.2912    0.1957    0.0466    0.3083   -0.0804   -0.4048
    ];

B = [
   -0.9738    0.5583    0.7618
   -1.6347    0.7752    1.2038
         0   -0.7823    1.5441
    3.2967   -0.5119   -1.0941
   -0.4837    0.3074    1.4009
    0.3281   -0.7417   -0.7114
     ];

N = 10000;
U = randn(3, N)';

tic
for i=1:50
    X = ltitr(A, B, U); 
end
timing = toc*1000/50;
disp(['Built-in took ', num2str(timing), ' ms.']);

for i = 1:50
    X = my_ltitr(A, B, U');
end
timing = toc*1000/50;
disp(['Non-built-in took ', num2str(timing), ' ms.']);

end

function X = my_ltitr(A, B, U)
    N = size(U, 2);
    X = zeros(size(A, 1), N);
    for i = 2:N
        X(:, i) = A*X(:, i-1) + B*U(:, i-1);
    end
end