Different behaviors of solve, solve_triangular, and solve_cholesky (possibly a bug)

[shiqingw]

I want to solve linear equations like AX = B where A is a n-by-n nonsingular square matrix, and B is a n-by-m matrix. When I tried to make use of the structure of the coefficient matrix A, I found that the behaviors of solve, solve_triangular, and solve_cholesky are different. It seems that solve does what I want, but solve_triangular and solve_cholesky only solve the first column of B.

Here is a simple example:

#include <iostream>
#include <xtensor/xarray.hpp>
#include <xtensor/xio.hpp>
#include <xtensor-blas/xlinalg.hpp>

int main()
{   xt::xarray<double> A {{4,0,0},
                          {0,4,0},
                          {0,0,4}};
    xt::xarray<double> L {{2,0,0},
                          {0,2,0},
                          {0,0,2}};
    xt::xarray<double> B {{1,0,0},
                          {0,2,0},
                          {0,0,3}};

    std::cout << xt::linalg::solve(A, B) << std::endl;
    std::cout << xt::linalg::solve_triangular(A, B) << std::endl;
    std::cout << xt::linalg::solve_cholesky(L, B) << std::endl;
}

The output is:

{{ 0.25,  0.  ,  0.  },
 { 0.  ,  0.5 ,  0.  },
 { 0.  ,  0.  ,  0.75}}
{{ 0.25,  0.  ,  0.  },
 { 0.  ,  2.  ,  0.  },
 { 0.  ,  0.  ,  3.  }}
{{ 0.25,  0.  ,  0.  },
 { 0.  ,  2.  ,  0.  },
 { 0.  ,  0.  ,  3.  }}

Can we align solve_triangular and solve_cholesky with solve?