New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 230 - ARPACK support
ARPACK support
 Status: NEW None Eigen Unclassified Sparse (show other bugs) 3.0 All All --- enhancement Nobody

 Reported: 2011-03-22 00:16 UTC by Gael Guennebaud 2014-08-27 10:21 UTC (History) 3 users (show) brix dharmon jim

Attachments
Patch for initial version of ArpackSelfAdjointEigenSolver (20.01 KB, application/octet-stream)
2012-07-24 06:28 UTC, David Harmon
no flags Details
Test driver for ArpackSelfAdjointEigenSolver, -larpack needed to build. (2.24 KB, application/octet-stream)
2012-07-24 06:29 UTC, David Harmon
no flags Details
ARPACK support (31.30 KB, application/octet-stream)
2012-10-07 01:21 UTC, David Harmon
no flags Details
New patch fixing some tiny bugs including EIGEN_STATIC_ASSERTs (3.68 KB, patch)
2013-06-18 08:06 UTC, Kolja Brix
no flags Details | Diff

 Gael Guennebaud 2011-03-22 00:16:47 UTC One way to get eigen decomposition on sparse matrices is to wrap ARPACK++, here is a start: http://forum.kde.org/viewtopic.php?f=74&t=84238&p=192268#p192198 and here is a copy of schizofrenikh's message from the forum: #include #include #include class TMatrixForArpackpp { private: Eigen::SparseMatrix< double > &M; public: TMatrixForArpackpp( Eigen::SparseMatrix< double > &_M ) : M(_M) { ; }; void MultMv( double *v, double *w ) { Eigen::Map< Eigen::VectorXd > v_map = Eigen::VectorXd::Map( v, M.rows() ); Eigen::Map< Eigen::VectorXd > w_map = Eigen::VectorXd::Map( w, M.cols() ); w_map = M.marked()*v_map; } }; void Lanczos_symm( Eigen::SparseMatrix< double > &M, std::vector< double > &eigenValues, std::vector< Eigen::VectorXd > &eigenVectors ) { int n = M.cols(); TMatrixForArpackpp M_arp( M ); ARSymStdEig< double, TMatrixForArpackpp > prob( n, n-1, &M_arp, &TMatrixForArpackpp::MultMv, "SM" ); //"SM" to get smallest eigenvalues prob.FindEigenvectors(); //prob.FindEigenvalues(); int m = prob.ConvergedEigenvalues(); eigenValues.resize(m); eigenVectors.resize(m); eigenValues.assign( prob.RawEigenvalues(), prob.RawEigenvalues()+m ); for( int i=0; i. Is this ok? (Rows is a typedef of the nbr of rows, and can be Dynamic as well) 2) If an unsupported Scalar type is used, is it better to have a runtime assertion tripped or compile-time error? The advantage I see of a runtime assertion is that I can explicitly say "That type is not supported", but a compile-time error will be more cryptic when it doesn't find the necessary template implementation. 3) How do I call the Eigen SimplicialLLT code from a header file? I need a linear solver for generalized eigen problems. I tried to use SimplicialLLT, but calling from a header file gave a bunch of errors. Are there known ways around this? 4) Is there an interface for cleanly allowing the user to choose between linear solver types that I can copy? I would like to give the option of built-in solver, CHOLMOD, UMFPACK, etc. I'm sure there's other code that does this, and I can just copy what it does if someone points me to it. 5) Do we trust what the user gives us? For example, in the generalized problem Av=\lambda B v, different choices must be made based on if B is PD, PSD, etc. Do we trust what the user says or do some error checking? The choices I see are a) always trust what the user says, and if we get garbage, then so be it. b) never trust the user, check the rank of B. c) check the rank of B if NDEBUG is not defined. That's all my questions for now. Gael Guennebaud 2012-07-21 22:22:27 UTC 1) It should simply be Matrix, the size of the input matrix is assumed to be large, so Dynamic. 2) Compile-time error. To this end you can use a static assertion in the general not supported template instantiation: STATIC_ASSERT(false,NUMERIC_TYPE_MUST_BE_REAL); If NUMERIC_TYPE_MUST_BE_REAL is not appropriate, we can add one in Core/util/StaticAssert.h 3) I don't really understand your issue. You must include Eigen/SparseCholesky, and then you can use SimplicialLDLT. I guess that ARPACK requires a function pointer. ou have to write it as a wrapper around SimplicialLDLT. Ideally, this wrapper function could be a template function supporting any solver, and the solver class should be specified by the user with a default to SimplicialLDLT. 4) Half answered in 3). You could make it a template parameter of ArpackEigenSolver: template > class ArpackEigenSolver { ... }; 5) We trust what the user say, but we need to report to the user when ARPACK failed. Our standard is to report errors via the .info() method which returns Success or NumericalIssue, NoConvergence, InvalidInput. David Harmon 2012-07-22 19:02:50 UTC Thanks for the info. 2) Small issue, but just calling EIGEN_STATIC_ASSERT(false, NUMERIC_TYPE_MUST_BE_REAL); gives me a "NUMERIC_TYPE_MUST_BE_REAL is not member of static_assertion" error. I explicitly put in static_assertion::NUMERIC_TYPE_MUST_BE_REAL and the error went away. I didn't seen anywhere else in Eigen doing this, so I must be doing something wrong. 3) Dumb mistake on my part. I wasn't following the Eigen convention of letting the caller be responsible for including necessary header files. if I put include Eigen/SparseCholesky in my ArpackSelfAdjointEigenSolver header file, crazy things happen, but if my caller simply says #include #include and uses the arpack solver, then everything is fine. 4) No wrapper is needed. Arpack computes stuff, returns it to you, you perform the linear solve, then call again, repeat until convergence. I seem to recall the old linear solver module using a template to choose the solver (Cholmod, Umfpack, etc.). I'll follow that format. I'll have the symmetric and generalized symmetric code ready in the next week or two. It's a lot of repeated code, but since the dense solver does it, I can create two separate classes for SelfAdjoint problems and GeneralizedSelfAdjoint problems. Gael Guennebaud 2012-07-23 20:58:28 UTC "lot of repeated code"-> there must be something wrong. Perhaps you could share your code even if it is not ready at all. Better to get early feedbacks to avoid useless work! David Harmon 2012-07-24 06:25:55 UTC I'm attaching a patch of the class ArpackSelfAdjointEigenSolver that solves for the specified number of eigenvalues/eigenvectors for both dense and sparse matrices, with either float or double Scalar type. Error checking is included as well. I've put the code in unsupported/, mirroring the dense setup, except naming the module "ArpackSupport" instead of "Eigenvalues". I'm also attaching a test driver that generates random dense/sparse symmetric matrices and computes the eigenvalues/vectors, and then prints out || A v - \lambda v || for each one. David Harmon 2012-07-24 06:28:47 UTC Created attachment 288 [details] Patch for initial version of ArpackSelfAdjointEigenSolver David Harmon 2012-07-24 06:29:43 UTC Created attachment 289 [details] Test driver for ArpackSelfAdjointEigenSolver, -larpack needed to build. David Harmon 2012-10-07 01:19:58 UTC I've had a chance to work on this some more and have added support for the generalized eigenvalue problem. I've tested it for a variety of problems and have been using it in my own research. I'd like it to go "live" as an unsupported module so others can start using/testing it. Currently it only works with self-adjoint problems, and does not support the arpack functionality of finding all eigenvalues around a given value (still trying to get it to work). I created a new file in unsupported/Eigen/ArpackSupport. The source is in unsupported/Eigen/src/Eigenvalues/. The code is pretty well commented (I think). The patch with the new files will follow, or should I commit to the repository directly? David Harmon 2012-10-07 01:21:37 UTC Created attachment 298 [details] ARPACK support David Harmon 2012-12-04 02:08:41 UTC I've had some people contact me after applying the patch manually themselves and using this extension. Could somebody with write access to the repository please just check it in so people can start using it properly? If there is some other way I should submit the patch, please let me know. Thanks. Kolja Brix 2013-06-18 08:06:41 UTC Created attachment 350 [details] New patch fixing some tiny bugs including EIGEN_STATIC_ASSERTs Kolja Brix 2013-06-18 08:09:02 UTC Unfortunately, there were still some tiny bugs in the ARPACK wrapper, i.e. wrong a include file name. These issues are fixed with the new patch including the EIGEN_STATIC_ASSERT issue. Thanks to David and Gael for proving the ARPACK wrapper, special thanks to Gael for helping with fixing the EIGEN_STATIC_ASSERT issue. Gael Guennebaud 2013-06-18 09:45:28 UTC https://bitbucket.org/eigen/eigen/commits/c5057358673a/ Changeset: c5057358673a User: Kolja Brix Date: 2013-06-18 09:44:40 Summary: Bug 230, fix compilation issues and wrong static assertions Jim Garrison 2013-08-16 02:11:12 UTC The code in ArpackSelfAdjointEigenSolver.h says "by default, the upper triangular part is used, but can be changed through the template parameter." However, the documentation for the regular (dense) SelfAdjointEigenSolver says "Only the lower triangular part of the input matrix is referenced." To keep consistency, should ArpackSelfAdjointEigenSolver be changed to instead look at the lower triangular part of the matrix by default? Gael Guennebaud 2013-08-20 13:26:23 UTC I agree with you that the default should be to read the lower part, but that's unclear how to do it, and the default linear solver is SimplicialLLT which already reads the lower part. So, it might be that the doc is wrong. Jim Garrison 2014-07-05 08:51:31 UTC Also, I have found that ArpackGeneralizedSelfAdjointEigenSolver fails to compile if used with a dense matrix type -- it only seems to work with sparse matrices. It would be nice if it worked with any type of matrix. Actually, even more important: it would be nice if we did not have to store the entire sparse matrix in memory at all if it can be easily computed on the fly. For instance, I have a function that takes a row index, then calls a callback for each nonzero column entry in that row. I would love to be able to use this directly with Eigen::VectorXd. Gael Guennebaud 2014-07-09 14:17:54 UTC You can already write a small class inheriting from SparseMatrixBase and implementing an InnerIterator to loop through the non-zeros. A good example is SparseView which iterate over the non-zeros of a dense matrix. Jim Garrison 2014-08-26 06:16:11 UTC I've been working on some improvements to Eigen's ARPACK support at . I made some minor cleanups to ArpackSelfAdjointEigenSolver.hpp, and I also have created ArpackEigenSolver.hpp. Currently ArpackEigenSolver works only for the type std::complex. It should be easy to also extend to work with std::complex. Working with regular float and double will be a bit tricky and involve a bit of refactoring, since ARPACK takes a slightly different set of arguments when called for real types. There are a few suboptimal things with the current version, which can be noticed by doing a diff against ArpackSelfAdjointEigenSolver.hpp. I tried to keep the differences to a minimum, since much of the code has been copy-pasted from there and tweaked. Also, the patch series contains a patch similar to the one provided for bug #858, a patch which I made before noticing that report. Gael Guennebaud 2014-08-26 13:26:02 UTC I had a quick look to your changesets, and I don't really understand why the sigma variable has to be of type Scalar while it is really real. Regarding the new class, would it make sense to introduce a common base class to reduce duplicated code? Jim Garrison 2014-08-26 18:40:53 UTC Gael, you are right about sigma. I tried to make it consistent in both versions but failed. In ArpackSelfAdjointEigenSolver sigma should be type RealScalar, but in ArpackEigenSolver it should be type ComplexScalar. I am very unfamiliar with hg so I am unsure how to modify my revision history to fix this. (Actually, I remember hearing once that hg discourages such modifications, so I'm unsure what to do instead.) For background, there are three classes problems as far as ARPACK is concerned (see http://www.caam.rice.edu/software/ARPACK/UG/node32.html): (i) real symmetric, (ii) real non-symmetric, and (iii) complex. Each is supported for both single and double precision. Notice that for complex types, the user is supposed to call the same routine for both Hermitian and non-symmetric problems. There's really two possible solutions as far as Eigen's API is concerned: make them aware of this difference, and tell them to use ArpackEigenSolver any time the matrix is complex (even if it is Hermitian); or, have ArpackSelfAdjointEigensolver abstract over this difference by accepting complex matrices, taking a triangular view, and then in the end converting the eigenvalues to real before returning. If the latter option is chosen, we would also need to abstract over the "which" field with an enum, because e.g. if one wants the smallest algebraic eigenvalue using ArpackSelfAdjointEigensolver, the "which" string would need to be "SA" for real problems and "SR" for complex problems. Jim Garrison 2014-08-27 04:26:40 UTC FWIW, scipy chooses the latter option I mentioned: if eigsh (sparse hermitian) is called with a complex (Hermitian) matrix, it calls znaupd and simply discards the (small) imaginary part of each returned eigenvalue. Gael Guennebaud 2014-08-27 10:21:07 UTC yes I'd also go the scipy way. Regarding code factorization, I was thinking about something like our CholmodSupport module (https://bitbucket.org/eigen/eigen/src/a1f5fed0d734c2993474868950af76a6889ba9f2/Eigen/src/CholmodSupport/CholmodSupport.h?at=default) where a CholmodBase class implements all the common code to interface cholmod while the derived classes (CholmodSimplicialLLT, CholmodSimplicialLDLT, CholmodDecomposition) simply provide the end-user specific interface. That being said, I don't known if there is much to factorize or not as I did not look carefully to your wrapper code.

 Note You need to log in before you can comment on or make changes to this bug.