Clean up fmatrixev
This incorporates !277 (merged) (because bug fixes are good) but not !278 (closed) (because we might want to rethink an interface before we extend it).
So far, this addresses at least one bug in addition to !277 (merged) but leaves the old interface intact. This might also be a good opportunity to give fmatrixev and dynmatrixev a unified interface.
Merge request reports
Activity
mentioned in merge request !278 (closed)
A question worth addressing: Did
eigenValuesNonSymever work? Because it appears to assume that as template parameterCone passes in a class for instancescof which one can doc.real = ...; c.imag = ...;because that's what's done here:
eigenValues[i].real = eigenR[i]; eigenValues[i].imag = eigenI[i];This, however, is not possible for
std::complex<T>! Yes, withstd::complex<double> cyou can access the real and imaginary part by calling(!)c.real()andc.imag(), but to set them you need to do something likec = { ..., ... }; so is there a class that people are using that it would work for? Or was the code just never tested?Edited by Elias Pipping-
As mentioned before, one advantage of the native implementation over lapack is the ability to work with more general fields; I tried the
GMPFieldclass and hit a few obstacles, though:-
GMP does not come with trigonometric functions, in particular not
cosoracos(MPFR provides them, albeit not through such a plain name; maybe through boost's gmp/mpfr wrapping one would end up with something workable) -
Consequently,
Dune::MathematicalConstants::piis not and cannot be defined forDune::GMPField -
And finally the lines
real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() ); real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() );in densematrix.hh do not compile for a matrix with
GMPFieldentries because the first argument is aGMPFieldin each case but the latter is an expression which one'd have to pass to astatic_castfirst (in addition, one's need to switch tousing std::maxand use plainmax()to allow for ADL). Not sure anyone cares about that, though, because it's unclear if one'd actually get reasonable accuracy for aGMPFieldthrough the native implementation.
Edited by Elias Pipping -
Hi Elias, you may want to consult b9e712ac on the question of whether the Lapack eigenvalue stuff ever worked.
-
@oliver.sander That's both shocking and eye-opening, thank you!
-
This MR is still incomplete in the sense that a new interface is missing; a new interface could
- eliminate the ambiguity of the current interface (whether the native c++ dune implementation or lapack is used should be controllable)
- expose the new functionality of finding eigenvalues of matrices with complex entries (be the matrices hermitian or not)
That said, the functionality is there, it's tested and it works. I prefer to act as a worker rather than a designer so I consider my share of the work here done for now, hence I've removed the WIP label that might keep others from commenting.
Edited by Elias Pipping -
Note to self: Judging by 7cb10fb1 it will be necessary to move the fortran routines back into the library.
I haven't checked in detail what you do in this MR, but I like the general idea. Please ping me, if you feel of merging it. If nobody comments, wait a couple of days and I am willing to merge.
Edited by Christoph Grüninger-
@gruenich: thank you for commenting! My current impression is that in light of 7cb10fb1, I might actually be introducing a regression with this MR in its current form. I should try and fix that; at the same time i would greatly appreciate it if somebody else could take over this subtask — not because I expect it to be particularly hard or time consuming but because my health makes it impossible to predict if I’ll get around to it myself (be that in the near future or ever).
As for what this MR currently does: it mostly removes a hardcoded tolerance in one place, and instead of real double matrices you can now have real or complex float or double matrices. Also, fmatrixev was previously entirely untested if I’m not mistaken. Now there are quite a few tests.
Edited by Elias Pipping -
(Sorry for the noise - can someone more knowledgable on the topic chime in maybe?) - the Fortran function declarations might be safe already because they’re in a namespace; @robert.kloefkorn can you comment on what type of issues you were seeing prior to 7cb10fb1 and whether you would expect them to resurface the way this MR is now?
-
@christi Is this of help or connected to your recent work with eigen values?
