Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 8 down vote favorite
1

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!

share|improve this question
possible duplicate of Recommendations for a usable, fast C++ matrix library? – GertVdE Oct 4 '12 at 15:24
Dear jep, welcome to the forum. Your question is very similar to the one here: scicomp.stackexchange.com/questions/351/… – GertVdE Oct 4 '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 4 '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 4 '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 5 '12 at 1:46
show 3 more comments

3 Answers

up vote 6 down vote

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.

share|improve this answer
+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 5 '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 5 '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 5 '12 at 14:11
up vote 3 down vote

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.

share|improve this answer
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 5 '12 at 0: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 5 '12 at 1: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 5 '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 6 '12 at 3:17
up vote 2 down vote

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.

share|improve this answer
1  
Michael - Please consider this a friendly warning that you need to follow the rule in the faq regarding disclosing affiliation. – Aron Ahmadia Oct 5 '12 at 0: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 5 '12 at 0: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 5 '12 at 0: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 5 '12 at 0:47
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 5 '12 at 1:33
show 8 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.