generic linear-algebra functions for Number types (e.g. BigFloat)

[stevengj]

It would be nice to have linear-algebra functions for BigFloat (possibly as an external package?) and other generic Number types.

Currently, \ on a BigFloat converts to Float64; eig and qr don't work at all.eig doesn't work.

The good news is that LINPACK-style triple-loop textbook implementations (straight out of Trefethen or Golub & Van Loan) should mostly suffice. Since BigFloat operations are CPU-limited rather than memory-limited, there is no reason to implement a LAPACK/BLAS-style blocked algorithm.

cc: @andrioni

Updated 21 August 2015 by @andreasnoack

[mlubin]

This would also help for #3.

[ViralBShah]

On a wacky note - I have often wondered if we can translate LAPACK into julia code, and make it type independent.

Realistically, just writing a few loops out of the textbooks may be the quickest way forward. Marking this up for grabs.

[JeffBezanson]

Yes that would be especially good for rational matrices.

[stevengj]

@ViralBShah, realize that the whole structure of LAPACK is dictated by the assumption that the performance of non-blocked algorithms is memory-limited, not arithmetic-limited. It makes no sense to use this structure, with all of the complications that memory-optimization implies, in situations that are arithmetic-limited (i.e. as soon as you depart from hardware-accelerated numeric types).

[ViralBShah]

Yes, LAPACK is certainly designed to optimize memory access. My comment was largely tongue-in-cheek, in that if LAPACK was translated to julia, we would have fast linear algebra for the hardware-accelerated numeric types, and the linear algebra for user-defined numeric types for free.

Undoubtedly, the right thing here would be to do the trivial implementations.

Matrix multiplication of bigfloat matrices should be easy to test, since we already have julia implementations.

[andrioni]

I pretend to start working on this soon, especially since I already have some MPFR matrix code in C, but a generic Number implementations will be extremely slow on BigFloats, thanks to immutability + the immense amount of temporary variables created.

[raichu]

Are BigFloats immutable? I haven't seen it in the docs.
If it is the case, this should be related to JuliaLang/julia#4169

[andrioni]

For all purposes and functions in the standard library, they are technically immutable, but you can mutate them by calling MPFR directly.

[raichu]

How can I operate on BigFloats using MPFR directly?

Also, this sounds horrible. How can Julia be efficient by making them immutable? I thought something is immutable only when you explicitly define them to be immutable. It would result in generation of endless garbage, ctor, dtor and copy calls (which I believe is what is responsible for JuliaLang/julia#4169).

Is this going to be fixed? Otherwise, I feel like I'm better off continue using C++ for bignum operations.

[andrioni]

Just use ccall, here's the manual section about it. Every function that takes a mpfr_t will take a Ptr{BigFloat}.

[lsorber]

+1 for this feature, I would love to be able to define my own Number types and use them to compute matrix factorizations. This really hard to do in other languages such as MATLAB.

[raichu]

As I mentioned in #19, mplapack could solve this problem.

[ViralBShah]

That page was last updated in 2011!