3

For dense matrices, the following code solves x^T A = b^T just fine.

Matrix3d A;
RowVector3d bT, xT;

A << 1, 2, 3, 
     0, 5, 6, 
     0, 0, 9;

bT << 1, 2, 3;

xT = A.triangularView<Upper>().solve<OnTheRight>(bT);
printf("(%g, %g, %g)", xT(0), xT(1), xT(2));

I cannot continue this approach to sparse matrices, however.

SparseMatrix<double> spA = A.sparseView();
spA.triangularView<Upper>().solve<OnTheRight>(bT);  // COMPILE ERR!
spA.triangularView<Upper>().solve<OnTheRight>(bT.sparseView());  // COMPILE ERR!

The compile errors are

no matching function for call to ‘Eigen::SparseTriangularView<Eigen::SparseMatrix<double, 0>, 2>::solve(Eigen::RowVector3d&) const’
no matching function for call to ‘Eigen::SparseTriangularView<Eigen::SparseMatrix<double, 0>, 2>::solve(const Eigen::SparseView<Eigen::Matrix<double, 1, 3> >) const’

candidate is:

template<class OtherDerived> typename Eigen::internal::plain_matrix_type_column_major<OtherDerived>::type Eigen::SparseTriangularView::solve(const Eigen::MatrixBase<OtherDerived>&) const [with OtherDerived = OtherDerived, MatrixType = Eigen::SparseMatrix<double, 0>, int Mode = 2, typename Eigen::internal::plain_matrix_type_column_major<OtherDerived>::type = Eigen::internal::plain_matrix_type_column_major<T>::type]

I could not find the answer in the documentation, can anyone figure out how to do this?

EDIT SparseTriangularView::solve accepts neither OnTheLeft nor OnTheRight as template argument, but I just tried neglecting the argument and it seems to compile. My guess is that it is a missing feature and have reported it to the developers as so. If they confirm, I will post their response as an answer.

user2646234
  • 197
  • 1
  • 8

1 Answers1

2

This is indeed a missing feature, but you can easily workaround by transposing everything:

xT.transpose() = spA.transpose().triangularView<Lower>().solve(bT.transpose());

or, if you directly deal with column vectors:

x = spA.transpose().triangularView<Lower>().solve(b);
ggael
  • 28,425
  • 2
  • 65
  • 71
  • Is there a similar trick for the lack of QR.solve(b) where QR is a SparseQR object factorizing A? QR.matrixR().triangularView() is not square for general A, so I guess I have to take the QR.rank() x QR.rank() top-left corner of it before solving? – user2646234 Jun 27 '14 at 21:08
  • The only reason I am not just using QR.compute(A.transpose()) is because I am already using QR.compute(A) for other purposes and the mathematics seems so easy: xT A = bT is solved with AP = QR in two steps. Find yT in yT R = bT P. Compute xT = yT QT. – user2646234 Jun 27 '14 at 21:24