3
$\begingroup$

I am working on a project with electrical circuits, where I am trying to compute the voltages at all the nodes of an electrical circuit. I know that the electrical circuit is a perfect grid, so each node only touches at most 8 other nodes. Which means that I endup trying to solve a system: L v = i where L is at least a 1000x1000 matrix where each row has only 9 non-zeros. So, it's really really sparse.

I have tried solving it using SoPlex, LaPACK and SuperLU (these last 2 through Armadillo). But all are too slow. On a 10000x10000 the best I have is 18s.

I know there is another software that does the same task as mine, written in Python that uses Scipy (scipy.sparse.linalg to be precise) that is ridiculously fast (can solve those systems in less than one second).

Is there a library that is equivalent to Scipy or a way of porting Scipy to C++? I need to write the software in C++, for other reasons...

EDIT: My code to call SuperLU/LaPACK through Armadillo is simply: voltages[i] = spsolve(laplacians[i],iflow[i],"lapack"); Or voltages[i] = spsolve(laplacians[i],iflow[i],"superlu"); No options have been given before.

Improve this question
$\endgroup$
21
  • $\begingroup$ you say scipy.sparse.linalg.spsolve is much faster? what is the storage format of your sparse matrix? do you re-order rows/columns to minimize fill-in? what's the sparsity pattern (and memory usage) of your LU factors? do you use a preconditioner? have you considered using an iterative method instead of a direct one (GMRES)? $\endgroup$
    – GoHokies
    May 11, 2017 at 12:37
  • $\begingroup$ also, the scipy sparse solvers use either SuperLU or UMFPACK so, for a $\rm 10k \times 10k$ sparse matrix with $\tilde 10^{6}$ nonzeros, a good C++/SuperLU implementation should match (or at least come very close to) scipy's performance. $\endgroup$
    – GoHokies
    May 11, 2017 at 12:44
  • $\begingroup$ Oh, that is good to know. I amcomparing to someone else's code so I dont know if they used something tricky. There is someyhing called pyamg dangling around the code that I'm not sure what they use it for. How would I reorder the columns to minimize fill in? My diagonal is always nonzero, so should I put other nonzeros around it? Not sure what LU factors are..., Nor a preconditioner. Sorry, first time trying to do this kind of things $\endgroup$
    – excalibur1491
    May 11, 2017 at 12:51
  • 1
    $\begingroup$ @excalibur1491 could you share a minimal working code sample for your Armadillo/SuperLU implementation? have you changed any of the default SuperLU options before calling spsolve? $\endgroup$
    – GoHokies
    May 11, 2017 at 16:05
  • 1
    $\begingroup$ @excalibur1491 I feel that your original question has been answered extensively here and by the other question's answers. It's increasingly difficult to satisfy all the off-topic branches, because now you even mix recommendations for methods with their implementations (i.e. solvers). Nevertheless, to answer one more question, PETSc has multigrid preconditioners, including an interface for the other packages listed under "Large Distributed Iterative Solver Packages" in the other question. $\endgroup$
    – Jakub Klinkovský
    May 12, 2017 at 15:27

5 Answers 5

Reset to default
6
$\begingroup$

I second the idea of using Eigen, which is pretty efficient, but also very simple to include. If you need a lot more performance, you could try to use PETSc or Trilinos. They are very powerful libraries to store and solve sparse systems, they allow for a large number of iterative or direct solvers and are compatible with MPI for added performance. However, I think that they would be highly efficient even in serial.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Hi, I will take a look at those. Thanks! Could you develop more on that MPI comment? I keep comming across that, but I have no clue what it is (I am fsr from being an expert of scientific computing) $\endgroup$
    – excalibur1491
    May 11, 2017 at 12:03
  • 1
    $\begingroup$ MPI means message passing interface, I am far from being an expert on the topic, but briefly, it is a library that can be used to allow communication between computing processes. It is a way to allow for parallel computing on small to very large systems (the processes can be on a single computer or across nodes, like in a cluster). What this means is that both these libraries are made with supercomputing in mind and they can utilize the maximum of the ressources available. $\endgroup$
    – BlaB
    May 11, 2017 at 12:38
  • 1
    $\begingroup$ As a word of warning, "parallel" and "direct solve" aren't always things that play well together. $\endgroup$
    – origimbo
    May 11, 2017 at 19:45
  • $\begingroup$ I agree, but I feel like the best place to find an adequate parallel direct solver would be within PETSc and Trilinos. I mean, to me these libraries are at the top of the game when it comes to linear solvers in parallel... $\endgroup$
    – BlaB
    May 13, 2017 at 12:44
3
$\begingroup$

the size is 1k, it is not large-scale. I think many sparse direct solvers can be used. First of all, you need to make sure that your matrices are symmetric (or not)? , then you can choose suitable packages for handling them, such as PARDISO, MUMPS, UMFPACK, SuperLU and so on. If you also want to try the iterative solvers with preconditioning, you can try to use ILUPACK, ARMS, etc.

By the way, can you share your matrices with us? Recently, we collect some test matrices arsing from real applications for studying the numerical algorithms. We plan to establish a database of test matrices.

Improve this answer
$\endgroup$
2
$\begingroup$

I would recommend SuiteSparse or cuSPARSE.

Improve this answer
$\endgroup$
1
  • $\begingroup$ perhaps useful to note that Eigen can be coupled with both SuiteSparse and MKL Pardiso $\endgroup$
    – GoHokies
    May 11, 2017 at 16:07
0
$\begingroup$

Try Eigen? ... basic advice, I know.

Eigen has added several sparse solvers in the last time, and also offers benchmarks with established methods such as UMFPACK or SUPERLU against which it seems to perform quite well.

Improve this answer
$\endgroup$
0
$\begingroup$

10000x10000 isn't that large for such a sparse matrix, so direct solvers should still work fine. In addition to the Eigen recommendation already made, here's another option: download Intel MKL and then use the MKL PARDISO direct solver, which is in my experience much, much faster than SuperLU.

Improve this answer
$\endgroup$

Your Answer

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

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