6

Where can I get a fast linear system solver written in D? It should be able to take a square matrix A and a vector b and solve the equation Ax = b for b and, ideally, also perform explicit inversion on A. I have one I wrote myself, but it's pretty slow, probably because it's completely cache-naive. However, for my use case, I need something with the following absolute, non-negotiable requirements, i.e. if it doesn't meet these, then I don't otherwise care how good it otherwise is:

  1. Must be licensed public domain, Boost license, or some similar permissive license. Ideally it should not require attribution in binaries (i.e. not BSD), though this point is somewhat negotiable.

  2. Must be written in pure D or easily translatable to pure D. Inscrutable Fortran code (i.e. LAPACK) is not a good answer no matter how fast it is.

  3. Must be optimized for large (i.e. n > 1000) systems. I don't want something that's designed for game programmers to solve 4x4 matrices really, really fast.

  4. Must not be inextricably linked to a huge library of stuff I don't need.

Edit: The reason for these seemingly insane requirements is that I need this code for a permissively licensed open source library that I don't want to have any third-party dependencies.

dsimcha
  • 67,514
  • 53
  • 213
  • 334
  • 1
    What properties does A have? Symmetric? Positive? Positive definite? Banded? Sparse? You ask for something fast, public domain and exclude LAPACK! It doesn't sound like you actually want a solution! Are you aware that all performant linear solvers are derived from The Handbook? – David Heffernan Feb 08 '11 at 20:01
  • @Baxissimo "The Handbook for Automatic Computation: Linear Algebra" originally published in 1972, authored by J. H. Wilkinson and C. Reinsch – David Heffernan Feb 09 '11 at 10:44

1 Answers1

3

If you don't like Fortran code, one reasonably fast C++ dense matrix library with modest multi-core support, well-written code and a good user-interface is Eigen. It should be straightforward to translate its code to D (or to take some algorithms from it).

And now my "think about your requirements": there is a reason why "everyone" (Mathematica, Matlab, Maple, SciPy, GSL, R, ...) uses ATLAS / LAPACK, UMFPACK, PARDISO, CHOLMOD etc. It is hard work to write fast, multi-threaded, memory-efficient, portable and numerically stable matrix solvers (trust me, I have tried). A lot of this hard work has gone into ATLAS and the rest.

So my approach would be to write bindings for the relevant library depending on your matrix type, and link from D against the C interfaces. Maybe the bindings in multiarray are enough (I haven't tried). Otherwise, I'd suggest looking at another C++ library, namely uBlas and the respective bindings for ideas.

stephan
  • 10,104
  • 1
  • 51
  • 64
  • The problem with LAPACK, etc. is I need this for one small section of code and don't want a dependency. Also, I can't use Eigen because it's LGPL, so my D translation would be covered by the copyleft. – dsimcha Feb 09 '11 at 13:39
  • @dsimcha: Makes sense. Maybe looking at Eigen's partial pivot LU and triangularization source might be helpful. I doubt that the algorithms are covered by the LGPL (given that they are most likely based on previous academic research), but I am not a lawyer. As for the LAPACK dependency: it might not be as big as it looks. LAPACK is pre-installed on OS X, either pre-installed or one apt-get away on Unix / Linux, and Windows binaries are readily available. But I don't want to flog a dead horse... – stephan Feb 09 '11 at 14:29
  • All I know is I never managed to get LAPACK working on Windows last time I tried to for an unrelated reason. (I don't remember the details.) Maybe if I figure out how, I'll keep my slow pure D routines around as the default and add a -version=LAPACK flag so if you have LAPACK installed and need the speed, you can get it. (The performance of this solver only matters in the relatively unusual case that N is large. For small N, naive algorithms are fast enough.) – dsimcha Feb 09 '11 at 18:21
  • The algorithms may be public domain, but the implementation can be copyright. – David Heffernan Feb 09 '11 at 19:19
  • 1
    @dsimcha It's really not hard to compile a handful of LAPACK routines, there's never many dependencies. You'd be crazy to try anything else. I get the Fortran, push it through `f2c -a` and then compile the C. Bob's your uncle! – David Heffernan Feb 09 '11 at 19:19