[!711] eigenValuesNonSym: optionally also compute eigenvectors

Merge branch 'eigenValuesNonSym-compute-eigenvectors' into 'master'

See merge request [core/dune-common!711]

  [core/dune-common!711]: Nonecore/dune-common/merge_requests/711

dune/common/dynmatrixev.hh

@@ -3,6 +3,7 @@
 #ifndef DUNE_DYNMATRIXEIGENVALUES_HH
 #define DUNE_DYNMATRIXEIGENVALUES_HH
 
+#include <algorithm>
 #include <memory>
 
 #include "dynmatrix.hh"
@@ -26,17 +27,20 @@ namespace Dune {
         \param[in]  matrix matrix eigenvalues are calculated for
         \param[out] eigenValues FieldVector that contains eigenvalues in
                     ascending order
+        \param[out] eigenVectors (optional) list of right eigenvectors
 
         \note LAPACK::dgeev is used to calculate the eigen values
      */
     template <typename K, class C>
     static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
-                                  DynamicVector<C>& eigenValues)
+                                  DynamicVector<C>& eigenValues,
+                                  std::vector<DynamicVector<K>>* eigenVectors = nullptr
+      )
     {
       {
         const long int N = matrix.rows();
         const char jobvl = 'n';
-        const char jobvr = 'n';
+        const char jobvr = eigenVectors ? 'v' : 'n';
 
 
         // matrix to put into dgeev
@@ -55,11 +59,13 @@ namespace Dune {
         // working memory
         auto eigenR = std::make_unique<double[]>(N);
         auto eigenI = std::make_unique<double[]>(N);
-        auto work = std::make_unique<double[]>(3*N);
+
+        const long int lwork = eigenVectors ? 4*N : 3*N;
+        auto work = std::make_unique<double[]>(lwork);
+        auto vr = eigenVectors ? std::make_unique<double[]>(N*N) : std::unique_ptr<double[]>{};
 
         // return value information
         long int info = 0;
-        const long int lwork = 3*N;
 
         // call LAPACK routine (see fmatrixev_ext.cc)
         eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, matrixVector.get(), &N,
@@ -75,6 +81,15 @@ namespace Dune {
         eigenValues.resize(N);
         for (int i=0; i<N; ++i)
           eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
+
+        if (eigenVectors) {
+          eigenVectors->resize(N);
+          for (int i = 0; i < N; ++i) {
+            auto& v = (*eigenVectors)[i];
+            v.resize(N);
+            std::copy(vr.get() + N*i, vr.get() + N*(i+1), &v[0]);
+          }
+        }
       }
     }