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.