3

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.

excalibur1491
  • 141
  • 1
  • 1
  • 7
  • 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)? – GoHokies May 11 '17 at 12:37
  • 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. – GoHokies May 11 '17 at 12:44
  • 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 – excalibur1491 May 11 '17 at 12:51
  • (efficient) sparse linear algebra comes with a steep learning curve. you'll need to get acquainted with these concepts and understand how the various libraries work if you want to get good performance out of your implementation. if that is not an option, you can code up something close to a black-box solver using Eigen's sparse matrix interface, and tweak the most obvious solver parameters (ask a colleague with domain knowledge for some help on this one). – GoHokies May 11 '17 at 13:35
  • Would you mind giving a pointer to the other circuit simulator in Python? Is it ahkab? – AlexE May 11 '17 at 14:54
  • 1
    @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? – GoHokies May 11 '17 at 16:05
  • With superlu, make sure you use a permutation that optimizes the sparsity of the factors ( see superlu doc), it has a huge impact on performances ( depending on the structure of your matrices) – BrunoLevy May 11 '17 at 17:10
  • Are you sure the python code is actually solving the linear system with scipy.sparse.linalg solvers? If you see the word "pyamg", then probably it is solving the linear system iteratively with algebraic multigrid (PyAMG is a python algebraic multigrid package), which is a completely different approach from sparse direct factorization. – Nick Alger May 11 '17 at 18:20
  • Possible duplicate of Libraries for solving sparse linear systems – Jakub Klinkovský May 11 '17 at 18:26
  • @NickAlger: You are absolutely right! They use PyAMG to solve it in an iterative way (not sure what that means, I am reading wikipedia page for "Multigrid method"). Here is a tutorial of PyAMG. Starting from that code I got a quick demo of solving a matrix very similar to the one I want of 10000x10000 and it solve almost immidiatly. Is there an equivalent library in C++? – excalibur1491 May 12 '17 at 00:35
  • @GoHokies See comment above where I reply to Nick Alger. You were rights, I should consider an iterative method. Any suggestion of a fast c++ library for that? – excalibur1491 May 12 '17 at 00:36
  • @AlexE The software I was looking at was this: https://pypi.python.org/pypi/Circuitscape/ it's used for habitat conservation, but deep inside it uses an electrical circuit (cf compute.py and compute_base.py). What they do would be one of the features of my software, so I am doing the same thing + many others, but I still need that little piece. – excalibur1491 May 12 '17 at 00:38
  • @GoHokies I added in the original post the call to spsolve. I did not change any option of superlu. Now that I think of it, my matrices are symetric, so probably I should change that option. Ideally I would like my software to work with different libraries, so that the user can choose depending on their system, so I am still keen to get superlu to be fast, but I guess I am more interested on those iterative methods you suggested as they seem to give better results for my type of matrices. – excalibur1491 May 12 '17 at 00:44
  • @jlk : The question you marked as duplicate does not ask for the library to be in C++ like I do. And in fact, the 3 suggestions for iterative methods are NOT in c++. Please remove the tag – excalibur1491 May 12 '17 at 07:17
  • @excalibur1491 Except for Eigen and Armadillo, none of the libraries suggested so far is a pure C++ library. And Armadillo is mostly a wrapper for external C libraries. So in the end, any fast library will be usable from C++ one way or the other. Although the questions are formulated differently, they yield basically the same answers, so I think the flag is appropriate. – Jakub Klinkovský May 12 '17 at 07:43
  • @jkl I don't think I can call Matlab or Mathematica from C++ and I'd rather not do those horrible tricks to call Python from C++. Those are the 3 lamguages for iterative methods in the other question. A C or C++ libs are fine, other languages, not really. – excalibur1491 May 12 '17 at 07:46
  • @excalibur1491 Python, Matlab and Mathematica are under "Interactive Environments" in the other question/answer, not "Iterative Methods". Python uses C solvers internally, which you can use directly. (Matlab and Mathematica have the same principle but are proprietary, so you can't really use their routines.) In any case, if you ignore "Interactive Environments", you will find your answers as a subset of the answers to the other question. – Jakub Klinkovský May 12 '17 at 07:59
  • @jkl Sorry, my bad. What I meant is, so far in this question I have been offerece two iterative method solutions: PyAMG and scipy (the cg function). Both are in Python. I don't know how to use their inside routines from C/C++. If you can tell me that, that would solve all my problems. If possible, I'd rather avoid having to go through the Python C++ interface, but rather call the inside function in PyAMG or Scipy directly. – excalibur1491 May 12 '17 at 08:27
  • 1
    @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. – Jakub Klinkovský May 12 '17 at 15:27
  • @jlk Thanks for the comment. All I wanted is to know which one of those libs can be used from C++ and contains iterative methods. PETSc looks promissing, I am trying a few examples now. Thanks a lot! – excalibur1491 May 13 '17 at 07:13
  • Comments are not for extended discussion; this conversation has been moved to chat. – Paul May 17 '17 at 20:20

5 Answers5

6

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.

BlaB
  • 1,157
  • 6
  • 17
  • 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) – excalibur1491 May 11 '17 at 12:03
  • 1
    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. – BlaB May 11 '17 at 12:38
  • 1
    As a word of warning, "parallel" and "direct solve" aren't always things that play well together. – origimbo May 11 '17 at 19:45
  • 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... – BlaB May 13 '17 at 12:44
3

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.

Hsien-Ming Ku
  • 129
  • 1
  • 9
2

I would recommend SuiteSparse or cuSPARSE.

Juan M. Bello-Rivas
  • 3,994
  • 16
  • 24
0

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.

davidhigh
  • 3,127
  • 14
  • 15
0

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.

cfh
  • 586
  • 2
  • 6