Re: [eigen] Solving a tridiagonal linear system

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

2010/11/26 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
> actually this class is in BandMatrix.h file which itself is in Core.
>
> Both BandMatrix and the TridiagonalMatrix specialization are currently
> useless since they support only storage but no operation are provided.
> Consequently I doubt that anybody use them.
>
> BandMatrix is only used internally by UpperBidiagonalization, and so
> I'm suggesting to make these classes internal as well.

Notice that I just made UpperBidiagonalization internal because it
wasn't used by anything else (and I expect the new SVD to use a
dynamic-upper-or-lower bidiagonalization like LAPACK does).

So I'm 100% OK to make UpperBidiagonalization and BandMatrix internal.
For good measure we could also do it for Hessenberg.

I remember that someone earlier objected against removing them (he was
using Hessenberg IIRC) but making them internal should be OK, it's
just a matter for him of doing internal::

Benoit

>
>
>
> gael
>
>
> On Fri, Nov 26, 2010 at 2:03 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>> 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
>>>>
>>>>
>>>
>>>
>>>
>>
>>
>>
>
>
>