0
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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!

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.