compute Eigenvalues and -vectors in 1-3d without LAPACK

When dealing with small matrices calling LAPACK to calculate Eigenvalues and Eigenvectors can be quite slow. 1d- and 2d-specializations are quite straightforward, for the 3d-specialization I´m using the approach presented in https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf.

added feature label

@christi The main work should be done here, however, there are still two design choices to be made: I implemented the new Eigenvector functions with the signature

eigenValuesVectors(const FieldMatrix<K, dim, dim>& matrix,
                  FieldVector<K, dim>& eigenValues,
                  FieldVector<FieldVector<K, dim>, dim>& eigenVectors)

while the existing Lapack-binding has the signature

eigenValuesVectorsLapack(const FieldMatrix<K, dim, dim>& matrix,
                        FieldVector<K, dim>& eigenValues,
                        FieldMatrix<K, dim, dim>& eigenVectors,
                        const char& jobz)

First, I propose to adopt the jobz-tag (which controls whether to only calculate eigenvalues or eigenvalues and eigenvectors) and get rid of the eigenValues-function. Second, I feel that FieldVector<FieldVector<K, dim>, dim> is the more natural type for a set of eigenvectors, so I´d like to change the return type in the Lapack-binding. Both of those changes would most likely lead to subsequent changes in dependend code. How do you feel about this?

Second, I feel that FieldVector<FieldVector<K, dim>, dim> is the more natural type for a set of eigenvectors, so I´d like to change the return type in the Lapack-binding. Both of those changes would most likely lead to subsequent changes in dependend code. How do you feel about this?

I'd stick to the current interface. It is not uncommon that you want to use the eigenvector basis to map an operator or a vector. In this case FieldMatrix is the right type. If you want to access individual eigenevectors, you can still do so in the current interface, as FieldMatrix::operator[] returns a FieldVector.

added 1 commit

• c35c18f4 - resolve signature difference with lapack-binding

Compare with previous version

Okay, I wasnt aware of the operator[] functionality - then there is no need to change anything. What about my other question?

What about my other question?

I'm undecided and have to think about it.

An additional parameter sounds good, but I don't like the char argument

Just a small side note: Please don't use FieldVector<FieldVector<...>> at all. Since FieldVector<...> is intended to contain entries from a field we dropped all occurrences long ago. If this should represent a linear operator or tensor, use FieldMatrix<...>, it it's just a statically sized collection of vectors then use std::array<FieldVector<...>>.

added 2 commits

• f1a47aa5 - [temp] reenable tests for dim>3
• 08fda1d6 - restructure interface and update tests

Compare with previous version

There are currently two issues with this MR:

1. if constexpr only works for newer compilers, so you have to switch to the updated .gitlab-ci.yml (easiest is to rebase on current master)
2. you removed the eigenValues functions. We still need them for backwards compatibility and can think about deprecating them...

mentioned in merge request !724 (closed)

added 23 commits

• 08fda1d6...cd4c780f - 15 commits from branch master
• b3b21f87 - [test] test multiplicity of 3x3 EV calculation
• c7786333 - add 1d and 2d functionality
• 5a3526e1 - add 3d specialization
• a8f44b95 - change interface, minor bugfixes
• 91d8d961 - update tests
• ebb8e965 - resolve signature difference with lapack-binding
• 6dd1a09a - [temp] reenable tests for dim>3
• 9a4fb1f1 - restructure interface and update tests

Compare with previous version

added 2 commits

• b958a4c0 - re-add eigenValues for compatibility
• f0cbbe33 - deprecate 'eigenValues()'

Compare with previous version

adressing 2.: Yes, I'm still thinking about how to solve this best. I'm also not happy with the ComputationJob-tag yet - in practice this can get really long (Dune::FMatrixHelp::eigenValuesVectors<Dune::FMatrixHelp::ComputationJob::EigenvaluesEigenvectors>(...)).

How about the following...

• we keep the two interfaces fr eigenvalues and eigenvalues+vectors
• these forward then to the same implementation using a particular tag

Thats probably the easiest solution. I was also thinking about keeping eigenValues and somehow making ComputationJob::EigenvaluesEigenvectors the "default-tag" for eigenValuesVectors- passing a tag called EigenvaluesEigenvectors to a function eigenValuesVectors (which is expected to compute exactly that) is also kind of counterintuitive...

added 1 commit

• d4d36dc5 - simplify interface: hide ComputationTag internally

Compare with previous version

I suggest to move this array to eigenValuesVectorsLapackImpl, as it is only used there... actually, you wouldn't even need an array for this.

as you want to reuse the enum as an index, you should explicitly assign values.

you can simply write

  const char jobz = "nv"[Tag];

Okay. I added some documentation so the usage is clear if this will ever be needed again in the future.

mentioned in merge request !727 (closed)