Optimization on the multinomial manifolds of stochastic non-square matrices

Thanks for note! So I have an optimization problem with simple form but the decision variable is a large-scale matrix. My problem is similar to a existing problem here about multinomial manifolds and my objective function is simpler than the counterpart in that problem. The decision variable $mathbf T in {bf R}^{1300times 3100}$ is a stochastic matrix (each row sums to 1) with non-negative elements. $mathbf S in {bf R}^{600times 1300}$ is a known constant matrix with non-negative elements. $mathbf L in {bf R}^{3100times 9000}$ is a known constant sparse matrix. There are exactly one $1$ and one $-1$ in each column of $mathbf L $, so the column sum of $mathbf L $ is equal to $0$ for every column. Let $mathbf X$ denote the product of above three matrices, i.e. $mathbf X=mathbf Smathbf Tmathbf L$. The optimization problem is below.
$$
text{minimize}hspace{3mm}f(mathbf T)=sum_{i=1}^{600}sum_{j=1}^{9000}|mathrm X_{ij}|
$$
$$
text{subject to}hspace{19mm}mathbf T{bf 1}=bf{1}
$$
$$
hspace{39mm}mathrm T_{ij}>0
$$
The optimization function means I expect to obtain a matrix $mathbf X$ as similar to zero matrix as possible, so it can be replaced by minimizing Frobenius norm of $mathbf X$ or other similar forms. As far as I know, my problem resembles optimal transport in some ways. In addition, there may be some trick to reduce the amount of computation and cost of storage given the special nature of matrix $mathbf L$. I’m doubting whether the manopt toolbox in Matlab here is capable of dealing with this large-scale optimization problem (speed is not so important). I have thought that a python package named pymanopt here might work but the multinomial manifolds is not currently supported by pymanopt. I’m planning to employ a feasible algorithm or a computation environment to solve this optimization task. If there are ideas to make my mission possible, what would be the best approach here? Thanks very much for any advice or comments.