diff --git a/cholesky.h b/cholesky.h
index 5830c71..4fda38d 100644
--- a/cholesky.h
+++ b/cholesky.h
@@ -21,8 +21,20 @@ Authors:
     Caner Candan <caner.candan@thalesgroup.com>
 */
 
+#include <stdexcept>
+
+#include <boost/numeric/ublas/symmetric.hpp>
+using namespace boost::numeric;
+
 namespace cholesky {
 
+class NotDefinitePositive : public std::runtime_error
+{
+public:
+    NotDefinitePositive( const std::string & what ) : std::runtime_error( what ) {}
+};
+
+
 /** Cholesky decomposition, given a matrix V, return a matrix L
  * such as V = L L^T (L^T being the transposed of L).
  *
@@ -31,7 +43,7 @@ namespace cholesky {
  * Thus, expect a (lower) triangular matrix.
  */
 template< typename T >
-class CholeskyBase
+class Cholesky
 {
 public:
     //! The covariance-matrix is symetric
@@ -58,12 +70,12 @@ public:
     }
 
     /** Compute the factorization and cache the result */
-    virtual void operator()( const CovarMat& V ) = 0;
+    virtual void factorize( const CovarMat& V ) = 0;
 
     /** Compute the factorization and return the result */
-    virtual const FactorMat& operator()( const CovarMat& V ) 
+    const FactorMat& operator()( const CovarMat& V )
     {
-        (*this)( V );
+        this->factorize(V);
         return decomposition();
     }
 
@@ -123,19 +135,19 @@ protected:
  * will raise an assert before calling the square root on a negative number.
  */
 template< typename T >
-class CholeskyLLT : public CholeskyBase<T>
+class LLT : public Cholesky<T>
 {
 public:
-    virtual void operator()( const CovarMat& V )
+    virtual void factorize( const typename Cholesky<T>::CovarMat& V )
     {
         unsigned int N = assert_properties( V );
 
         unsigned int i=0, j=0, k;
-        _L(0, 0) = sqrt( V(0, 0) );
+        this->_L(0, 0) = sqrt( V(0, 0) );
 
         // end of the column
         for ( j = 1; j < N; ++j ) {
-            _L(j, 0) = V(0, j) / _L(0, 0);
+            this->_L(j, 0) = V(0, j) / this->_L(0, 0);
         }
 
         // end of the matrix
@@ -143,32 +155,36 @@ public:
             // diagonal
             double sum = 0.0;
             for ( k = 0; k < i; ++k) {
-                sum += _L(i, k) * _L(i, k);
+                sum += this->_L(i, k) * this->_L(i, k);
             }
 
-            _L(i,i) = L_i_i( V, i, sum );
+            this->_L(i,i) = L_i_i( V, i, sum );
 
             for ( j = i + 1; j < N; ++j ) { // rows
                 // one element
                 sum = 0.0;
                 for ( k = 0; k < i; ++k ) {
-                    sum += _L(j, k) * _L(i, k);
+                    sum += this->_L(j, k) * this->_L(i, k);
                 }
 
-                _L(j, i) = (V(j, i) - sum) / _L(i, i);
+                this->_L(j, i) = (V(j, i) - sum) / this->_L(i, i);
 
             } // for j in ]i,N[
         } // for i in [1,N[
     }
 
-
     /** The step of the standard LLT algorithm where round off errors may appear */
-    inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const
+    inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
     {
         // round-off errors may appear here
-        assert( V(i,i) - sum >= 0 );
+        if( V(i,i) - sum < 0 ) {
+            std::ostringstream oss;
+            oss << "V(" << i << "/" << V.size1() << ")=" << V(i,i) << " - sum=" << sum << "\t== " << V(i,i)-sum << " < 0 ";
+            throw NotDefinitePositive(oss.str());
+        }
         return sqrt( V(i,i) - sum );
     }
+
 };
 
 
@@ -181,10 +197,10 @@ public:
  * hack to avoid crash.
  */
 template< typename T >
-class CholeskyLLTabs : public CholeskyLLT<T>
+class LLTabs : public LLT<T>
 {
 public:
-    inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const
+    inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
     {
         /***** ugly hack *****/
         return sqrt( fabs( V(i,i) - sum) );
@@ -200,10 +216,10 @@ public:
  * hack to avoid crash.
  */
 template< typename T >
-class CholeskyLLTzero : public CholeskyLLT<T>
+class LLTzero : public LLT<T>
 {
 public:
-    inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const
+    inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
     {
         T Lii;
         if(  V(i,i) - sum >= 0 ) {
@@ -224,18 +240,18 @@ public:
  * The factorized matrix is (L D^1/2), because V = (L D^1/2) (L D^1/2)^T
  */
 template< typename T >
-class CholeskyLDLT : public CholeskyBase<T>
+class LDLT : public Cholesky<T>
 {
 public:
-    virtual void operator()( const CovarMat& V )
+    virtual void factorize( const typename Cholesky<T>::CovarMat& V )
     {
         // use "int" everywhere, because of the "j-1" operation
         int N = assert_properties( V );
         // example of an invertible matrix whose decomposition is undefined
         assert( V(0,0) != 0 ); 
 
-        FactorMat L = ublas::zero_matrix<T>(N,N);
-        FactorMat D = ublas::zero_matrix<T>(N,N);
+        typename Cholesky<T>::FactorMat L = ublas::zero_matrix<T>(N,N);
+        typename Cholesky<T>::FactorMat D = ublas::zero_matrix<T>(N,N);
         D(0,0) = V(0,0);
 
         for( int j=0; j<N; ++j ) { // each columns
@@ -257,19 +273,19 @@ public:
             } // for i in rows
         } // for j in columns
         
-        _L = root( L, D );
+        this->_L = root( L, D );
     }
 
 
-    inline FactorMat root( FactorMat& L, FactorMat& D )
+    inline typename Cholesky<T>::FactorMat root( typename Cholesky<T>::FactorMat& L, typename Cholesky<T>::FactorMat& D )
     {
-        // now compute the final symetric matrix: _L = L D^1/2
+        // now compute the final symetric matrix: this->_L = L D^1/2
         // remember that V = ( L D^1/2) ( L D^1/2)^T
 
         // fortunately, the square root of a diagonal matrix is the square 
         // root of all its elements
-        FactorMat sqrt_D = D;
-        for(int i=0; i<N; ++i) {
+        typename Cholesky<T>::FactorMat sqrt_D = D;
+        for(unsigned int i=0; i<D.size1(); ++i) {
             sqrt_D(i,i) = sqrt(D(i,i));
         }
 
diff --git a/test.cpp b/test.cpp
index 70fbd2a..37de315 100644
--- a/test.cpp
+++ b/test.cpp
@@ -1,3 +1,29 @@
+/*
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this library; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+
+Copyright (C) 2010 Thales group
+*/
+/*
+Authors:
+    Johann Dréo <johann.dreo@thalesgroup.com>
+*/
+
+#include <iostream>
+#include <limits>
+#include <iomanip>
+#include <numeric>
 
 #include "cholesky.h"
 
@@ -79,7 +105,7 @@ MT error( const MT& M1, const MT& M2 )
 template<typename MT >
 double trigsum( const MT& M )
 {
-    double sum;
+    double sum = 0;
     for( unsigned int i=0; i<M.size1(); ++i ) {
         for( unsigned int j=i; j<M.size2(); ++j ) { // triangular browsing
             sum += fabs( M(i,j) ); // absolute deviation
@@ -96,52 +122,77 @@ double sum( const T& c )
 }
 
 
-typedef double real;
 
 
 int main(int argc, char** argv)
 {
     srand(time(0));
 
-    unsigned int N = 4; // size of matrix
-    unsigned int R = 1000; // nb of repetitions
+    unsigned int M = 10; // sample size
+    unsigned int N = 12; // variables number
+    unsigned int F = 10; // range factor
+    unsigned int R = 1; // nb of repetitions
 
     if( argc >= 2 ) {
-        N = std::atoi(argv[1]);
+        M = std::atoi(argv[1]);
     }
     if( argc >= 3 ) {
-        R = std::atoi(argv[2]);
+        N = std::atoi(argv[2]);
+    }
+    if( argc >= 4 ) {
+        F = std::atoi(argv[3]);
+    }
+    if( argc >= 5 ) {
+        R = std::atoi(argv[4]);
     }
 
-    std::cout << "Usage: t-cholesky [matrix size] [repetitions]" << std::endl;
-    std::cout << "matrix size = " << N << std::endl;
-    std::cout << "repetitions = " << R << std::endl;
+    std::clog << "Usage: test [sample size] [variables number] [random range] [repetitions]" << std::endl;
+    std::clog << "\tsample size  = " << M << std::endl;
+    std::clog << "\tmatrix size  = " << N << std::endl;
+    std::clog << "\trandom range = " << F << std::endl;
+    std::clog << "\trepetitions  = " << R << std::endl;
 
-    typedef cholesky::Cholesky::CovarMat CovarMat;
-    typedef cholesky::Cholesky::FactorMat FactorMat;
+    typedef double real;
+    typedef cholesky::Cholesky<real>::CovarMat CovarMat;
+    typedef cholesky::Cholesky<real>::FactorMat FactorMat;
 
-    cholesky::LLT     llt();
-    cholesky::LLTabs  llta();
-    cholesky::LLTzero lltz();
-    cholesky::LDLT    ldlt();
+    cholesky::LLT<real>     llt;
+    cholesky::LLTabs<real>  llta;
+    cholesky::LLTzero<real> lltz;
+    cholesky::LDLT<real>    ldlt;
 
+    unsigned int llt_fail = 0;
     std::vector<double> s0,s1,s2,s3;
     for( unsigned int n=0; n<R; ++n ) {
 
-        // a variance-covariance matrix of size N*N
-        CovarMat V(N,N);
-
-        // random covariance matrix
-        for( unsigned int i=0; i<N; ++i) {
-            V(i,i) = std::pow(rand(),2); // variance should be >= 0
-            for( unsigned int j=i+1; j<N; ++j) {
-                V(i,j) = rand();
+        // a random sample matrix
+        ublas::matrix<real> S(M,N);
+        for( unsigned int i=0; i<M; ++i) {
+            for( unsigned int j=0; j<N; ++j) {
+                S(i,j) = F * static_cast<real>(rand())/RAND_MAX;
             }
         }
+        
+        // a variance-covariance matrix of size N*N
+        CovarMat V = ublas::prod( ublas::trans(S), S );
+        assert( V.size1() == N && V.size2() == N );
 
-        FactorMat L0 = llt(V);
-        CovarMat V0 = ublas::prod( L0, ublas::trans(L0) );
-        s0.push_back( trigsum(error(V,V0)) );
+        if( R == 1 ) {
+            std::cout << std::endl << "Covariance matrix:" << std::endl;
+            std::cout << format(V) << std::endl;
+        }
+
+        FactorMat L0; 
+        try {
+            L0 = llt(V);
+            CovarMat V0 = ublas::prod( L0, ublas::trans(L0) );
+            s0.push_back( trigsum(error(V,V0)) );
+        } catch( cholesky::NotDefinitePositive & error ) {
+            llt_fail++;
+#ifndef NDEBUG
+            std::cout << "LLT FAILED:\t" << error.what() << std::endl;
+#endif
+        }
 
         FactorMat L1 = llta(V);
         CovarMat V1 = ublas::prod( L1, ublas::trans(L1) );
@@ -157,8 +208,16 @@ int main(int argc, char** argv)
     }
 
     std::cout << "Average error:" << std::endl;
-    std::cout << "\tLLT:  " << sum(s0)/R << std::endl;
+    
+    std::cout << "\tLLT:  ";
+    if( s0.size() == 0 ) {
+        std::cout << "NAN";
+    } else {
+        std::cout << sum(s0)/R;
+    }
+    std::cout << "\t" << llt_fail << "/" << R << " failures" << std::endl;
+
     std::cout << "\tLLTa: " << sum(s1)/R << std::endl;
     std::cout << "\tLLTz: " << sum(s2)/R << std::endl;
     std::cout << "\tLDLT: " << sum(s3)/R << std::endl;
-
+}