2

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.

Thanks!!

pkaran
  • 21
  • 2

1 Answers1

4

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.

rchilton1980
  • 4,862
  • 13
  • 22