0
Answered

Poor MCT performance on simple arithmetic

Ahmad 2 years ago updated by Pavel Holoborodko 2 years ago 3

I have the following piece of code which takes 100 times longer than simple double precision calculation on my laptop. Is this expected?

s = 0.00001;

delta = 0.9;

e = linspace(0,2,1015);

gamma = 1e-05;

mp.Digits(50);

tic;

for i = 1:100

m = mp((delta/s)^2 - 1); b = mp(-2*(1i*e-gamma)./s);

p1 = m - b.^2/8; p2 = - b*(m + 2)/2;
h1 = (b.^2/4 + m).^2;
h2 = b.^6/32+(3*b.^4*m)/8+(3*b.^2*m^2)/2+27*b.^2*m+27*b.^2+2*m^3;
S = (1/2*(h2+sqrt(h2.^2-4*h1.^3))).^(1/3);
S = 1/2*sqrt(-2/3*p1+1/3*(S+h1./S));
sol = -b./4+S+1/2*sqrt(-4*S.^2-2*p1-p2./S);

sol = double(sol);

end

toc;

P.S. this is the application of Cardano formula to find the roots of a quartic polynomial when coefficients vary.

I have bolded the parts that need to be changed if run using MATLAB double precision arithmetic.

Answer

Answer
Answered

By the way, you need to compute constants (1/3, 2/3) using extended precision: mp('1/3'), mp('2/3').

Otherwise MATLAB evaluates them using double precision. Of course, this must be pre-calculated once before the loop (same as m).

There is still a room for improvement in your code. 



+1
Under review

Naturally, extended precision computations are slower compared to double (double run in hardware, extended precision is emulated in software). 

However, how much "slower" depends entirely on the particular code and how it is implemented.

Another important slow-down comes from the MATLAB overhead. Each operation in the M-code results in data copy and passing it to toolbox. This adds non-negligible overhead.

Thus, for higher speed, you might consider re-writing your code to avoid redundant calculations (e.g. b.^2 computed many times) and to use vectorization (when one operation works with vectors & matrices instead of crunching through scalars in the loop). This will reduce MATLAB overhead and allow toolbox to handle data in parallel. You will see nice speed-up.

Thanks Pavel. I think I have used vectorization in this particular example to the extent it is possible. However you are completely right about the two other factors, namely: natural slow-down and MATLAB overhead.

Answer
Answered

By the way, you need to compute constants (1/3, 2/3) using extended precision: mp('1/3'), mp('2/3').

Otherwise MATLAB evaluates them using double precision. Of course, this must be pre-calculated once before the loop (same as m).

There is still a room for improvement in your code.