Fast Quadruple Precision Computations in MATLAB

by Pavel Holoborodko on January 20, 2013

In order to improve performance we have made important changes to the toolbox recently.

Now toolbox is heavily optimized for computations with quadruple precision.

In more details, new core engine of toolbox (written in Assembler/C++) has two different variants for every routine:

  • The first variant is truly precision independent, capable of computing results with any required precision (same as before).
  • The other is highly optimized for quadruple precision (128-bit or 34 decimal digits) computations.

Toolbox automatically switches to the fast quadruple routines if the user asks toolbox to operate with the precision of 34 decimal digits. No additional actions are required.

Quadruple precision routines have an advantage of knowing size of numbers in advance thus avoiding expensive memory allocations. Also this allows more efficient usage of hardware capabilities of modern CPUs, even exploiting some of the advanced instructions, not to mention possibilities of harnessing power of the GPU!

Quadruple precision numbers and operations are fully compliant with IEEE 754-2008 standard.

Comparison with Symbolic Math Toolbox

We compare optimized quadruple precision mode in our toolbox with Symbolic Math Toolbox from MathWorks. We use 64-bit MATLAB 2012b, Windows 7 and PC based on Core 2 Quad CPU.

At first we do initialization for both toolboxes to operate with quadruple precision:

% Symbolic Math Toolbox uses 9 guard digits. 
% Hence we should use 25 digits of precision 
% to get 'true' quadruple precision (34 = 25+9)
>> digits(25);           
 
% Setup 34-digits mode in Multiprecision Computing Toolbox.
>> mp.Digits(34);
>> mp.GuardDigits(0); 
 
% Show all digits of the numbers
>> format longG
 
% Prepare test matrices
>> A = rand(100);       % 100 x 100 full rank test matrix
>> X = mp(A);           
>> Y = vpa(A);

Then we do the tests for svd, eig and qr functions (much more detailed comparisons can be found here: 500×500 and 100×100):

Singular Value Decomposition

>> tic; [U,S,V] = svd(Y); toc;       % Symbolic Math Toolbox
Elapsed time is 59.357723 seconds. 
>> norm(Y-U*S*V',1)
3.061145275981987186191661e-31
 
>> tic; [U,S,V] = svd(X); toc;       % Multiprecision Computing Toolbox
Elapsed time is 0.348560 seconds.    % x170 times faster 
>> norm(X-U*S*V',1)                 
5.718086004885703247723198511534649e-31

Eigenvalues and Eigenvectors

>> tic; [V,D] = eig(Y); toc;         % Symbolic Math Toolbox
Elapsed time is 108.474148 seconds. 
>> norm(Y*V - V*D,1)
1.664003471950571777040365e-30
 
>> tic; [V,D] = eig(X); toc;         % Multiprecision Computing Toolbox
Elapsed time is 1.114424 seconds.    % x100 times faster 
>> norm(X*V - V*D,1)
3.520599938339867139386327679045247e-31

QR Decomposition

>> tic; [Q,R] = qr(Y); toc;      % Symbolic Math Toolbox
Elapsed time is 7.269699 seconds. 
>> norm(Y-Q*R,1)
3.608711233295583179663639e-31
 
>> tic; [Q,R] = qr(X); toc;          % Multiprecision Computing Toolbox
Elapsed time is 0.100289 seconds.    % x72 times faster 
>> norm(X-Q*R,1)
4.758852521959335590660505497505944e-32

Large Matrices Computing

Fast ‘quadruple’ precision mode even makes possible processing sufficiently large matrices in reasonable time:

>> A = mp(rand(10000,400));               % 10000 x 400 test matrix
 
>> tic; [U,S,V] = svd(A,'econ'); toc;
Elapsed time is 113.552919 seconds.       % ~2 minutes
 
>> norm(A - U*S*V',1)                     % Check accuracy of the result
8.332401597108828366965178092631891e-29

By our estimation, Symbolic Math Toolbox requires 7-8 hours to compute the same result (in fact, we stopped it after 5 hours of computations).

Check our User Manual for more detailed comparison with Symbolic Math Toolbox.

{ 5 comments… read them below or add one }

Maarten J. van der Burgt March 12, 2019 at 7:58 pm

I would like to calculate an integral in quadruple precision. I have an old Mathcad 7 program on my old computer but cannot transfer it to my new computer a HP Pavillion with an Intel core 15 processor. What to I have to do in order to do my calculation in a free trial? Please comment.

Best regards

Maarten J. van der Burgt

Reply

Pavel Holoborodko March 12, 2019 at 8:25 pm

Our toolbox supports MATLAB only. Unfortunately you will not be able to run Mathcad program in MATLAB. You need to rewrite your program in MATLAB. Then toolbox can be used to run program in extended precision.

Reply

maarten van der Burgt March 13, 2019 at 1:52 am

Thanks a lot for the prompt reply,

Maarten

Reply

Alan September 5, 2021 at 3:01 am

Is there a routine that takes advantage of a matrix being hermitian? I am looking for a quad precision eigensolver for complex hermitian matrices and wondering if this will work.

Reply

Pavel Holoborodko September 5, 2021 at 10:47 am

The toolbox’s eigensolver analyses input matrix and uses the most suitable solver for it.
Yes, of course, it uses special solver for Hermitian matrices. Just call “eig” function, it will check matrix properties and will use Hermitian eigensolver automatically.

This post outlines architecture of our eigensolver (and special cases it handles):
https://www.advanpix.com/2016/10/20/architecture-of-eigenproblem-solver/

Reply

Leave a Comment

Previous post:

Next post: