[sage-devel] multiprecision linear algebra

138 views
Skip to first unread message

Jason Grout

unread,
May 20, 2010, 7:48:03 PM5/20/10
to sage-...@googlegroups.com
A while ago someone talked about incorporating ALGLIB into Sage to
provide, among other things, multiprecision linear algebra. I just saw
another package that does this:

http://mplapack.sourceforge.net/

Does anyone have experience with mplapack? How does it compare to
ALGLIB in linear algebra (for algorithms that are implemented in both)?

Thanks,

Jason

--
To post to this group, send an email to sage-...@googlegroups.com
To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org

Sergey Bochkanov

unread,
May 21, 2010, 2:36:03 AM5/21/10
to Jason Grout
Hello, Jason.

You wrote 21 мая 2010 г., 3:48:03:

> A while ago someone talked about incorporating ALGLIB into Sage to
> provide, among other things, multiprecision linear algebra.

It was me :) I am working on this issue. I expect that interface
between double precision ALGLIB and Python will be ready in 2-3 weeks,
with multiple precision interface following several weeks later.

> Does anyone have experience with mplapack? How does it compare to
> ALGLIB in linear algebra (for algorithms that are implemented in both)?

I haven't benchmarked mplapack, but I can say something about it and
can compare it with ALGLIB:


== ALGORITHMS ========================================================

MPLAPACK includes all LAPACK algorithms, especially wide set of
generalized EVD/SVD solver. Nothing beyond LAPACK is included.

ALGLIB includes subset of LAPACK: real/complex/SPD/HPD triangular
factorizations, inverses, solvers; real/complex/Hermitian EVD, real
SVD. However, it also includes algorithms beyond linear algebra:
numerical integration, interpolation/fitting, optimization, many other
algorithms (including those which won't benefit from multiple
precision).


== MULTIPLE PRECISION SUPPORT ========================================

MPLAPACK supports several multiple precision types:
* mpf_t (GMP)
* mpfr_t (MPFR)
* double-double and quad-double (DD/QD package)
* double (standard 64-bit floating point type)

However, only mpfr_t (and double) support is important. Other multiple
precision types are almost useless due to IEEE-nonconformance:
* mpf_t is not recommended for numerical work (it leads to
non-reproducible results and have very odd rounding rules; even GMP
manual recommends to use MPFR)
* DD/QD may look attractive (very high performance for 106 and 212
bits when compared with mpfr_t), but it is no good for numerical work
too (except for BLAS internals; see page 31 of
http://www.netlib.org/lapack/lawnspdf/lawn149.pdf ).

ALGLIB supports mpfr_t and double.


== PROGRAMMING LANGUAGE ==============================================

MPLAPACK is written in C++ and uses C++ to handle multiple precision
types like built-in floating point types (i.e. it uses "a=b+c/d"
instead of directly calling MPFR functions). It leads to unnecessary
allocation of temporaries and some performance penalty (about 25-30%
for 128-bit precision).

Multiple precision ALGLIB 2.5.0 is written in C++ and suffers from the
same drawbacks. But what I develop for SAGE is a faster pure C library
which caches all temporaries (it is really easy to develop - I just
have to modify code generator which is used to generate ALGLIB
source).


== OPTIMIZATIONS =====================================================

Both packages rely on highly optimized MPFR for multiple precision
computations.

As for double precision, ALGLIB uses better optimized linear algebra
code.


== MULTITHREADING ====================================================

Both project don't support multithreading yet. Both want to :)


== LICENSE ===========================================================

MPLAPACK is LGPL 3+, ALGLIB is GPL 2+.


--
With best regards,
Sergey mailto:sergey.b...@alglib.net

Sergey Bochkanov

unread,
May 21, 2010, 5:17:01 AM5/21/10
to Sergey Bochkanov
correction to my previous message:

> It leads to unnecessary allocation of temporaries and some
> performance penalty (about 25-30% for 128-bit precision).

Sorry, I made a mistake when estimating MPLAKACK performance. "25-30%"
are ALGLIB's, not MPLAPACK's.


I've never measured MPLAPACK performance directly, but looking at its
RGEMM code and measuring allocation penalty for code like this:

> mpfr_t tmp1, tmp2;
> mpfr_init2(tmp1, Precision);
> mpfr_init2(tmp2, Precision);
> mpfr_mul(tmp1, mp_b, mp_c, GMP_RNDN);
> mpfr_add(tmp2, mp_a, mp_tmp, GMP_RNDN);
> mpfr_set(mp_d, tmp1, GMP_RNDN);
> mpfr_clear(tmp1);
> mpfr_clear(tmp2);

(which reproduces allocations done by MPLAPACK during "a=a+b*c" update
in the RGEMM call) and comparing it with

> mpfr_mul(mp_tmp, mp_b, mp_c, GMP_RNDN);
> mpfr_add(mp_d, mp_a, mp_tmp, GMP_RNDN);

which does the same thing, but more efficiently, I can conclude that
MPLAPACK allocation penalty will be 100% for 128-bit precision code
(i.e. it will run twice longer), 50% for 256-bit precision.


As for ALGLIB, it really has 25-30% allocation penalty for 128-bit
code. The reason behind such difference from MPLAPACK is that ALGLIB
tries to allocate new multiple precision variables from internal
cache, but still relies on OOP, which slows it down.

So when I said that "ALGLIB ... suffers from the same OOP drawbacks",
I was a bit incorrect - it suffers from OOP drawbacks, but from
_other_ drawbacks.

Fredrik Johansson

unread,
May 21, 2010, 5:39:38 AM5/21/10
to sage-...@googlegroups.com

MPFR has a fused multiply-add function, so d = a+b*c can be written mpfr_fma(d,a,b,c, MPFR_RNDN). This is almost certainly what you'd want to base MPFR linear algebra around.

Fredrik

Bill Hart

unread,
May 21, 2010, 9:26:38 AM5/21/10
to sage-devel
Oh my, I never noticed this! I'd been looking for mpfr_addmul in line
with GMP!! I see it is slightly different in that it doesn't do a = a
+ b*c.

I can actually speed some code up with this!!

Bill.

On 21 May, 10:39, Fredrik Johansson <fredrik.johans...@gmail.com>
wrote:
> >  Sergey                          mailto:sergey.bochka...@alglib.net

Bill Hart

unread,
May 21, 2010, 9:29:37 AM5/21/10
to sage-devel


On 21 May, 07:36, Sergey Bochkanov <sergey.bochka...@alglib.net>
wrote:
+1

> *  DD/QD  may  look  attractive (very high performance for 106 and 212
> bits  when compared with mpfr_t), but it is no good for numerical work
> too    (except    for    BLAS    internals;    see    page    31    ofhttp://www.netlib.org/lapack/lawnspdf/lawn149.pdf ).
>

+1 - very important point

> ALGLIB supports mpfr_t and double.
>
> == PROGRAMMING LANGUAGE ==============================================
>
> MPLAPACK  is  written in C++ and uses C++ to handle multiple precision
> types  like  built-in  floating  point  types  (i.e. it uses "a=b+c/d"
> instead  of  directly calling MPFR functions). It leads to unnecessary
> allocation  of  temporaries and some performance penalty (about 25-30%
> for 128-bit precision).
>
> Multiple precision ALGLIB 2.5.0 is written in C++ and suffers from the
> same drawbacks. But what I develop for SAGE is a faster pure C library
> which  caches  all  temporaries (it is really easy to develop - I just
> have  to  modify  code  generator  which  is  used  to generate ALGLIB
> source).
>
> == OPTIMIZATIONS =====================================================
>
> Both  packages  rely  on  highly optimized MPFR for multiple precision
> computations.
>
> As  for  double precision, ALGLIB uses better optimized linear algebra
> code.
>
> == MULTITHREADING ====================================================
>
> Both project don't support multithreading yet. Both want to :)
>
> == LICENSE ===========================================================
>
> MPLAPACK is LGPL 3+, ALGLIB is GPL 2+.
>
> --
> With best regards,
>  Sergey                          mailto:sergey.bochka...@alglib.net
>
> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group athttp://groups.google.com/group/sage-devel
> URL:http://www.sagemath.org

Fredrik Johansson

unread,
May 21, 2010, 9:46:37 AM5/21/10
to sage-...@googlegroups.com
On Fri, May 21, 2010 at 3:26 PM, Bill Hart <goodwi...@googlemail.com> wrote:
Oh my, I never noticed this! I'd been looking for mpfr_addmul in line
with GMP!! I see it is slightly different in that it doesn't do a = a
+ b*c.

I can actually speed some code up with this!!

Bill.

Nice. Indeed, I specified the function incorrectly; it computes a*b+c, not a+b*c (a+b*c wouldbe more logical for a function named "addmul").

Fredrik

--
To post to this group, send an email to sage-...@googlegroups.com
To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com

Bill Hart

unread,
May 21, 2010, 9:56:19 AM5/21/10
to sage-devel
I just took a look at their website:

http://mplapack.sourceforge.net/

* I'm personally confused about the licensing (may be a reflection on
me, not mplapack). The MBLAS part is BSD, but the MPLAPACK part is
LGPL v3+, but the authors say they are going to replace this later.
But with what? They list "Including in a package system like FreeBSD
ports" in their future goals, though this doesn't necessarily imply
they are moving towards BSD licensing.

* Some nonsense about GMP (and incidentally no mention of MPIR - note
how many other bignum libraries *are* mentioned).

* From a code quality point of view, their code looks pretty good to
me, commented well, in english. I didn't do an extensive survey
though.

* They implement a lot of functionality, however I wouldn't personally
use this as a metric. Sergey Bochkanov's automatic code generation is
exceptionally clever and should allow him to be streets ahead of
mplapack very quickly

* I would benchmark carefully

I became acquainted with ALGLIB and Sergey through his interaction
with the MPIR project. Even if at some point mplapack is better than
ALGLIB, I reckon Sergey would be exceptionally helpful in addressing
that. He's a smart guy and very supportive of what Sage is doing!

Bill.

On 21 May, 00:48, Jason Grout <jason-s...@creativetrax.com> wrote:
> A while ago someone talked about incorporating ALGLIB into Sage to
> provide, among other things, multiprecision linear algebra.  I just saw
> another package that does this:
>
> http://mplapack.sourceforge.net/
>
> Does anyone have experience with mplapack?  How does it compare to
> ALGLIB in linear algebra (for algorithms that are implemented in both)?
>
> Thanks,
>
> Jason
>
> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group athttp://groups.google.com/group/sage-devel
> URL:http://www.sagemath.org

Bill Hart

unread,
May 21, 2010, 9:57:20 AM5/21/10
to sage-devel
Ah, I missed that. Thanks for the clarification!

On 21 May, 14:46, Fredrik Johansson <fredrik.johans...@gmail.com>
wrote:
> On Fri, May 21, 2010 at 3:26 PM, Bill Hart <goodwillh...@googlemail.com>wrote:
>
> > Oh my, I never noticed this! I'd been looking for mpfr_addmul in line
> > with GMP!! I see it is slightly different in that it doesn't do a = a
> > + b*c.
>
> > I can actually speed some code up with this!!
>
> > Bill.
>
> > Nice. Indeed, I specified the function incorrectly; it computes a*b+c, not
>
> a+b*c (a+b*c wouldbe more logical for a function named "addmul").
>
> Fredrik
>
> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group athttp://groups.google.com/group/sage-devel
> URL:http://www.sagemath.org

Jason Grout

unread,
May 21, 2010, 10:41:18 AM5/21/10
to sage-...@googlegroups.com
On 05/21/2010 08:56 AM, Bill Hart wrote:

> I became acquainted with ALGLIB and Sergey through his interaction
> with the MPIR project. Even if at some point mplapack is better than
> ALGLIB, I reckon Sergey would be exceptionally helpful in addressing
> that. He's a smart guy and very supportive of what Sage is doing!


Thanks for everyone's answers!

This last phrase is also an important comparison point:

ALGLIB: Sergey has already been writing code and is willing to do work
to integrate this with Sage

MPLAPACK: ???


Thanks! I'm really excited about this.

Jason

--
To post to this group, send an email to sage-...@googlegroups.com
To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com

Bill Hart

unread,
May 21, 2010, 11:09:35 AM5/21/10
to sage-devel
For comparison, I just had a look at the MPLAPACK author's previous
project:

http://sdpa.indsys.chuo-u.ac.jp/sdpa/

I *barely* know anything about this stuff - once did a basic course on
this sort of thing as an undergrad.

However, SDPA might be an indication of what to expect from the same
author.

Note there are parallel versions of this package and the whole thing
looks very mature (should we be asking the engineers whether they'd
like this in Sage?) The whole thing appears to achieve uberdom in the
Open Source world by being one of the few projects that looks
relatively "complete" instead of "in development". In short, the
author seems very competent indeed.

I still give a tentative +1 to incorporating ALGLIB rather than
MPLAPACK in Sage for the reasons given.

Bill.
> For more options, visit this group athttp://groups.google.com/group/sage-devel
> URL:http://www.sagemath.org

Maho NAKATA

unread,
May 22, 2010, 2:52:23 AM5/22/10
to sage-devel
Hi Jason

I'm Nakta, Maho, a principal developer of MPACK.

I've never used ALGLIB before...

thanks,

On 5月21日, 午前8:48, Jason Grout <jason-s...@creativetrax.com> wrote:
> A while ago someone talked about incorporating ALGLIB into Sage to
> provide, among other things, multiprecision linear algebra.  I just saw
> another package that does this:
>
> http://mplapack.sourceforge.net/
>
> Does anyone have experience with mplapack?  How does it compare to
> ALGLIB in linear algebra (for algorithms that are implemented in both)?
>
> Thanks,
>
> Jason
>
> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group athttp://groups.google.com/group/sage-devel
> URL:http://www.sagemath.org

Maho NAKATA

unread,
May 22, 2010, 2:58:54 AM5/22/10
to sage-devel
Hi Bill
thanks for your comments

On 5月21日, 午後10:56, Bill Hart <goodwillh...@googlemail.com> wrote:
> I just took a look at their website:
>
> http://mplapack.sourceforge.net/
>
> * I'm personally confused about the licensing (may be a reflection on
> me, not mplapack). The MBLAS part is BSD, but the MPLAPACK part is
> LGPL v3+, but the authors say they are going to replace this later.

sorry, what i meant is I'll replace all the parts by mBSD style
license.

> But with what? They list "Including in a package system like FreeBSD
> ports" in their future goals, though this doesn't necessarily imply
> they are moving towards BSD licensing.

I'm also a FreeBSD ports committer, so I'd like to add it for FreeBSD
ports tree
for easier compilation.

> * Some nonsense about GMP (and incidentally no mention of MPIR - note
> how many other bignum libraries *are* mentioned).
>
> * From a code quality point of view, their code looks pretty good to
> me, commented well, in english. I didn't do an extensive survey
> though.
>
> * They implement a lot of functionality, however I wouldn't personally
> use this as a metric. Sergey Bochkanov's automatic code generation is
> exceptionally clever and should allow him to be streets ahead of
> mplapack very quickly
>
> * I would benchmark carefully

You should, recently I wrote a paper in Japanese, and my result
is something like

Raxpy
DD: 136MFlops (250MFlops)
QD: 14.2MFlops (41.2MFlops)
GMP 14.2MFlops (38.8MFlops)
GotoBlas 1.5GFlops

Rgemm
DD: 141M
QD 13.9M
GMP 14.4M
GotoBLAS 42.5G :-O

Maho NAKATA

unread,
May 22, 2010, 2:59:48 AM5/22/10
to sage-devel
Hi Bill,

Of course, I'd like to provide MScaLAPACK etc.

Thanks,

Maho NAKATA

unread,
May 22, 2010, 3:02:18 AM5/22/10
to sage-devel
Hi Bill,

Currently usable MLAPACK routines are only 90 among total about 660.
Just fortran to C conversion is completed.

Purpose of MPACK
I'd like to provide reference build. based on that others can provide
faster implementations like GotoBLAS, ATLAS, intel MKL etc.

For Benchmark issue...
There is an implementation of Rgemm DD by Mukunoki-san and using
GPGPU. It's 30 times faster than my reference routine.

Thanks

On 5月21日, 午後10:56, Bill Hart <goodwillh...@googlemail.com> wrote:

Maho NAKATA

unread,
May 22, 2010, 3:04:26 AM5/22/10
to sage-devel
Hi Bill

DD 136Mflops (250MFlops)
means that 250MFlops using OpenMP parallelism.

Note that above results are using 0.6.4. It would be bit different
using 0.6.5.
Also seems to depend which compiler you use.

thanks

Maho NAKATA

unread,
May 22, 2010, 3:29:44 AM5/22/10
to sage-...@googlegroups.com
Hi Bill again (sorry)

With 0.6.5, you can easily benchmark by just building and installing.
you'll find benchmarks in ${PREFIX}/share/mpack/bench/
in my case,
Raxpy.dd Raxpy.dd_ref ...
Thanks

Subject: [sage-devel] Re: multiprecision linear algebra
Date: Sat, 22 May 2010 00:04:26 -0700 (PDT)

Maho NAKATA

unread,
May 31, 2010, 3:23:26 AM5/31/10
to sage-...@googlegroups.com
Reply all
Reply to author
Forward
0 new messages