Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
I have a sparse lower-triangular matrix $L$, and a right-hand side $b$, and I'd like to solve the linear system
$$L^T x = b$$
but without explicitly creating $L^T$. Ideally, I could write something like:
L.selfadjointView<Lower>().triangularView<Upper>().solve(b).
It seems that this would be valid if $L$ were dense, but it is sparse. I am trying to be memory efficient here, so I don't want to create $L^T$ explicitly. Is there any simple way to do this?
EDIT: I just noticed a glitch in your code. You used selfadjoint and triangular. Try this:
L.triangularView<Lower>().transpose().solve(b);
The rest of my answer may still be useful, so I won't delete it.
Presumably, $L$ uses the column-major storage order. If not, just swap row for column everywhere in this post.
$L^T$ can be represented by casting $L$ as a row-major sparse matrix. This will store $L^T$ in the upper half.
VectorX<double> b;
Eigen::SparseMatrix<double,Eigen::ColMajor> L;
Eigen::SparseMatrix<double,Eigen::RowMajor>::Map LT(L.cols(), L.rows(), L.nonZeros(), L.outerIndexPtr(), L.innerIndexPtr(), L.valuePtr());
VectorX<double> x = LT.triangularView<Upper>().solve(b);
This should avoid any copies as it just passes pointers to L's data. Although Eigen doesn't support RowMajor matrices for all its solvers, it does support the triangular solver. Also, this trick is very handy for when you want to use something like MKL which only supports RowMajor sparse matrices. If your matrix is symmetric...who cares if you're storing A or A^T!
