Made the cholesky decomposition code use the xPOTRF routines in LAPACK

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

dlib/matrix/matrix_cholesky.h

@@ -10,6 +10,10 @@
 #include "matrix_subexp.h"
 #include <cmath>
 
+#ifdef DLIB_USE_LAPACK
+#include "lapack/potrf.h"
+#endif
+
 namespace dlib 
 {
 
@@ -105,6 +109,32 @@ namespace dlib
             << "\n\tthis:      " << this
             );
 
+#ifdef DLIB_USE_LAPACK
+        L_ = A_;
+        const type eps = max(abs(diag(L_)))*std::sqrt(std::numeric_limits<type>::epsilon())/100;
+
+        // check if the matrix is actually symmetric
+        bool is_symmetric = true;
+        for (long r = 0; r < L_.nr() && is_symmetric; ++r)
+        {
+            for (long c = r+1; c < L_.nc() && is_symmetric; ++c)
+            {
+                // this is approximately doing: is_symmetric = is_symmetric && ( L_(k,j) == L_(j,k))
+                is_symmetric = is_symmetric && (std::abs(L_(r,c) - L_(c,r)) < eps ); 
+            }
+        }
+
+        // now compute the actual cholesky decomposition
+        int info = lapack::potrf('L', L_);
+
+        // check if its really SPD
+        if (info == 0 && is_symmetric && min(abs(diag(L_))) > eps*100)
+            isspd = true;
+        else
+            isspd = false;
+
+        L_ = lowerm(L_);
+#else
         const_temp_matrix<EXP> A(A_);
 
 
@@ -153,6 +183,7 @@ namespace dlib
                 L_(j,k) = 0.0;
             }
         }
+#endif
     }
 
 // ----------------------------------------------------------------------------------------

dlib/matrix/matrix_la.h

@@ -13,6 +13,10 @@
 #include "matrix_cholesky.h"
 #include "matrix_eigenvalue.h"
 
+#ifdef DLIB_USE_LAPACK
+#include "lapack/potrf.h"
+#endif
+
 namespace dlib
 {
 
@@ -1223,8 +1227,14 @@ convergence:
             << "\n\tA.nr(): " << A.nr()
             << "\n\tA.nc(): " << A.nc() 
             );
-
         typename matrix_exp<EXP>::matrix_type L(A.nr(),A.nc());
+
+#ifdef DLIB_USE_LAPACK
+        L = A;
+        lapack::potrf('L', L);
+        // mask out upper triangular area
+        return lowerm(L);
+#else
         typedef typename EXP::type T;
         set_all_elements(L,0);
 
@@ -1273,6 +1283,8 @@ convergence:
         }
 
         return L;
+#endif
+
     }
 
 // ----------------------------------------------------------------------------------------

dlib/matrix/matrix_la_abstract.h

@@ -181,6 +181,9 @@ namespace dlib
             - else
                 - returns a matrix with the same dimensions as A but it 
                   will have a bogus value.  I.e. it won't be a decomposition.
+
+            - If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF 
+              is used to compute the cholesky decomposition.
     !*/
 
 // ----------------------------------------------------------------------------------------
@@ -387,6 +390,9 @@ namespace dlib
                 If the matrix is not symmetric or positive definite, the function
                 computes only a partial decomposition.  This can be tested with
                 the is_spd() flag.
+            
+                If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF 
+                is used to compute the cholesky decomposition.
         !*/
 
     public:

dlib/test/matrix_chol.cpp

@@ -96,7 +96,7 @@ namespace
         DLIB_TEST_MSG(equal(m*test.solve(col), col,eps),max(abs(m*test.solve(m2)- m2)));
 
         // now make us a non-spd matrix
-        if (m.nr() > 1)
+        if (m.nr() > 2)
         {
             matrix<type> sm(lowerm(m));
             sm(1,1) = 0;