Re: [eigen] Solving a tridiagonal linear system

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

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

I think we should decide on an "owner" for each module, who's
responsible for making sure that we're ready to finalize the 3.0 API.

For the Eigenvalues module, it's either you or Jitse.

At this late stage, I would not make API changes unless there's a
really really serious reason to. People have started relying on our
APIs. Whether it's a template parameter or 2 separate classes seems
like mostly a matter of taste, so I'm inclined to leave things as they
are.

Benoit

>
> 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
>>
>>
>
>
>