I have a somewhat peculiar Jx=R system that I need to solve. The matrix J is 2N -by- 2N. The first N rows have all entries filled. The next N rows are banded in two places, i.e. for the (N+k)th row, the (k-1,k,k+1)th and the N+(k-1,k,k+1)th entries are filled, and the rest are zero. There is absolutely no symmetry in the values, in the bands or the dense region, because I have a non-uniform grid. Also, J is the Jacobian matrix for a multi-variable Newton-Raphson implementation, and the entries of J get changed at each iteration. So, any method that imparts advantage due to repetition with different Rs but same J isn't of utility here.

Right now, I am solving using lapack's dgesv subroutine. I am happy with the accuracy of the results, but the speed is somewhat slow. I wanted to ask that given the nature of the matrix I am solving, would it be judicious to look for a faster method (which would require me to invest some time), or am I better off just continuing with dgesv.


  • 3
    $\begingroup$ It might help to show a picture of the sparsity structure. $\endgroup$ Nov 7 '19 at 15:08

I don't see any obvious way to factor $\mathbf J$ more efficiently. Even if you attacked the banded part with sparsity-exploiting solvers, you'd still have the burden of Schur-complementing it back onto the dense part, then factoring the dense part itself (so I would expect no real/asymptotic gains would be made).

Perhaps this is not answering the question you asked, but you might consider a jacobian-free-newton-krylov (JFNK) method since solving by your $\mathbf J$ is so problematic. They solve $\mathbf J \mathbf v= \mathbf r$ using Krylov techniques (most likely GMRES in your case), which means they don't need $\mathbf J$ explicitly, only its action upon a vector $\mathbf J \mathbf v$. And since the $\mathbf J$ matrix arises from differentiating some underlying objective function $\mathbf f(\mathbf x)$, the action $\mathbf J\mathbf v$ can be approximated by taking a suitable finite difference like $\mathbf J\mathbf v \approx \left[ \mathbf f(\mathbf x+\epsilon \mathbf v)-\mathbf f(\mathbf x) \right] / \epsilon$. (If the accuracy of the finite difference proves insufficient, you can also apply automatic-differentiation approaches to $\mathbf f$ instead).

See Knoll, Dana A., and David E. Keyes. "Jacobian-free Newton–Krylov methods: a survey of approaches and applications." Journal of Computational Physics 193.2 (2004): 357-397 for details.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.