I am solving a large (~1e5 equations & unknowns) set of nonlinear equations using Newton-Raphson iterations. Currently I am using the GPU accelerated Krylov methods implemented in ViennaCL to solve the linear system to get the update increment. I am solving the system on a single 24 core, 64GB workstation with a NVIDIA Quadro K4000 GPU. However as the number of unknowns exceeds 1e5 and/or the Jacobian matrix gets more dense, there is not enough memory on the GPU to fit the compressed Jacobian matrix on the device. Using the CPU cores allows the Jacobian matrix to fit into memory, however the compute time is very long.
I would like to use a cluster of the above mentioned workstations to solve either the nonlinear system or the linear update increment system, however I am not sure how to go about decomposing the system of equations into pieces that can be tackled via a distributed memory approach. The equations include radiative transfer, and therefore it is not easy to decompose the domain geometrically as the Jacobian matrix is dense.
Does anyone know about distributed memory approaches to solving large dense linear equation systems?
All help is greatly appreciated!
PS I have looked into the parallelised nonlinear solvers implemented in PETSc, however this is a clunky c library and one must write there entire problem in terms of the PETSc interface which I would like to avoid if possible. I would be more interested in understanding the details as to how a distributed memory parallelised nonlinear or linear solver works...