Re: [eigen] Solving a tridiagonal linear system

[Andrea Arteaga <yo.eres@xxxxxxxxx>]

Many thanks, Gael, for your code. There is still a problem: it does not compile on my computer. I have the last version of eigen from bitbucket. The gcc (version 4.4.5) error is:

$ g++ tridiag_llt.cpp -I /home/spiros/Develop/eigen/ -o tridiag

In file included from /home/spiros/Develop/eigen/Eigen/Core:318,

                 from /home/spiros/Develop/eigen/Eigen/Dense:1,

                 from tridiag_llt.cpp:2:

/home/spiros/Develop/eigen/Eigen/src/Core/BandMatrix.h: In constructor ‘Eigen::TridiagonalMatrix<Scalar, Size, Options>::TridiagonalMatrix(typename Eigen::BandMatrix<Scalar, Size, Size, 1, ((Options & Eigen::SelfAdjoint) ? 0 : 1), (Options | Eigen::RowMajor)>::Index) [with Scalar = double, int Size = -0x00000000000000001, int Options = 16]’:

tridiag_llt.cpp:44:   instantiated from here

/home/spiros/Develop/eigen/Eigen/src/Core/BandMatrix.h:228: error: type ‘Eigen::BandMatrix<double, -0x00000000000000001, -0x00000000000000001, 1, 0, 17>’ is not a direct base of ‘Eigen::TridiagonalMatrix<double, -0x00000000000000001, 16>’

/home/spiros/Develop/eigen/Eigen/src/Core/BandMatrix.h: In member function ‘typename Eigen::BandMatrix<Scalar, Size, Size, 1, ((Options & Eigen::SelfAdjoint) ? 0 : 1), (Options | Eigen::RowMajor)>::DiagonalIntReturnType<-0x00000000000000001>::Type Eigen::TridiagonalMatrix<Scalar, Size, Options>::sub() [with Scalar = double, int Size = -0x00000000000000001, int Options = 16]’:

tridiag_llt.cpp:46:   instantiated from here

/home/spiros/Develop/eigen/Eigen/src/Core/BandMatrix.h:235: error: cannot call member function ‘typename Eigen::BandMatrix<_Scalar, Rows, Cols, Supers, Subs, Options>::DiagonalIntReturnType<Index>::Type Eigen::BandMatrix<_Scalar, Rows, Cols, Supers, Subs, Options>::diagonal() [with int N = -0x00000000000000001, _Scalar = double, int Rows = -0x00000000000000001, int Cols = -0x00000000000000001, int Supers = 1, int Subs = 0, int Options = 17]’ without object

I had the same problem yesterday trying to instantiate an object of TridiagonalMatrix with the SelfAdjoint option tuned on. I had a look of the actual error, and it seems to be here (BandMatrix.h, lines 222-228):

template<typename Scalar, int Size, int Options>


class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor>


{

	typedef BandMatrix<Scalar,Size,Size,1,Options&SelfAdjoint?0:1,Options|RowMajor> Base;


	typedef typename Base::Index Index;


  public:

	TridiagonalMatrix(Index size = Size) : Base(size,size,1,1) {}

The TridiagonalMatrix has some BandMatrix as base, then the BandMatrix type is setted ("typedef-ined") as Base type, but the 4th and the 5th template arguments are swapped.




The constructor of TridiagonalMatrix calls the constructor of Base, but the type Base is not exactly the base of the class, because of the swapped template parameters (well, I think so, but maybe you have a meaningful reason for doing that).

I changed the order of the template parameters to match the typedef and vice-versa and in both cases the code compiles, but during the execution an assert fails:

$ ./tridiag 

tridiag: /home/spiros/Develop/eigen/Eigen/src/Core/PlainObjectBase.h:493: void Eigen::PlainObjectBase<Derived>::_init2(typename Eigen::internal::traits<Derived>::Index, typename Eigen::internal::traits<Derived>::Index, typename Eigen::internal::enable_if<(Eigen::internal::dense_xpr_base::type::SizeAtCompileTime != 2), T0>::type*) [with T0 = long int, T1 = long int, Derived = Eigen::Matrix<double, 2, -0x00000000000000001, 1, 2, -0x00000000000000001>]: Assertion `rows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows) && cols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols)' failed.

Aborted

I'm not very expert with this type of exception, so I just report this issue to you, but I'm very unsure about my deductions.

If the SelfAdjoint-ness of the matrix is neglected, so if I instead write:

TridiagonalMatrix<double,Dynamic,0> t(n);

  t.diagonal().setRandom();

  t.sub().setRandom();

  t.super() = t.sub();

then your code works very well. Since I can do so without to get troubles, I'm very happy with your code.

Cheers.

Andrea Arteaga

2010/11/26 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