What cuda sovling libraries and APIs are recommended to use for BBD(Bordered block diagonal) linear equations?

The given BBD sparse linear equations are as the following eq.(1), where coefficient matrix is a 2x2 block sparse matrix, x and b are both vectors.

The principle of solving the eq.(1) is as follow(eq.(1)(2)(3)(4) are shown in the picture):

  1. Factor sparse coefficient matrix of eq.(1) into a lower and an upper triangular matrix, as coefficient matrix of eq.(2).
  2. Forward substitute variable y of eq.(3) to get solution y.
  3. Backward substitute variable x of eq.(4) to get solution x.

Could you help recommend available libraries and APIs to solving eq.(1) according to the principle mentioned above?