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*: "Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>*Date*: Mon, 15 Dec 2008 04:20:55 +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=8gYrzZpPTRW3TNdc0bgALuAXDG3ufeqIItq2JgbHLX4=; b=jndrNaRwCX9FERxZXgKAZ0SMpGC7lFyATBCpe5OrTgmNe8o42VisxFXPkQahGPmQP3 kp4I2Xp0TEAG2Z8uWPn/apTL9cvQ72rATx3vPZrBV9IYmT/XnKk5hP1S2qjfucSZ615u p2YnVXsGFzr58/L6Bkniu+76miMPNZIaYvNSw=*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=P9mRSVhbLYPpWR1qZTEPr3OY2OdBTjLQt8bqwYks33XGDBD18PKtZP5byTHB0hhHM1 MZm4GNP8aqrfGOfzw0lnfk5HBq4MJVQ7YghL5Aw8tf1j5E2b5HkIoNNMG7M80R0zDmwL EjbZcMPslAiasLCuldQnrGJCSh+VIWtfjhdfE=

Hi Mauro, Thank you for your interest in Eigen, and your patch. It's very nice to see how you managed to understand and appropriate yourself the code -- reminds me of when Gael joined. Your patch is actually very timely because we were just a few days ago discussing fixes and improvements in diagonal product and Gael just committed some fixes. Basically, the performance of diagonal product could be much improved as it is not yet vectorized; we were still wondering whether we should use a hardcoded function like Gael's CacheFriendlyProduct, or vectorize the diagonal Product expression. [[ Gael: i just saw that you put packet() methods in DiagonalProduct, so does this mean that at least some cases of diagonal product are already getting vectorized? ]] Now your patch looks very good so I suggest that we quickly commit it and hereafter we take tridiagonal product together with diagonal product in consideration when doing future improvements. > To get this I had to apply the following patch. What it does: > - define bits for upper/lower hessenberg matrices, and a Tridiagonal flag > - define a TridiagonalMatrix class, copied from DiagonalMatrix and > define TridiagonalProduct accordingly sounds good > - 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'. > With this changes it was faster than lapack for that task. No actual > number since I did not benchmark properly. In any case even if it was > on par it would still be a very good result at this stage. Impressive, as this part is not yet vectorized. > Maybe some of this changes could be applied to the tree. Sure, i want all of these improvements in the tree, and sooner rather than later! As soon as you can apply the minor improvements that I request below, and Gael says OK (he is more a specialist than me of the "products" files in Eigen). > I think that > the possibility to define new matrix types (block diagonal comes to > mind) could be very useful, so making this easy would be very nice. 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? If you mean new storage types: > One could use tridiagonals to test how difficult it is to add a new > storage type to 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... Now some comments on your patch: 1) please use 2 spaces for indentation 2) please try stay close to our "coding style" except if you really hate it, e.g. we don't put spaces around parentheses. (though there is sometimes some inconsistency between Gael and me already) 3) please note that my current e-mail address is jacob.benoit.1@xxxxxxxxx 4) give yourself a copyright line on any new file you create. + /* FIXME: I don't know what it should do + template<int LoadMode> + const PacketScalar packet(int row, int col) const + { + if (RhsIsTridiagonal) + { + ei_assert((_LhsNested::Flags&RowMajorBit)==0); + return ei_pmul(m_lhs.template packet<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col))); + } + else + { + ei_assert(_RhsNested::Flags&RowMajorBit); + return ei_pmul(ei_pset1(m_lhs.coeff(row, row)), m_rhs.template packet<LoadMode>(row, col)); + } + } + */ You did the right thing: let that commented out for now. As I said, we have yet to decide how to best vectorize these diagonal / tridiagonal products. The packet(row,col) method should return the 128-bit packet (so 4 floats or 2 doubles) starting at the given (row,col) (so it depends on the storage order). Gael: for diagonal we were hesitating between a C-style function or expressions for vectorization; now for tridiagonal product it seems like the expressions way would be a nightmare, no? So that's one more argument in favor of C-style functions. Plus, as I noted, for the diagonal product, a C-style function could work in-place.... not sure but i think that the tridiagonal case could still partly work in place in the sense that the amount of temporary buffer memory needed would be small compared to the memory used by whole matrices... + 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... Then m_coeffs could become m_diagonalCoeffs (and then we'd have m_upperCoeffs and m_lowerCoeffs). -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. - 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. The solution is to have a meta-selector, you can see tons of examples throughout Eigen, just grep for '_selector" or "_impl". Basically it's just using template specialization as a "early compile-time if" that really allows the compiler to only spend time in the path that's actually chosen. Don't worry if this is unclear, we can handle that for you. And, the line if ( (Flags & Tridiagonal) == Tridiagonal ) { can be written as just if(Flags & Tridiagonal) { it's actually less error-prone in case of future changes (more local, only 1 place to change)... Cheers, Benoit ---

**Follow-Ups**:**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Mauro Iazzi

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Usage of Eigen2 in a closed source application** - Next by Date:
**Re: [eigen] Usage of Eigen2 in a closed source application** - Previous by thread:
**[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/ |