Re: [eigen] Solving a tridiagonal linear system |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Solving a tridiagonal linear system*From*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Fri, 26 Nov 2010 12:58:24 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:in-reply-to :references:from:date:message-id:subject:to:content-type :content-transfer-encoding; bh=E0jvoYgmaNKrW/1SnfP79d/v5NWqPQIzb01/oGjWZOo=; b=p8V9qKZSev6xmOJ+0n1IzLiNzvC2ewMtr2Faey0KNHH5HRYO/ouAmzKrbcMW0QUwxf 8U5KdfmolPoWDX6OF2EQscN2TU7pXOnsnqt7cuOB3C59GePer6v2UHYZPVLHM2xe75dx tZd7rq8wVoVK0uGZaeW+cUKXUoM/wU5XLx2hk=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; b=jOguMC9++J33XyfmiMOP2AufqsA6LyLpOjji8EIk7rMrc8pAzPqghPUkIJk8USxeNN nb8JJm8zM3GTjFay8VKexkJUU32IFtBT9tjpyJ8QFd0NQA+i2wTj8fU466qMu3enLwHF FLPPM9QFHCfGv/VwjXnqMgK1hBph/ngEkDIbk=

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

**Follow-Ups**:**Re: [eigen] Solving a tridiagonal linear system***From:*Benoit Jacob

**References**:**[eigen] Solving a tridiagonal linear system***From:*Andrea Arteaga

**Re: [eigen] Solving a tridiagonal linear system***From:*Gael Guennebaud

**Re: [eigen] Solving a tridiagonal linear system***From:*Andrea Arteaga

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Solving a tridiagonal linear system** - Next by Date:
**Re: [eigen] Solving a tridiagonal linear system** - Previous by thread:
**Re: [eigen] Solving a tridiagonal linear system** - Next by thread:
**Re: [eigen] Solving a tridiagonal linear system**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |