# PETSc-like library for Julia

I want to build an application for Material Point Method (and probably other meshfree methods too) in Julia and I am looking for library for direct and iterative solvers that can help me with it. One of the most common library used for such applications is PETSc and I do know that it has bindings to Julia. But I am also looking for alternatives because of lack of documentation.

Some of the alternatives I found include: Elemental.jl (Again, no documentation), ArrayFire (seems promising especially for GPU computations).
Another alternative is to use CUDA.jl for writing the relevant kernels for the computation for GPU and using PETSc for multi-CPU support, but that means the app won't support other vendors. To be honest, this is the stack I am leaning towards. So, can anyone suggest libraries like PETSc for Julia that has ample documentation that works with GPUs and CPUs?

Julia is built in such a way that you will never see a full PETSc-like library, and that's on purpose. PETSc is not a single thing: it is an HPC library with some utility functions, linear solvers, nonlinear solvers, ... a whole can of soup that works within its own world but not outside of it.

Julia's package ecosystem is build on generics and composability. For this context it's best to explain it in an example. CUDA.jl provides the CuArray type for CUDA-based array operations. If you do x .+ y .* z with a CuArray, it uses the .ptx JIT compiler to generate the function (x,y,z) -> x + y * z element-wise and apply the map via CUDA. It has high level operations so that A*x is on the GPU, and A\x is on the GPU.

But is CUDA.jl incomplete compared to PETSc? Where are the iterative solvers like GMRES? Well, if someone writes a Julia library for iterative solvers that also uses generics, then it is "the" iterative solver for CUDA.jl. Krylov.jl for example has such CG/GMRES/etc. implementations. For ODE solvers, DifferentialEquations.jl. Nonlinear solvers, NonlinearSolve.jl, and the list goes on. The generics are all sufficiently good that if you use a GPU array, all computations in the package recompile to be on the GPU, and if it's distributed it won't do memory transfers without you asking, etc.

So this changes the question to the following: what are the good array types and what can they be used with? The SciML ecosystem works to be generally compatible with generic array types, so that's where you'll find extra functionality (which does a lot more than PETSc due to a larger labor force, separation of concerns, etc.). So the real glut is the array types. There are a lot of good ones, but for distributed it's less complete. Here's the general summary:

• CUDA.jl- CUDA-based arrays. Very good, very complete, uses codegen to allow easily generating new kernels so it supports most operations. Works with multi-GPU in many ways. NVIDIA GPU only.
• ROCArrays.jl- Similar to CUDA.jl but ROC-based support (AMD GPUs). Less complete than CUDA.jl, tends to miss a bit more linear algebra, but fairly usable.
• ArrayFire.jl- ArrayFire-based for "general GPUs". It loses performance compared to the other two because it doesn't include kernel generation for fusing GPU calls, but it's still a good choice. Much less used than the other two.
• DistributedArrays.jl - for distributed HPC stuff, but it's slow and missing a lot of linear algebra. It would rely on Elemental.jl for most linear algebra operations.
• MPIArrays.jl- a great prototype of a MPI-based array type which sadly isn't maintained much anymore.
• PETSc.jl- there is a wrapper being worked on, mostly to get the array type Mat sufficiently complete to do what I mentioned with the generic package usage.

As you can see, there needs to be more work in distributed array types, but that's what you're looking for. Not necessarily something PETSc-like, but something that just provides a solid distributed array to be used with downstream functionality.