Make PETSc solvers available

[stevengj]

Rather than rolling our own sparse-matrix library, it seems like it would be much better to use PETSc, which is the de-facto standard here.

PETSc would give us access to a huge number of sparse-direct back-ends (#2632), iterative solvers, sparse eigensolvers (#1573) via the SLEPc library, optimized implementations (#942), multiple output formats, and parallel support (currently completely missing and quite nontrivial to implement well).

Re-inventing the wheel here is a little like rewriting LAPACK.

[stevengj]

If you only change ABI on the decade scale, that is easy enough to keep up with. I've maintained plenty of C libraries (e.g. you are using my FFTW library), and there is no reason for things like enums or basic function signatures to change unless there is a major API upheaval, and in that case we would have to modify our wrappers anyway.

We can certainly parse the header file to generate the corresponding C interfaces (we even have a whole Clang package to do just that using the C compiler). But if you don't maintain a modicum of ABI stability, even this is extremely fragile — we cannot distribute pre-generated code since it might not match the Petsc on the user's system, and it is dangerous not to regenerate at runtime in case the library ABI has changed out from under them.

Moreover, I don't want to just munge your header files and present the raw C interface to Julia users. That would be do-able with Clang in a fully automated way. I want a high-level interface that integrates nicely with the rest of Julia. But this has to be hand-written, and if you have any major API changes it will have to be hand-modified to keep up. As long as you don't have ABI changes for no reason (i.e. not unless there is a corresponding API change), then it is not really any more work to keep with ABI changes than it is to keep up with API changes.

[stevengj]

At the very least, we have to be cautious at runtime about the possibility that PetscReal and PetscScalar have changed due to a library recompile (since this changes the ABI even without changing Petsc versions), which is why it would be good to be able to query these dynamically.

[jedbrown]

@stevengj I humbly submit that your comment under-estimates the size of the design space for scalable algebraic solvers by some orders of magnitude. FFTW is a great library, but it restricts itself to doing essentially one well-studied thing very well. The number of papers per year on new algorithmic variants, often with order-of-magnitude impact on problem classes appearing in various applications is one way to measure the difference in scope.

PETSc's base APIs rarely change, but the advanced APIs will probably continue to change occasionally until our field stops finding better algorithms and different forms of exploitable structure. I understand that the Julia interface will need updates in these cases (if the Julia interface tries to be "complete"; i.e., if your target audience includes people developing new methods of the sort that might "normally" go inside PETSc) and I think that compiler type checking of some sort is extremely valuable in this setting. If you can use clang to verify that the calls make by the Julia interface is compatible with the interface in the header, perhaps even optionally, I think it would be very worthwhile. This is the sort of type checking that Cython offers to petsc4py, which also exposes a more pythonic (and complete) interface to PETSc. (The first generation of petsc4py used generated bindings and another layer to provide a pythonic interface.)

We'll put Barry's patch into the next maintenance release (3.4.3) so you can use it.

[ViralBShah]

The approach I usually prefer its to wrap the C api using clang if possible giving a raw interface and then writing a Julian api on top of that. Users can then work at either layer.

[stevengj]

@ViralBShah, the thing to realize is that the Petsc C interface is huge, and the vast majority of it consists of things that most Julia users will never call directly if we have a native high-level interface. Until we get precompiled modules, wrapping the entire thing with Clang will impose a big time penalty for little benefit on anyone loading PETSc.jl .... the tiny number of very advanced users who need some inner Petsc functionality can always use ccall to call that routine directly.

(An analogous thing would have been for my PyCall library to first use clang to wrap the entire Python C API, which almost no one would use directly from Julia, and then write the higher-level interface on top of that.)

@jedbrown, don't your API changes consist mainly of new APIs for new structures etcetera, rather than changing the old ones? The thing to realize is that writing wrappers is different from using the C API directly. In writing wrappers, I only call the C function once in the code, and in that case having an automatic translation of the C header is at best a minor convenience; I have to look carefully at your documentation (and occasionally your source code) in order to figure how best to call it, and typing in the parameters is the least of the difficulties.

[stevengj]

(And there is absolutely no technical excuse for ever changing the value of an enum, which is where this discussion started.) (Aren't your advanced API changes mainly the addition of new APIs, rather than changing the old APIs?)

[stevengj]

Also, writing wrappers is different from using the C API directly. When one writes wrappers, one only calls each C API function once. Compared to the effort required to look in the Petsc manual and carefully scrutinize the documentation in order to determine the best way to call it, the effort of typing in the parameters is minimal; automatic translation of the headers is only a minor convenience.

And keeping up with occasional API changes in wrappers is not really a burden, as long as you aren't changing existing APIs gratuitously (e.g. renumbering enums on every release).

[stevengj]

But since people seem to feel strongly about this, I'll run Clang on the Petsc header files and include an auto-generated C API in addition to the high-level wrappers.

[jedbrown]

@stevengj I can think of at least one case where arguments to a function changed because the old interface implied that the caller would need to have done a non-scalable operation. Occasionally, a function is renamed or arguments are changed to be consistent with similar operations elsewhere. (We value consistency and simplicity above long-term binary stability.) One reason for an enum value to change is if making the values bitwise ORable makes the interface simpler and more powerful, or if it is converted to dynamic registration (the macros become string constants, users/plugins can register more). These changes aren't common, but we're not going to promise to never make such changes.

If the wrapper is distinctly higher level than the underlying interface, how is the wrapper going to maintain a more stable API while still being feature-complete? As you say, PETSc's interface is large, but most users start out using only a few basic routines, and as the demands of their application grow, they incrementally use more. If the wrapper only provides a conservative subset of the functionality, chances are that the interfaces it touches won't ever change, but a large fraction of users will need to dip into the lower level interface.

I look forward to seeing how the Clang-generated API is done. That approach may be very useful for us in other projects.

[BarrySmith]

@stevengj Note that PETSc already contains string representations for all enum types that you can use in Julia so you don't need to re-represent them all in your Julia bindings. For example

/*E KSPCGType - Determines what type of CG to use

Level: beginner

.seealso: KSPCGSetType() E_/ typedef enum {KSP_CG_SYMMETRIC=0,KSP_CG_HERMITIAN=1} KSPCGType; PETSC_EXTERN const char const KSPCGTypes[]; const char const KSPCGTypes[] = {"SYMMETRIC","HERMITIAN","KSPCGType","KSPCG",0}; /_@C PetscEListFind - searches enum list of strings for given string, using case insensitive matching

Not Collective

Input Parameters:

• enumlist - list of strings to search, followed by enum name, then enum prefix, then NUL

• str - string to look for

  Output Parameters:

• value - index of matching string (if found)

• found - boolean indicating whether string was found (can be NULL)

  Notes: Not for use in Fortran

  Level: advanced @/ PetscErrorCode PetscEnumFind(const char const enumlist,const char str,PetscEnum value,PetscBool found)

If these don't suite your need they can always be extended/modified for more general usability.

[ViralBShah]

Steve, you are right that auto generated headers may take forever to load.

[ihnorton]

It's about 6 seconds, which is not great, but better than I expected. (should be a bit faster with proper exclusion criteria so that it only pulls in the public API). See https://github.com/ihnorton/Clang.jl/commit/1489ab3e05227ca60bc4d785b5afc25446fd707b#commitcomment-3937825

edit: also, this doesn't help with the issue of PetScScalar and PetScReal sizes (unless you run Clang.jl with exactly the same config arguments as when the library was built, which just begs the question)

[ViralBShah]

6 seconds is actually ok.

[timholy]

If the work to reduce loading times of Julia itself can also be extended to packages (presumably by building them along with Julia itself, at the user's discretion), then such issues will become much less important.

And I'd agree that 6 seconds for a single package is tolerable, as long as there are few such packages you need to load to get work done. It doesn't take many such packages to yield 30s load times, which is why (currently) it's important to think about keeping things lean when possible.

[stevengj]

Note that on top of this 6 seconds you will have to add the time to load an MPI package, if you want to do any kind of extensive parallel work with PETSc (the basic PETSc interface is usable without directly calling MPI, but in my experience real applications often need a few direct MPI calls outside the PETSc interface).

In the long run, it is really essential that Julia gain the ability to cache precompiled copies of modules, but looks like we have a ways to go. (Precompiled Julia is on the horizon with #3892, but it seems a big jump from there to precompile modules in such a way as to still be able to load them in arbitrary order at runtime.)

[hiemstar]

Hi everyone. Great initiative interfacing PETSc with Julia. Since all posts on the subject are more then a year old, I was wondering if any progress has been made since?

[stevengj]

I'm afraid that I haven't had a chance to work on it.

[hiemstar]

Thanks for the update Steven.

On Mon, Dec 1, 2014 at 7:40 PM, Steven G. Johnson notifications@github.com wrote:

I'm afraid that I haven't had a chance to work on it.

— Reply to this email directly or view it on GitHub https://github.com/JuliaLang/julia/issues/2645#issuecomment-65170762.

[jiahao]

@JaredCrean2 and @stevengj have started work on PETSc.jl, so let's move further discussion over there. PETSc.jl still needs a lot of work since PETSc is such a massive library, and we can use all the help we can get!

[stevengj]

Note that PETSc.jl seems to be progressing nicely. @JaredCrean2 and @ihnorton put together auto-generated Clang wrappers of the entire PETSc API, so the low-level C API is completely wrapped and available, and we are slowly progressing towards a nicer high-level API built on top of that.

Also note that we are wrapping/linking the double, single, and complex-double versions of PETSc simultaneously, so that you don't need to recompile to switch precisions or from real to complex.

[BarrySmith]

On Nov 18, 2015, at 11:51 AM, Steven G. Johnson notifications@github.com wrote:

Note that PETSc.jl seems to be progressing nicely. @JaredCrean2 and @ihnorton put together auto-generated Clang wrappers of the entire PETSc API, so the low-level C API is completely wrapped and available, and we are slowly progressing towards a nicer high-level API built on top of that.

Also note that we are wrapping/linking the double, single, and complex-double versions of PETSc simultaneously, so that you don't need to recompile to switch precisions or from real to complex.

Nice, how are you managing the fact that in our C code we use the same symbol for, for example, VecAXPY() independent of he data type. Are you somehow automatically mapping to different symbols? Maybe we should use that same technology of PETSc as well?

Does this mean I can remove the small amount of .jl code I put in PETSc a few years ago as a test prototype?

Barry

— Reply to this email directly or view it on GitHub.

[JaredCrean2]

Hello Barry, We build 3 versions of Petsc, and add a static parameter to all the Julia objects that is the datatype of the library they belong to. Multiple dispatch does the rest for us. We do have to run the auto-generation process 3 times, generating 3 versions of the low level wrappers.

I think its safe to remove the .jl code in Petsc at this point.

[tkelman]

the small amount of .jl code I put in PETSc a few years ago as a test prototype

Fun fact: that little bit of code was the first place I ever heard of Julia. I didn't actually start using it for another 1.5-2 years, but oh well.

[stevengj]

@BarrySmith, to clarify, when you call a C function in Julia, you supply both the symbol and the library name, via ccall((symbol,libname), ....args...). So, there is no obstacle to calling the same symbolic name for functions in different shared libraries (here, 3 different Petsc libraries).

(The auto-generated wrappers on top of this mean that the user never needs to see the low-level ccall, and gets the full type-checked PETSc API, automatically dispatched to the correct underlying library.)

The analogous thing in C would be to dlopen all three Petsc libraries, call dlsym to lookup a symbol in the desired library, and call the function via that function pointer. But because C lacks overloading, doing this in a type-checked way would be awkward.