Re: [eigen] Solving a tridiagonal linear system

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Indeed, currently we only provide a storage class for band and
tridiagonal matrices and the operations like products, triangular
solving and more general solver are missing.

However writing an efficient dedicated solver is quite straightforward
since there is no need to consider special optimization such as
blocking and/or vectorization.

In the meantime, I've attached an example of a Cholesky solver for
selfadjoint tridiagonal matrices which uses the subdiagonal only.

Adapting it for a tridiagonal LU factorization should be easy.

gael

On Fri, Nov 26, 2010 at 2:03 AM, Andrea Arteaga <yo.eres@xxxxxxxxx> wrote:
> Hi.
> First of all thank you for this fantastic library!
> I have a TridiagonalMatrix resulting from a cubic natural spline
> interpolation problem. Now I have to solve a system with this matrix and a
> rhs vector, but I don't find any method to do this in an efficient way.
> Actually the only method which seems to exist is doing A.toDenseMatrix()....,
> but in this way the optimization due to the tridiagonal pattern gets lost (I
> suppose). So, how can I do that?
> Thanks.
> Andrea Arteaga

#include <iostream>
#include <Eigen/Dense>

namespace Eigen {
namespace internal {
  
  template<typename Scalar,int Size, int Options> bool tridiag_llt_inplace(TridiagonalMatrix<Scalar,Size,Options>& mat) {
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef DenseIndex Index;
    const Index size = mat.cols();
    for(Index k = 0; k < size; ++k)
    {
      RealScalar x = real(mat.diagonal().coeff(k));
      if (k>0) x -= abs2(mat.sub().coeff(k-1));
      if (x<=RealScalar(0)) return false;
      mat.diagonal().coeffRef(k) = x = sqrt(x);
      if (k+1<size) mat.sub().coeffRef(k) /= x;
    }
    return true;
  }
  
  template<typename Scalar,int Size, int Options, typename Vector> void tridiag_llt_solve_inplace(TridiagonalMatrix<Scalar,Size,Options>& mat, Vector& vec) {
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef DenseIndex Index;
    const Index size = mat.cols();
    for(Index k = 0; k < size; ++k)
    {
      if(k>0) vec.coeffRef(k) -= vec.coeff(k-1) * mat.sub().coeff(k-1);
              vec.coeffRef(k) /= mat.diagonal().coeff(k);
    }
    for(Index k = size-1; k >= 0; --k)
    {
      if(k<size-1) vec.coeffRef(k) -= vec.coeff(k+1) * conj(mat.sub().coeff(k));
                   vec.coeffRef(k) /= mat.diagonal().coeff(k);
    }
  }
}
}

using namespace Eigen;
int main()
{
  int n = 10;
  TridiagonalMatrix<double,Dynamic,SelfAdjoint> t(n);
  t.diagonal().setRandom();
  t.sub().setRandom();
  MatrixXd d = t.toDenseMatrix();
  VectorXd x0(n), x1;
  x0.setRandom();
  x1 = x0;
  
  internal::tridiag_llt_inplace(t);
  tridiag_llt_solve_inplace(t,x0);
  x1 = d.selfadjointView<Lower>().llt().solve(x1);
  
  std::cerr << x0.transpose() << "\n" << x1.transpose() << "\n";
  
  return 0;
}