Advanpix vs. Maple – Direct Sparse Solvers Comparison

by Pavel Holoborodko on October 3, 2013

MATLAB doesn’t support sparse matrices computing with high precision.

As a result, we have no other way but compare our performance with famous mathematical packages Maple and Mathematica, which do have support for sparse matrices with arbitrary precision elements.

Table below show timing comparisons of solving sparse linear systems of moderate size (min 1000x1000) by Advanpix toolbox and Maple 17 using quadruple precision.

Advanpix vs. Maple: Sparse LU, 34 digits (quadruple precision).
Matrix NameSizeNNZ Advanpix (sec) Maple (sec)Speed Gain (times)
nnc13741374 x 1374 8606 0.2009 20.858103.8
orsreg_12205 x 2205 14133 1.4409 22.94815.9
watt__11856 x 1856 11360 0.6724 649.308965.7
watt__11856 x 1856 11360 0.6723 342.282509.1
lns_39373937 x 3937 25407 1.3501 432.482320.3
lnsp39373937 x 3937 25407 1.2554 433.324345.2
sherman11000 x 1000 3750 0.0686 2.18431.8
sherman21080 x 1080 23094 1.6683 231.724138.9
sherman35005 x 5005 20033 2.0002 335.824167.9
sherman41104 x 1104 3786 0.0583 1.77830.5
sherman53312 x 3312 20793 0.8152 81.15299.5

Read More

{ 0 comments }

Direct Solvers for Sparse Matrices

by Pavel Holoborodko on September 5, 2013

Newest version of toolbox includes sparse linear system solver – function mldivide or operator "\".

Solver analyzes input matrix and automatically uses the most suitable decomposition:

  1. Cholesky LDLT factorization is applied if input matrix is symmetrical/hermitian and has real positive diagonal.
  2. LU factorization is applied if LDLT failed or directly if matrix is obviously non-SPD.
    (Algorithm is based on ideas from SuperLU and uses column approximate minimum degree permutation)
  3. QR factorization for rectangular matrices to handle least squares problems (beta).
    (Left-looking QR decomposition, also uses colamd pre-ordering)

All methods are capable to compute with arbitrary precision, in real and complex cases. Tables below present timing for solving sparse linear systems with quadruple precision (34 decimal digits).
See the timing tests

{ 2 comments }

Performance of Array Manipulation Operations

by Pavel Holoborodko on August 2, 2013

Besides sparse matrices support latest version of toolbox also includes significant changes in performance optimization of the matrix / array manipulation operations. This post is attempt to highlight some of the key points & reasons of the development.

As we noted before [1, 2], MATLAB has somewhat limited object oriented capabilities. The good feature that MATLAB has in this area – it allows elements of dense matrix / n-dim array to be of custom-type (e.g. multiprecision number). In this case MATLAB automatically takes care of basic manipulation operations (indexing, assigning, concatenation, deletion, etc.).

This nice feature removes huge amount of work from toolbox’s developers – so that we can concentrate on mathematical functionality rather than spending time re-inventing the wheel.

Although we have been using this functionality from the start, we decided to abandon it in the recent version. The reason is simple: MATLAB’s implementation of array manipulation operations (with custom objects as elements) suffers from enormous loss in performance!
Read More

{ 0 comments }

Basic Support for Sparse Matrices

by Pavel Holoborodko on July 31, 2013

Starting from version 3.5.0 toolbox has support for sparse matrices with multiple precision elements.

Please note, this is the first beta version with limited functionality – at the moment toolbox can only create and display multiprecision sparse matrices:

>> S = zeros(7,3);
>> S(2,1) = eps;
>> S(5,1) = sqrt(2);
>> S(3,2) = pi;
>> S(2,3) = exp(1);
>> S(5,3) = log(2);
>> S(6,3) = sqrt(-1);
>> S = sparse(S)
 
S =
 
   (2,1)          2.22044604925031e-16                         
   (5,1)               1.4142135623731                         
   (3,2)              3.14159265358979                         
   (2,3)              2.71828182845905                         
   (5,3)             0.693147180559945                         
   (6,3)                             0 +                     1i
 
>> A = mp(S)    % Convert S to multiprecision sparse matrix of 34 digits (default)
                % Note accuracy of all constants are extended automatically.
 
A =
 
   (2,1)    1.925929944387235853055977942584927e-34
   (5,1)    1.414213562373095048801688724209698
   (3,2)    3.141592653589793238462643383279503
   (2,3)    2.718281828459045235360287471352662
   (5,3)    0.693147180559945309417232121458177
   (6,3)    0 +     1i

From now on we will extend sparse matrices support to the point of full transparency & flexibility as we have done for dense matrices.

Next steps are: element-wise arithmetic and mathematical functions, basic matrix-matrix and matrix-vector operations. Major goal after that will be to add iterative Krylov subspace methods.

Although such (currently) basic support may seem to be trivial, its implementation is very far from elementary. In fact, sparse matrices support requires by order of magnitude more efforts than dense matrices. Below we outline the major reasons why it is so.
Read More

{ 0 comments }

Citations Digest v.2

by Pavel Holoborodko on July 22, 2013

Recently I have got new notifications from researchers about toolbox usage in their works. Pre-prints of the papers are available by links below:

Thank you for letting me know.

{ 0 comments }

Undocumented MEX API in MATLAB

by Pavel Holoborodko on July 19, 2013

It is not a secret that MATLAB has many undocumented (or deliberately hidden) features and commands. There are even excellent website & book devoted specifically to this topic: Undocumented Matlab

However most of the findings are related to MATLAB language itself and investigations on undocumented MEX API seems to be missing (or scarce at least).

During development of our toolbox we have found lots of hidden functions which can be helpful for creating speed-efficient extensions for MATLAB using native languages.

Here we want to explain some of them in details and provide complete list of undocumented MEX API functions.

Read More

{ 9 comments }

User Feedback Forum

by Pavel Holoborodko on July 11, 2013

After almost 2 years of email based communication with users we have decided to open publicly visible forum. Where anyone can leave feedback, ask questions, request new features or report bugs (there are no bugs in toolbox of course :)).

Hopefully this will make question resolution more efficient and less time consuming for everybody. Please follow this link (or use “Support” menu in navigation bar above): Advanpix Support Forum

Don’t forget to write something there, vote for requests or give us your unique ideas on further toolbox improvement.

{ 0 comments }

Dear readers of the Advanpix blog,

My name is Denys Dutykh, a guest blogger today, and I would like to share my experience with the Multiprecision toolbox.

My research interests lie in the fields of applied mathematics and hydrodynamics. More precisely, with my collaborators we study mathematically free surface flows or, in other words, water waves.

Read More

{ 0 comments }

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.

Read More

{ 5 comments }

World Record – 200,000,000th Bernoulli Number

by Pavel Holoborodko on July 17, 2012

Recently we have improved our previous record by a factor of 2. This time we have computed 200,000,000th Bernoulli number (the previous world record was the 110,000,000th).

We have used the same program of David Harvey and hardware. Computation took ~ 79 days.

The final result, B_{200\,000\,000} occupies 1.5GB in uncompressed, textual format. In its very shortened form:

    \[ B_{200\,000\,000} = \displaystyle{-\frac{3369394027155030175\dots38089741093921539129}{394815332706046542049668428841497001870}} \]

Contact us if you are interested in downloading all of the digits of the B_{200\,000\,000}.

{ 0 comments }