9

Most of my programming is one-off research codes in C for my own use. I have never distributed any code to other than close collaborators. I have developed an algorithm that I am publishing in a scientific journal. I want to provide the source code and perhaps executable code in the online supplement to the article. A colleague requested that I make a generalization to the algorithm which required me to write in C++ (ack!) and which requires that I solve small dense linear systems. If I succeed in getting a user base for the algorithm it will be partly because the entry bar to using it is low (like on the floor). Potential users won't install libraries, etc. in order to use the code. I want the code to be fully stand alone and unencumbered by any license at all. I might simply write my own solver by taking something out of Golub and van Loan but I'd rather use a vanilla solver that someone else has already written if there are any out there. Suggestions appreciated. Thanks!

Aron Ahmadia
  • 6,951
  • 4
  • 34
  • 54
jep
  • 193
  • 2
  • 7
  • possible duplicate of Recommendations for a usable, fast C++ matrix library? – GertVdE Oct 04 '12 at 15:24
  • Dear jep, welcome to the forum. Your question is very similar to the one here: http://scicomp.stackexchange.com/questions/351/recommendations-for-a-usable-fast-c-matrix-library – GertVdE Oct 04 '12 at 15:24
  • Library solvers tend to be complex and big for the sake of robustness, efficiency, and generality. If your problems are very small and reasonably well conditioned, I would suggest you to write your own mini-implementation. – Stefano M Oct 04 '12 at 21:52
  • @GertVdE, thanks for the quick response on this question. I'm uncomfortable linking to the "Recommendations..." question because both the question and the top answer are too general to provide any help in situations like these. If you'd like to discuss this further, I suggest we take it to the scicomp chat room. – Aron Ahmadia Oct 04 '12 at 23:29
  • @AronAhmadia: I think the only way to start settling some of these debates is to start implementing a computational science programming chrestomathy that is both language and library dependent. If the code is clear, and configuration issues can be taken care of (using a shell script, Chef, or Puppet), then debates about performance can be taken care of (or made concrete) by just running the code and timing it on a reference machine. Debates about clarity can be resolved (or at least, made more concrete) by looking at the code. Otherwise, we'll keep having the same arguments. – Geoff Oxberry Oct 05 '12 at 01:46
  • I don't think it's a duplicate, since the focus is quite different (lightweight/ease of installation/wide availability instead of performance/convenient api). Maybe the question can be edited to make the requirements more explicit. – Christian Clason Oct 05 '12 at 08:55
  • I read the other post before making mine and while they are clearly related I didn't think that mine was a duplicate for the reasons Christian Clason mentioned. Performance is not much of an issue for me. I have to solve one of these low dimensional systems per iterations in an iterative procedure. – jep Oct 05 '12 at 12:57
  • The bulk of the work lies elsewhere. I was hoping to find a rudimentary solver (10-100 lines of code say LU or gaussian elimination) that I could simply paste into my code. – jep Oct 05 '12 at 13:26

3 Answers3

7

I would suggest to exactly duplicate the Lapack interface to the function that you need, most probably you just need dgesv. That way people that have Lapack installed can simply link to it and it will just work. For people that don't have Lapack installed, you provide your own simple implementation of this function, or possibly implement it using Eigen or FLENS as others suggested.

In the Fortran land, the Lapack library is such a standard, that most people simply use it and that's it, instead of providing their own implementations.

Ondřej Čertík
  • 2,930
  • 18
  • 40
  • +1 Add to it the fact that most Linux distributions (at least Debian based) have binary packages in the repository and all vendor supplied math libs (MKL, SunPerf, ACML, ESSL etc.) carry it. You should always use standard libs as much as possible though if you're on Windows/Mac you might be better off with something C based as installing a free Fortran compiler (gfortran) on them is some amount of work or so I have heard. – stali Oct 05 '12 at 12:59
  • I have used lapack many times but I am not currently in fortran land. I expect that the statistical distribution of platforms my user base run on would be similar to that of the world at large meaning mainly windows, a smaller percentage of macs and an even smaller percentage of *nix. My experience with windows is minimal and I prefer to keep it that way. This is the reason I want a stand alone C++ code. I figure I'll have to provide some of my users with help getting the code to compile and run. I need to minimize the work required to do that. – jep Oct 05 '12 at 13:41
  • If your user base is Windows/Macs then you're better of with a simple C based (perhaps even your own) implementation. A package that is difficult to install or depends on 5 other libs, specially when there is no first class binary package repository (like Debian) available, will turn your users off for a long time. Remember most Windows/Mac users are used to one click install. Ease of use triumphs everything else. – stali Oct 05 '12 at 14:11
5

A very early mistake that many people make when getting started in scientific computing is assuming that you need to write all of your code in the same language. I think this is due largely to historical reasons, when it wasn't clear how to make compiled programs communicate with each other across even versions of the same compiler. That said, in this case, if you are going to be using C++ anyway, there are several very good C++ header-only template libraries that might fit your needs.

Since you are distributing your code for academic reasons, and you would like to embed a dense linear algebra solver into your code, I would strongly recommend that you consider Eigen. Eigen has been licensed under the Mozilla Public License and is a header-only library. This means that you can distribute Eigen with your code in source form (this does not impose any licensing restrictions on your code), and you will receive access to its general capabilities, including extremely efficient dense linear solvers. As GertVdE mentions, you have several other options.

Community
  • 1
Aron Ahmadia
  • 6,951
  • 4
  • 34
  • 54
  • I was hoping for a single file. I've done scientific programming for quite a while. I have mixed languages like C and fortran quite a bit but for this project I really just want one file containing all my source code. I suppose I could put a C solver in the C++ code which wouldn't be a big deal. Mainly I want to keep the code as simple as possible. LU with pivoting should be adequate. I'll look at Eigen. Thanks! – jep Oct 05 '12 at 00:57
  • @jep, you could also try cherry-picking the routines you need from CLAPACK if you really don't care at all about performance. – Aron Ahmadia Oct 05 '12 at 01:12
  • There are good reasons for writing all the dependent code in the same language, in particular, in HPC environments, you have weird compiler/linking issues and 32/64-bit interface issues. For example, how do I know the width of an integer for built-in libraries? How do I know for sure what compiler was used for a built-in library, and can I link against it with this other compiler? Having everything in one language simplifies many of these issues. And yes, there should be documentation provided by the cluster maintainers, but most of the time there isn't. – Victor Liu Oct 05 '12 at 22:57
  • @VictorLiu - The issues you are referring to are more tightly coupled to implementations than languages. Comment space is a poor place to get into a serious discussion, but I am happy to engage you in chat or elsewhere if you'd like me to expand my thoughts on this. – Aron Ahmadia Oct 06 '12 at 03:17
4

If you want a reliable solver for systems of linear equations I would recommend FLENS. It contains an exact re-implementation of LAPACK (it even reproduces the same roundoff errors as LAPACK if a single-threaded BLAS implementation is used). This is true for all FLENS-LAPACK functions (together with the utility functions about 100 routines).

FLENS is under a BSD License and therefore allows to be incorporated into proprietary products.

FLENS is header only and if you only need a subset of FLENS I can give you a stripped-down version containing only those functions you need. FLENS comes with its own reference BLAS implementation. But optionally your users can link against optimized BLAS libraries like ATLAS, OpenBLAS or GotoBALS. For large matrices this gives a performance gain of about 40% compared to Eigen.

And yes, Eigen also uses the LAPACK test suite to check their results. They do this for 3 functions (Lu, Cholesky and Eigenvalues/-vectors of a symmetric matrix). However, their computation of eigenvalues/-vectors of a non-symmetric matrix would fail the LAPACK test suite.

Disclaimer: Yes, FLENS is my baby! That means I coded about 95% of it and every line of code was worth it.

Michael Lehn
  • 315
  • 2
  • 7
  • 1
    Michael - Please consider this a friendly warning that you need to follow the rule in the faq regarding disclosing affiliation. – Aron Ahmadia Oct 05 '12 at 00:30
  • Sure, but you also could re-phrase your posts from 'I would strongly recommend that you consider Eigen' to something like 'there is for example Eigen'. In this case I delete my remarks about Eigen (although they are all proven to be true) including this one. – Michael Lehn Oct 05 '12 at 00:39
  • 1
    Your remarks about Eigen are not at issue here (although they seem off-topic to me). You are a primary developer of FLENS, if you are going to recommend it in an answer here, you must disclose your affiliation as developer of the project. – Aron Ahmadia Oct 05 '12 at 00:43
  • Ah, ok then. I thought was was implicitly clear by '... I can give you ...'. Is the disclosure in this form ok? – Michael Lehn Oct 05 '12 at 00:47
  • Yes, it's better to be explicit about these things than implicit. Thank you for your revision. – Aron Ahmadia Oct 05 '12 at 00:56
  • I think I can get by with LU decomposition with pivoting . I'm solving Ax=b and ultimately I'll want x. I've forgotten how the lapack lu decomposition works but I guess there are several functions needed. One that does LU and one that takes the LU decomposition and b and gives back x. Anyway how do I contact you to get code? I thought we weren't allowed to post e-mail addresses on here? – jep Oct 05 '12 at 01:00
  • At the bottom of the FLENS site you find a link to the author (which is me). And you also can ask Google. or you can send an email to the FLENS email-list. – Michael Lehn Oct 05 '12 at 01:03
  • @jep - If you'd like your email address or other information to be publicly listed, put it in the About section of your profile. More information here. – Aron Ahmadia Oct 05 '12 at 01:17
  • 2
    I just want to say thanks for doing this; I had similar plans to re-implement a large part of Lapack in C++. However, it seems that for most of the advanced (eigenvalue) routines, you simply defer to calling into Lapack, so it's a bit of false advertising to say that you re-implement the whole thing. On the other hand, I have actually ported the ZGEEV source to C++ in RNP, albeit some parts are still in 1-based indexing from auto-conversion. – Victor Liu Oct 05 '12 at 01:33
  • I really completely re-implemted dgeev with all its dependent routines (and that was quite some effort). Just look how I compile the examples at http://apfel.mathematik.uni-ulm.de/~lehn/FLENS/flens/examples/tut04-page04.html I am not linking with any external LAPACK. It just needs a C++ compiler. BTW: I you are looking for a complete C++ port of LAPACK look at MPACK (http://mplapack.sourceforge.net). – Michael Lehn Oct 05 '12 at 01:38
  • My primary intension of FLENS-LAPACK was to show that the combination C++/FLENS makes it easy to implement something like LAPACK without sacrificing performance. Its like an Fortran 90 implementation of LAPACK. In my opinion FLENS is rather a tool than a concrete library.
    Also look at this simple example (http://apfel.mathematik.uni-ulm.de/~lehn/FLENS/flens/examples/tut01-page08.html) showing how FLENS simplifies the implementation of getf2/getrf. Analogously to FLENS-LAPACK I am working on something similar to MAGMA with ViennaCl as GPU backend.     – Michael Lehn Oct 05 '12 at 01:59
  • @VictorLiu: About your comment "you simply defer to calling into Lapack" and "false advertising" is that resolved? For each LAPACK routine that I re-implemented I also have an (optional) interface/binding to an external LAPACK routine. I mainly use this for checking my results. And on the FLENS_LAPACK page I marked each LAPACK routine that was not completely re-implemented. So I really do not think that I did any wrong advertising. But I want to clarify any concern you have. – Michael Lehn Oct 06 '12 at 12:06