Re: [eigen] Solving a tridiagonal linear system

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

indeed, I had to fix a couple of issue in TridiagonalMatrix, looks
like this class has never been used.

Actually, I'm not 100% sure about its current interface, in particular
regarding the SelfAdjoint option. Maybe it would be better to have two
different classes: TridiagonalMatrix, and
SelfAdjointTridiagonalMatrix.

gael

On Fri, Nov 26, 2010 at 12:43 PM, Andrea Arteaga <yo.eres@xxxxxxxxx> wrote:
> 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
>
>