3
$\begingroup$

I need to find the root of a nonlinear system (which comes out of collocation, so I will change the order to test). I will likely have about 50-300 variables, and the Jacobian is going to be somewhere between 1/9 and 1/2 dense, and calculated using auto-differentiation.

What do people suggest for the linear system of this type? I know about trilinos, PETSc, and sundials, but don't know the other alternatives or have exposure to them. I realize that the large solvers are overkill due to the relative lack of sparsity in the example, but can it hurt?

My preference is for a C++-style interface, but I can deal with C if necessary.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
5
$\begingroup$

I am a PETSc developer so take my suggestion with a grain of salt, but I would use PETSc because

  1. the problem sizes are large enough that execution overhead should be minimal,
  2. you can trivially switch between various sparse and dense solvers (1/2 sparse should be treated as dense, but it might pay off to use a sparse solver for 1/9 at your sizes),
  3. a suite of globalization techniques quick experimentation, which is important for difficult problems, and
  4. you get monitoring and debugging tools so you can debug and understand convergence without writing your own scaffolding.
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks. So using eigen + cppad +petsc is a defensible combination? Is there a well maintained c++ binding for petsc? I found a few but they don't seem maintained or "real" $\endgroup$
    – jlperla
    Jul 21, 2014 at 14:48
  • $\begingroup$ 1. Yes. 2. The (object-oriented and relatively dynamic) C interface can of course be called from C++. The obvious mappings of PETSc's object model to C++ would give up functionality so we do not support them. PETSc's interface is flexible enough that you should have no trouble hooking it up to the class organization you have chosen for your application. (Arguably, no C++ interface could be less opinionated about how you organize your classes.) $\endgroup$
    – Jed Brown
    Jul 22, 2014 at 0:36
5
$\begingroup$

What do people suggest for the linear system of this type? I know about trilinos, petsc, and sundials, but don't know the other alternatives or have exposure to them.

Generally speaking PETSc, Trilinos, and KINSOL (from SUNDIALS) are the best-of-breed when it comes to scalability. From an ecosystem standpoint, PETSc seems most flexible, since it does not try to take over main or attempt to be a framework; many packages (SLEPc, libmesh, deal.II, PyClaw, MOOSE, TAO (before it was folded into PETSc), FEniCS) build on top of PETSc successfully.

Trilinos has more of a framework/ecosystem feel, as Trilinos consists of many packages; packages are still built on top of Trilinos, like FEniCS and deal.II. From mailing list traffic, the PETSc interface to FEniCS tends to get more use than the Trilinos interface.

KINSOL is the least flexible. It's basically a Newton solver or a fixed-point iteration solver, with a small number of interfaces to iterative linear solvers. It's probably the easiest to learn completely, but doesn't have the runtime flexibility that something like PETSc does. Choosing between KINSOL and PETSc, I'd pick PETSc; the code maintenance and quality alone are better, and PETSc is way more powerful.

I realize that the large solvers are overkill due to the relative lack of sparsity in the example, but can it hurt?

Probably not. There's some overhead due to PETSc initializing data structures, for instance, but I've only run into two people who have actually complained about it with good reason (and in one case, PETSc developers did a good job of reducing the PETSc memory footprint). JedBrown would know better; I'm just a PETSc user.

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.