# Partially Banded Matrix

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

• It might help to show a picture of the sparsity structure. – Wolfgang Bangerth Nov 7 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).