Triangular linear system with triangular right hand side in python

I have to solve a linear system of equations with multiple right hand sides, A*X=B, where both, A and B are (upper) triangular, real, square matrices. The size is about 200 by 200. Is there a fast method for this in python/numpy?


I was considering looping over the columns,
n=A.shape[0]
X=zeros((n,n))
for i in range(n):
X[:i+1,i]=solve_triangular(A[:i+1,:i+1],B[:i+1,i])



But this does not use fast matrix-matrix operations.


I could also do all right hand sides simultaneously, X=solve_triangular(A,B), but this does not take into account the triangular structure in B.


Finally, I could invert A and multiply with B, X=inv(A)@B, but inverting matrices is usually discouraged from.