[LinearSolvers] new solver for (symmetric) generalized eigenproblems using the spectra library #8440

armingeiser · 2021-03-05T21:41:25Z

Description
Inclusion of the spectra lib and a interface to its more efficient solver for generalized eigensystems of real symmetric matrices.

What is Spectra

A header-only C++ library for large scale eigenvalue problems
License: MPL2 (same as Eigen)
From the spectra github page:

Spectra stands for Sparse Eigenvalue Computation Toolkit as a Redesigned ARPACK. It is a C++ library for large scale eigenvalue problems, built on top of Eigen, an open source linear algebra library.

Spectra is implemented as a header-only C++ library, whose only dependency, Eigen, is also header-only. Hence Spectra can be easily embedded in C++ projects that require calculating eigenvalues of large matrices.

In a very first test it seems to be faster then the eigen_eigensystem (except if only 1-2 eigenvalues are searched)

Eigenfrequency of solid model
with MKL support (both solvers use PardisoLU internally)

n-dofs	n-eigenvalues	time `eigen_eigensystem` [s]	time `spectra_sym_g_eigs_shift` [s]	speed up
~110k	1	2.77	5.0	0.5x
~110k	2	6.3	4.3	1.5x
~110k	5	18.513	5.8906	3x
~110k	20	418.165	17.3014	24x

The speed up becomes more prominent when more then just a few eigenvalues are searched, since that is a known limitation of the eigen_eigensystem solver - see #2286

Changelog

Added header only library Spectra to LinearSolversApplication
Added new solver for generalized sparse eigensystems: spectra_sym_g_eigs_shift
Removed -DEIGEN_DEFAULT_TO_ROW_MAJOR and explicitly specify ordering in Eigen types when interfacing with KRatos matrices.

…ctra

armingeiser · 2021-03-05T21:45:30Z

applications/LinearSolversApplication/CMakeLists.txt

@@ -16,7 +16,7 @@ message( "**** configuring KratosLinearSolversApplication ****" )
 ################### PYBIND11
 include(pybind11Tools)

-add_definitions( -DEIGEN_DEFAULT_TO_ROW_MAJOR -DEIGEN_MPL2_ONLY )
+# TODO add_definitions( -DEIGEN_MPL2_ONLY)


~~Spectra internally uses the Eigen::Cholesky, that however is not available with MPL2_ONLY. I will try to write an interface to another Eigen solver instead, ev, even MKL~~ Edit: Cholesky is not needed, the MPL2 flag can remain unchanged.

We have tried Spectra before without success: I think it was because it does not work with -DEIGEN_DEFAULT_TO_ROW_MAJOR. That does not have to be used however, as long as we specify the RowMajor ordering explicitly when interfacing Eigen matrices with Kratos matrices.

RiccardoRossi

dear @armingeiser i see myself obliged to block this PR since we are attempting to be restrictive with new libraries to be included.

Can u describe what is this library and why it is needed?

what is its license?

RiccardoRossi · 2021-03-05T22:10:16Z

just to add that adding a library to this Application is even more restricted since it is a default app used by everyone else

…with spectra" This reverts commit 69076e7.

…roblems with spectra"" This reverts commit f5c4472.

…olution

armingeiser · 2021-03-06T08:51:24Z

@RiccardoRossi i have updated the description, that should hopefully answer your questions.

armingeiser · 2021-03-06T08:53:15Z

@philbucher @emiroglu feel free to try this solver with your problems, would be interesting to see how it compares to the existing eigen_eigensystem solver. Note however that also this solver expects symmetric input matrices.

armingeiser · 2021-03-06T10:31:50Z

Does anyone have an idea why compilation works with gcc but not with clang or on windows? The error is always similar:

In file included from /__w/Kratos/Kratos/build/FullDebug/applications/LinearSolversApplication/cotire/KratosLinearSolversApplication_CXX_unity.cxx:4:
/__w/Kratos/Kratos/applications/LinearSolversApplication/custom_python/add_custom_solvers_to_python.cpp:148:43: error: no template named 'SpectraSymGEigsShiftSolver'
    using SpectraSymGEigsRealSolverType = SpectraSymGEigsShiftSolver<>;

philbucher · 2021-03-06T19:33:53Z

do you need the guards?

armingeiser · 2021-03-07T09:32:47Z

do you need the guards?

Yes.

Thanks for the hint. should be solved now.

I have not activted spectra for centos, windows, clang and nightly - since they also did not have MKL or Feast activated. Could be done however if desired.

oberbichler · 2021-03-07T11:06:06Z

@armingeiser do you remember the problems we had some years ago? Are they solved?

philbucher · 2021-03-29T12:20:27Z

applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_eigensolver.py

@@ -80,7 +80,7 @@ def _create_mechanical_solution_strategy(self):
        computing_model_part = self.GetComputingModelPart()

        solver_type = self.settings["eigensolver_settings"]["solver_type"].GetString()
-        if solver_type == "eigen_eigensystem":
+        if solver_type in ["eigen_eigensystem", "spectra_sym_g_eigs_shift"]: # TODO evaluate what has to be used for spectra


ping did you confirm this?

not yet, good catch.

philbucher

why did you call it spectra1?
spectra offers many more solvers right? Do you have plans to expose those too?

applications/LinearSolversApplication/tests/test_spectra_eigensystem_solver.py

philbucher · 2021-03-29T12:52:37Z

I just checked this with my shell turbine model (~415K dofs) and it is much faster that the current solver
for me when computing 10 eigenvalues it was more than 10 times faster, independent of number of threads

@KratosMultiphysics/technical-committee I recommend to approve this, this is a very valuable addition IMO

decision reached in @technicalcommittee

RiccardoRossi

approving on behalf of @KratosMultiphysics/technical-committee , on the assuption that it has the same license as Eigen

…tos matrices

armingeiser · 2021-04-14T15:52:59Z

@philbucher I have updated the PR with the things we discussed last time:

reuse the existing test
add a "define.h" where the Eigen matrix types are predefined using the row-major ordering. All other places in the app use these types now.

philbucher

nice work!

Please use the Kratos error macro to also get the stack trace when throwing, otherwise very nice!

philbucher · 2021-04-14T15:47:26Z

applications/LinearSolversApplication/linear_solvers_define.h

+template<typename _Scalar> using EigenDynamicMatrix = Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+template<typename _Scalar> using EigenDynamicVector = Eigen::Matrix<_Scalar, Eigen::Dynamic, 1>;
+
+template<typename _Scalar> using EigenSparseMatrix = Eigen::SparseMatrix<_Scalar, Eigen::RowMajor, int>;


not for this PR but why int and not size_t?

@oberbichler

Pardiso is expecting signed integers. Therefore the max size of a matrix is ‚only’ 2,147,483,647 which should be more then enough. Eigen is also using signed ints for dense matrices as unsigned ints are very error prone for algorithms.

that means that with shell elements we are limited to ~ 350 nodes (6dofs/node) which I agree is pretty big :)

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

philbucher · 2021-04-14T15:52:43Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

+        eigvecs = eigs.eigenvectors().transpose().colwise().reverse();
+
+        // --- normalization
+        // TODO seems to be normalized by Spectra already -> confirm!


it would be nice if you could take a look in a followup PR

philbucher · 2021-04-14T15:52:48Z

applications/LinearSolversApplication/custom_solvers/spectra_sym_g_eigs_shift_solver.h

+
+        OpType op(a, b);
+        BOpType Bop(b);
+        const int ncv = 3 * nroot;  // TODO find a good value


philbucher · 2021-04-14T15:54:39Z

* add a "define.h" where the Eigen matrix types are predefined using the row-major ordering.  All other places in the app use these types now.

this is very cool!

…ym_g_eigs_shift_solver.h Co-authored-by: Philipp Bucher <[email protected]>

armingeiser · 2021-04-14T16:03:51Z

* why did you call it spectra**1**?

Because this is the current master branch, which will soon become spectra version 1.0 (yixuan/spectra#113)

* spectra offers many more solvers right? Do you have plans to expose those too?

It has many more features, however, currently i don't have any plans to further extend.
But i think with this first version, it will be easy to extend in case anyone needs further features from spectra.

philbucher · 2021-04-14T20:16:19Z

applications/LinearSolversApplication/linear_solvers_define.h

+namespace Kratos
+{
+
+// IMPORTANT:


It is the other way round no?
Eigen default is column major while Kratos is row major

🙈 fixed the comment

applications/LinearSolversApplication/linear_solvers_define.h

RiccardoRossi · 2021-04-15T19:09:09Z

guys, our matrices are row_major... Riccardo

…

On Thu, Apr 15, 2021, 9:08 AM Philipp Bucher ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In applications/LinearSolversApplication/linear_solvers_define.h <#8440 (comment)> : > +{ + +// IMPORTANT: +// +// The matrices in Kratos use colum-major storage ordering. +// Eigen matrices, however, use row-major storage ordering by default. +// +// If these two worlds are combined, it is crucial that Eigen matrices are used +// with **column-major** storage order to match the Kratos matrices! +// +// To simplify things, the most commonly used types are defined below: + +template<typename _Scalar> using EigenDynamicMatrix = Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; +template<typename _Scalar> using EigenDynamicVector = Eigen::Matrix<_Scalar, Eigen::Dynamic, 1>; + +template<typename _Scalar> using EigenSparseMatrix = Eigen::SparseMatrix<_Scalar, Eigen::RowMajor, int>; that means that with shell elements we are limited to ~ 350 nodes (6dofs/node) which I agree is pretty big :) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#8440 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB5PWEKFMJ56NMYTIKHZC7LTI2GILANCNFSM4YV55D4A> .

philbucher · 2021-04-16T06:50:59Z

Seems like you didn’t check the latest changes

can you check master?

armingeiser added 5 commits March 5, 2021 21:56

remove -DEIGEN_DEFAULT_TO_ROW_MAJOR since it causes problems with spe…

69076e7

…ctra

add spectra lib master from 2021.03.05

46e994d

add interface to spectra SymGEigsShiftSolver and test

3bb9659

fix usage in structmech

b2b3350

activate spectra in CI

5e9eef7

armingeiser added the External library label Mar 5, 2021

armingeiser commented Mar 5, 2021

View reviewed changes

RiccardoRossi previously requested changes Mar 5, 2021

View reviewed changes

armingeiser added 4 commits March 6, 2021 08:22

Revert "remove -DEIGEN_DEFAULT_TO_ROW_MAJOR since it causes problems …

f5c4472

…with spectra" This reverts commit 69076e7.

Revert "Revert "remove -DEIGEN_DEFAULT_TO_ROW_MAJOR since it causes p…

33f912b

…roblems with spectra"" This reverts commit f5c4472.

add MPL2_ONLY flag and enable pardiso solver in spectra eigensystem s…

8052b4b

…olution

add spectra test to testsuite

e2585c5

armingeiser added the Discussion label Mar 6, 2021

spectra can be used without mkl

846cbbd

armingeiser changed the title ~~Linear solvers/spectra eigenvalue~~ [LinearSolvers] new solver for (symmetric) generalized eigenproblems using the spectra library Mar 6, 2021

cleaning

c28de0c

missing guards

09b1e06

skip test if spectra is not used

eb9339e

armingeiser marked this pull request as ready for review March 7, 2021 09:51

armingeiser requested review from a team as code owners March 7, 2021 09:51

philbucher reviewed Mar 29, 2021

View reviewed changes

applications/LinearSolversApplication/tests/test_spectra_eigensystem_solver.py Outdated Show resolved Hide resolved

applications/LinearSolversApplication/tests/test_spectra_eigensystem_solver.py Outdated Show resolved Hide resolved

RiccardoRossi previously approved these changes Apr 14, 2021

View reviewed changes

armingeiser added 2 commits April 14, 2021 16:19

Merge branch 'master' into LinearSolvers/spectra-eigenvalue

f37a5dd

explicitly define storage order for eigen matrices that interface kra…

7a8792e

…tos matrices

armingeiser dismissed RiccardoRossi’s stale review via 7a8792e April 14, 2021 14:40

armingeiser added 2 commits April 14, 2021 17:42

create types for easy interface of kratos matrices with eigen matrices

0d34a26

reuse exsiting test

6c56e3f

armingeiser requested a review from philbucher April 14, 2021 15:51

philbucher reviewed Apr 14, 2021

View reviewed changes

armingeiser and others added 2 commits April 14, 2021 17:57

Update applications/LinearSolversApplication/custom_solvers/spectra_s…

8000e49

…ym_g_eigs_shift_solver.h Co-authored-by: Philipp Bucher <[email protected]>

Update applications/LinearSolversApplication/custom_solvers/spectra_s…

9956676

…ym_g_eigs_shift_solver.h Co-authored-by: Philipp Bucher <[email protected]>

philbucher previously approved these changes Apr 14, 2021

View reviewed changes

philbucher reviewed Apr 14, 2021

View reviewed changes

armingeiser commented Apr 15, 2021

View reviewed changes

applications/LinearSolversApplication/linear_solvers_define.h Outdated Show resolved Hide resolved

Update applications/LinearSolversApplication/linear_solvers_define.h

89f0a2b

armingeiser dismissed philbucher’s stale review via 89f0a2b April 15, 2021 06:52

philbucher approved these changes Apr 15, 2021

View reviewed changes

armingeiser merged commit 261a3d5 into master Apr 15, 2021

armingeiser deleted the LinearSolvers/spectra-eigenvalue branch April 15, 2021 10:43

philbucher mentioned this pull request Jun 23, 2021

[Structral][EigenSolver] using sprecta by default #8907

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[LinearSolvers] new solver for (symmetric) generalized eigenproblems using the spectra library #8440

[LinearSolvers] new solver for (symmetric) generalized eigenproblems using the spectra library #8440

armingeiser commented Mar 5, 2021 •

edited

Loading

armingeiser Mar 5, 2021 •

edited

Loading

RiccardoRossi left a comment

RiccardoRossi commented Mar 5, 2021

armingeiser commented Mar 6, 2021

armingeiser commented Mar 6, 2021

armingeiser commented Mar 6, 2021

philbucher commented Mar 6, 2021

armingeiser commented Mar 7, 2021 •

edited

Loading

oberbichler commented Mar 7, 2021

philbucher Mar 29, 2021

armingeiser Mar 29, 2021

philbucher left a comment

philbucher commented Mar 29, 2021

RiccardoRossi left a comment

armingeiser commented Apr 14, 2021

philbucher left a comment

philbucher Apr 14, 2021

oberbichler Apr 15, 2021

philbucher Apr 15, 2021

philbucher Apr 14, 2021

philbucher Apr 14, 2021

philbucher commented Apr 14, 2021

armingeiser commented Apr 14, 2021

philbucher Apr 14, 2021

armingeiser Apr 15, 2021

RiccardoRossi commented Apr 15, 2021 via email

philbucher commented Apr 16, 2021

[LinearSolvers] new solver for (symmetric) generalized eigenproblems using the spectra library #8440

[LinearSolvers] new solver for (symmetric) generalized eigenproblems using the spectra library #8440

Conversation

armingeiser commented Mar 5, 2021 • edited Loading

armingeiser Mar 5, 2021 • edited Loading

Choose a reason for hiding this comment

RiccardoRossi left a comment

Choose a reason for hiding this comment

RiccardoRossi commented Mar 5, 2021

armingeiser commented Mar 6, 2021

armingeiser commented Mar 6, 2021

armingeiser commented Mar 6, 2021

philbucher commented Mar 6, 2021

armingeiser commented Mar 7, 2021 • edited Loading

oberbichler commented Mar 7, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

philbucher left a comment

Choose a reason for hiding this comment

philbucher commented Mar 29, 2021

RiccardoRossi left a comment

Choose a reason for hiding this comment

armingeiser commented Apr 14, 2021

philbucher left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

philbucher commented Apr 14, 2021

armingeiser commented Apr 14, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

RiccardoRossi commented Apr 15, 2021 via email

philbucher commented Apr 16, 2021

armingeiser commented Mar 5, 2021 •

edited

Loading

armingeiser Mar 5, 2021 •

edited

Loading

armingeiser commented Mar 7, 2021 •

edited

Loading