I have a system which is T * X = Y, T is tridiagonal. However, for numerical stability, I’ve used the bidiagonal factorization

(B’ * D * B) * X = Y, B bidiagonal with a unit diagonal, D diagonal matrix.

Obviously the easiest way to solve this would be 2 calls to gtsv, but this results in a lot of inefficiency in terms of useless data reads (the 1 main diagonal and 0 off-diagonal). Furthermore, it seems to me that there might be a simpler algorithm to solve bidiagonal matrices than PCR. Is there any algorithm that works only on Bidiagonal matrices that I should know about? Also, is there any other approach one might take to solve this system? Cheers.