5
$\begingroup$

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...

Improve this question
$\endgroup$
1
  • $\begingroup$ This is only borderline ontopic. Your problem seems to be with technology, not computer science concepts. I think the question would have been more at home on Computational Science. However, you received an (apparently helpful) answer that highlights concepts, so I'll let the question remain here. $\endgroup$
    – Raphael ♦
    Jun 23, 2015 at 12:27

2 Answers 2

Reset to default
4
$\begingroup$

The main point of Krylov-Newton is that it does not require computing and storing the whole Jacobian matrix. It only requires the ability to apply the matrix to a vector.

The stardard approach is to do a little bit of pencil and paper work to figure out how to express the action of the derivative without actually computing and storing the whole matrix. If this is not possible, one can try to figure out a compressed representation of the matrix.

For problems where the matrices are dense, but the long range interactions are low rank (I think many radiative transfer problems would fall in this category..?), fast multipole or H-matrix type methods tend to be very effective.

It's difficult to give further advice without knowing more specifics of the problem.

Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ Also, this question would probably get more and better answers at scicomp.stackexchange.com $\endgroup$
    – Nick Alger
    Jun 23, 2015 at 7:45
  • $\begingroup$ Thanks for the response Nick. Yes I think you are right, I need to implement a proper Newton-Krylov method. I was being lazy and trying to plug in linear solves from libraries such as ViennaCL to my own Newton solver. Perhaps I can still use the GPU acceleration of ViennaCL for the matrix vector products in the Newton-Krylov method... $\endgroup$
    – psilouette
    Jun 23, 2015 at 7:51
  • $\begingroup$ What this answer ignores is the solve. Unless you have 1) a good algorithm for the action of the approximate inverse (FMM, H-matrix, quasi-Newton, etc.), or 2) a good matrix-free preconditioner, then your Newton solve will fail if you naively use the action of the Jacobian. $\endgroup$
    – Matt Knepley
    Sep 27, 2018 at 14:47
  • $\begingroup$ @MattKnepley Agreed it's better to use a preconditioner if you have one, but properly implemented inexact newton-krylov with proper trust region or line search still works ok even without. It becomes roughly equivalent to nonlinear conjugate gradient. Worst case scenario you do one krylov iteration per Newton step, in which case it reduces to gradient descent. $\endgroup$
    – Nick Alger
    Sep 27, 2018 at 17:47
  • $\begingroup$ The statement is only true if you have an inherently well-conditioned problem when restricted to the trust region, such as what you would get from the Hessian of an optimization problem. This is definitely not true for general PDEs. $\endgroup$
    – Matt Knepley
    Oct 2, 2018 at 17:02
2
$\begingroup$

You could also consider nonlinear variants that do not need Jacobian inversions. For Jacobians that look like Hessians, there are quasi-Newton algorithms. You can also try Anderson Acceleration and the related nonlinear GMRES. We try to cover the taxonomy of these solvers in Composing Scalable Nonlinear Algebraic Solvers

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.