8
$\begingroup$

To give you some context,

I am currently implementing a simple finite element solver in Julia. I am getting run-times that are 70% of a Matlab code. (Both codes are essential equivalent in structure.) I've run some profiles on my Julia code and I find that most of the run-time is taken up by the linear system solver (essentially the "\" operator in $K \backslash F$ where $K$ is a sparse symmetric positive definite matrix and $F$ is a vector.)

I've tried looking up efficient solvers in Julia. I've come across the MUMPS package but I haven't used it yet. (I'm planning on giving it a go soon and will update this question when I do.)

More interestingly I came across this thread on the julia-users google-group which has a lot of content on how Matlab implements its (very efficient) solver.

My question is, what is the most efficient way of implementing a linear solver in Julia? Does the backslash operator choose the most efficient solution algorithm?

Improve this question
$\endgroup$
3
  • $\begingroup$ MATLAB is probably using UMFPACK (by Tim Davis) which you can compile and install yourself. I believe it is written in C so you'll have to figure out a way to call it from Julia. MUMPS is golden if you want distributed memory parallelism but there is little point in using it if the rest of the code is not parallel. $\endgroup$
    – stali
    Commented Jul 4, 2015 at 1:59
  • $\begingroup$ Thanks @stali. The rest of the code is not currently written to take advantage of parallel architecture. But once I've got a stable piece of code (there are other modules I'm currently working on) I will definitely be looking into parallelizing it. I've dabbled a little with MPI but I'm more curious to try Julia's own implementation of message passing. Perhaps at that point MUMPS will become an attractive option. $\endgroup$
    – Arjun Narayanan
    Commented Jul 4, 2015 at 6:30
  • 1
    $\begingroup$ You mean julia_time = 0.7 * matlab_time, so Julia is faster, or julia_time = 1.7 * matlab_time, so Julia is slower? $\endgroup$
    – Federico Poloni
    Commented Jan 31, 2016 at 9:28

2 Answers 2

Reset to default
9
$\begingroup$

In the context of finite element methods (and, especially, symmetric problems) the most common direct solution method is Cholesky factorization (plus following substitutions). MATLAB uses Tim Davis' CHOLMOD package to compute Cholesky factorization whenever the heuristics of backslash operator encounter a symmetric positive definite matrix.

In fact, Julia also interfaces Davis' CHOLMOD through its cholfact command. I have found that it is sufficient to call

u=cholfact(K)\F

where K is a sparse matrix.

Improve this answer
$\endgroup$
1
  • $\begingroup$ thanks for your reply. I suspected that Julia was attempting a Cholesky factorization when I saw calls to the cholfact routine while running some profiles recently. Julia seems to be doing this through a polyalgorithm, because I am only using the backslash operator, and Julia picks the type of factorization based on the type of the matrix. $\endgroup$
    – Arjun Narayanan
    Commented Jul 14, 2015 at 4:37
2
$\begingroup$

If you are using a Linux based OS, in a machine with Intel CPU and if you have access to Intel MKL (you can obtain a non-commercial copy for free), then the best you can have is PARDISO MKL, which is the Intel version of the Pardiso solver by Olaf Schenk. There is a Pardiso wrapper in Julia, and the performance is in average faster than UMFPACK and MUMPS. Furthermore, you can handle multiple right-hand-sides using BLAS 3 operations.

The syntax, however, is a bit more complicated than just a backslash, but you have enough examples to have a good idea how it works and you can take a look at the online documentation to see all the possible options. With the backslash, Julia will check the type of the matrix and it will use the best options, with PARDISO you need to specify the type of matrix yourself.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Note that Intel's interpretation of "non-commercial" no longer includes research done while employed anywhere. $\endgroup$
    – Christian Clason
    Commented Jan 31, 2016 at 9:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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