Setup the eigenvalue_decomposition to use LAPACK

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%403822

dlib/matrix/matrix_eigenvalue.h

@@ -12,6 +12,12 @@
 #include <complex>
 #include <cmath>
 
+#ifdef DLIB_USE_LAPACK
+#include "lapack/geev.h"
+#include "lapack/syev.h"
+#include "lapack/syevr.h"
+#endif
+
 namespace dlib 
 {
 
@@ -171,16 +177,32 @@ namespace dlib
         {
             V = A;
 
+#ifdef DLIB_USE_LAPACK
+            matrix<type,0,1,mem_manager_type, layout_type> work;
+            e = 0;
+            // I would use syevr but the last time I checked there was a bug in the 
+            // Intel MKL's implementation of syevr.
+            lapack::syev('V', 'L', V,  d, work);
+#else
             // Tridiagonalize.
             tred2();
 
             // Diagonalize.
             tql2();
+#endif
 
         } 
         else 
         {
+
+#ifdef DLIB_USE_LAPACK
+            matrix<type,0,0,mem_manager_type, column_major_layout> temp, vl, vr, work;
+            temp = A;
+            lapack::geev('N', 'V', temp, d, e, vl, vr, work);
+            V = vr;
+#else
             H = A;
+
             ort.set_size(n);
 
             // Reduce to Hessenberg form.
@@ -188,6 +210,7 @@ namespace dlib
 
             // Reduce Hessenberg to real Schur form.
             hqr2();
+#endif
         }
     }
 
@@ -222,11 +245,19 @@ namespace dlib
 
         V = A;
 
+#ifdef DLIB_USE_LAPACK
+        matrix<type,0,1,mem_manager_type, layout_type> work;
+        e = 0;
+        // I would use syevr but the last time I checked there was a bug in the 
+        // Intel MKL's implementation of syevr.
+        lapack::syev('V', 'L', V,  d, work);
+#else
         // Tridiagonalize.
         tred2();
 
         // Diagonalize.
         tql2();
+#endif
 
     }
 

dlib/matrix/matrix_la_abstract.h

@@ -632,6 +632,8 @@ namespace dlib
                     - V*trans(V) should be equal to the identity matrix.  That is, all the
                       eigenvectors in V should be orthonormal. 
                         - So A == V*D*trans(V)
+                    - If DLIB_USE_LAPACK is #defined then this object uses the xSYEV LAPACK
+                      routine.
 
                 On the other hand, if A is not symmetric then:
                     - Some of the eigenvalues and eigenvectors might be complex numbers.  
@@ -642,6 +644,8 @@ namespace dlib
                     - V*trans(V) won't be equal to the identity matrix but it is usually
                       invertible.  So A == V*D*inv(V) is usually a valid statement but
                       A == V*D*trans(V) won't be.
+                    - If DLIB_USE_LAPACK is #defined then this object uses the xGEEV LAPACK
+                      routine.
         !*/
 
     public:

dlib/test/matrix_eig.cpp

@@ -26,11 +26,6 @@ namespace
 
     dlib::rand::float_1a rnd;
 
-// ----------------------------------------------------------------------------------------
-
-    template <typename mat_type>
-    const matrix<typename mat_type::type> symm(const mat_type& m) { return m*trans(m); }
-
 // ----------------------------------------------------------------------------------------
 
     template <typename type>
@@ -65,8 +60,8 @@ namespace
 
 // ----------------------------------------------------------------------------------------
 
-    template <typename matrix_type>
-    void test_eigenvalue ( const matrix_type& m)
+    template <typename matrix_type, typename U>
+    void test_eigenvalue_impl ( const matrix_type& m,  const eigenvalue_decomposition<U>& test )
     {
         typedef typename matrix_type::type type;
         const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon());
@@ -74,8 +69,6 @@ namespace
         print_spinner();
 
 
-        eigenvalue_decomposition<matrix_type> test(m);
-
         DLIB_TEST(test.dim() == m.nr());
 
         // make sure all the various ways of asking for the eigenvalues are actually returning a
@@ -128,6 +121,18 @@ namespace
         }
     }
 
+// ----------------------------------------------------------------------------------------
+
+    template <typename matrix_type>
+    void test_eigenvalue ( const matrix_type& m )
+    {
+        eigenvalue_decomposition<matrix_type> test(m);
+        test_eigenvalue_impl(m, test);
+
+        eigenvalue_decomposition<matrix_type> test_symm(make_symmetric(m));
+        test_eigenvalue_impl(make_symmetric(m), test_symm);
+    }
+
 // ----------------------------------------------------------------------------------------
 
     void matrix_test_double()
@@ -143,14 +148,6 @@ namespace
         test_eigenvalue(10*randm<double,1,1>());
         test_eigenvalue(10*randm<double,2,2>());
         test_eigenvalue(10*randm<double,3,3>());
-
-        test_eigenvalue(10*symm(randm<double>(1,1)));
-        test_eigenvalue(10*symm(randm<double>(2,2)));
-        test_eigenvalue(10*symm(randm<double>(3,3)));
-        test_eigenvalue(10*symm(randm<double>(4,4)));
-        test_eigenvalue(10*symm(randm<double>(15,15)));
-        test_eigenvalue(10*symm(randm<double>(150,150)));
-
     }
 
 // ----------------------------------------------------------------------------------------
@@ -168,13 +165,6 @@ namespace
         test_eigenvalue(10*randm<float,1,1>());
         test_eigenvalue(10*randm<float,2,2>());
         test_eigenvalue(10*randm<float,3,3>());
-
-        test_eigenvalue(10*symm(randm<float>(1,1)));
-        test_eigenvalue(10*symm(randm<float>(2,2)));
-        test_eigenvalue(10*symm(randm<float>(3,3)));
-        test_eigenvalue(10*symm(randm<float>(4,4)));
-        test_eigenvalue(10*symm(randm<float>(15,15)));
-        test_eigenvalue(10*symm(randm<float>(50,50)));
     }
 
 // ----------------------------------------------------------------------------------------