Re: [eigen] patch for tridiagonal matrices in eigen |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] patch for tridiagonal matrices in eigen*From*: "Mauro Iazzi" <mauro.iazzi@xxxxxxxxx>*Date*: Mon, 15 Dec 2008 12:24:14 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:received:message-id:date:from:to :subject:in-reply-to:mime-version:content-type :content-transfer-encoding:content-disposition:references; bh=hQc5xo9nyf0VN+mA7Vebodb6rQOxhoVYr1w9vLLRCAY=; b=VFg+e1m+EeaxqRZS5ar14w6de49v8idXlwfaaL1cbafL4Li9g+nHhpdhJQdNudCDB7 rhz8aqeXmi8NHF7EciEtFFF/UXxQFt9ZVOIleyFlGo6lkWGEOXt/2LJvqYxfjLc2hk48 HplSiLjMcOd8DS65XBVUxGV8M5YGm2kj2VOgY=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=message-id:date:from:to:subject:in-reply-to:mime-version :content-type:content-transfer-encoding:content-disposition :references; b=nlxmG8hwYtDR2U9cJYamLWl+JwKvLQ7v5+WIKuWUMG7Jwpjohf+UQwEO+UzKshsMU/ 1/9OiuLyOeuGYuYNuUvt4/yf+1eimHPvcWE+m5Fo6xf+mBHjedkkSdnWQ6hqpcTBDEli F37XK/lUYuOMrA3vwENwX2g45n7N9bh06iBwo=

Hi Benoit, thanks for the comments. I'll apply the changes as soon as Gael comments on them and I adapt them to the correct coding style (do you have guidelines written somewhere?). >> - in SelfAdjointEigenSolver: use the _MatrixType::SquareMatrixType >> for storing the eigenvalues, so that it can be instantiated with e.g. >> _MatrixType == TridiagonalMatrix >> - skip tridiagonalization if the matrix is already tridiagonal (not >> very happy with this, since it's a runtime check if compiler does not >> optimize). > > I comment below, this should really use a meta-selector instead of a 'if'. or maybe there's a better way. > Do you mean adding new storage types, or just new expression / product > types as you just did? > If you mean new product types, sure it's always good to make Eigen > more extensible, but looking at Product.h it already looks quite > extensible to me. Maybe you can elaborate after your experience? I meant to define a new matrix class and specializing expressions, like product. It was very easy to do, so it's already very extensible, but I had to change the headers for that. Ideally I could just define without having to put new code in the library. > Actually I don't think that it should be difficult at all to add, > e.g., a storage version of TridiagonalMatrix. > The best solution, if possible, would be to just say that the storage > version is just a partial specialization of TridiagonalMatrix where > the three vector types are actual storage vector types (not abstract > expressions). I haven't looked closely enough at that to be sure that > it makes sense. That would probably require a flag forcing the coeffs > vectors to be stored by value. Then some neat typedefs like > TridiagonalMatrixXd so it's mostly transparent to the user. Or default > to storage by value, I don't know... I expressed myself badly, sorry. However I think that should work. > Now some comments on your patch: > + protected: > + const typename CoeffsVectorType::Nested m_coeffs; > + const typename UCoeffsVectorType::Nested m_ucoeffs; > + const typename DCoeffsVectorType::Nested m_dcoeffs; > > Hm... apparently you use the letter D for "down" ? Very slippery as it > could mean diagonal. > So I suggest more explicit names: Diagonal,Upper,Lower... very wise. I'll do that. > -const unsigned int Upper = UpperTriangularBit; > -const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit; > -const unsigned int Lower = LowerTriangularBit; > -const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit; > +const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit; > +const unsigned int StrictlyUpper = Upper | ZeroDiagBit; > +const unsigned int Lower = LowerTriangularBit | LowerHessenbergBit; > +const unsigned int StrictlyLower = Lower | ZeroDiagBit; > > here I don't really follow. These are parameters for part(). So unless > you want to add Hessenberg support in part (which would be a good > idea) no need to change that. If you want to add Hessenberg support, > you need to make the corresponding changes in Part.h. I saw that the product helper checked for diagonals using the flags, so I did the same for tridiagonals. > - TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, > computeEigenvectors); > + if ( (Flags & Tridiagonal) == Tridiagonal ) { > + for (int i=0;i<n-1;i++) { > + diag[i] = matrix(i, i); > + subdiag[i] = matrix(i, i+1); > + } > + diag[n-1] = matrix(n-1, n-1); > + m_eivec = EigenVectorMatrixType::Identity(n, n); > + } else { > + TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, > computeEigenvectors); > + } > > Here, the compiler will have to go through the whole logic of BOTH > cases, which means uselessly long compilation times. Yes the compiler > will always optimize the "if" itself (even in debug mode) so it's not > a runtime speed issue, nor a code size issue, it's really a > compilation-time issue, and we're extremely careful about that. understood. Maybe switching to MatrixType::SquareMatrixType also inside Tridiagonalization<> would allow to specialize it for the tridiagonal case. This would avoid storing and checking the flags and also the manual copy of the vectors, since we could use the ones that are already there (m_diagonalCoeffs, etc...). > And, the line > if ( (Flags & Tridiagonal) == Tridiagonal ) { > can be written as just > if(Flags & Tridiagonal) { I think that would incorrectly get Upper/Lower Hessenberg as Tridiagonals... I've seen (Flags & Diagonal) == Diagonal in other places for the same reason. cheers, mauro ---

**Follow-Ups**:**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Benoit Jacob

**References**:**[eigen] patch for tridiagonal matrices in eigen***From:*Mauro Iazzi

**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Usage of Eigen2 in a closed source application** - Next by Date:
**Re: [eigen] patch for tridiagonal matrices in eigen** - Previous by thread:
**Re: [eigen] patch for tridiagonal matrices in eigen** - Next by thread:
**Re: [eigen] patch for tridiagonal matrices in eigen**

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