diff --git a/include/boost/numeric/ublas/matrix_sparse.hpp b/include/boost/numeric/ublas/matrix_sparse.hpp
index 9324108e..b395a529 100644
--- a/include/boost/numeric/ublas/matrix_sparse.hpp
+++ b/include/boost/numeric/ublas/matrix_sparse.hpp
@@ -1359,7 +1359,7 @@ namespace boost { namespace numeric { namespace ublas {
 #endif
         typedef const matrix_reference<const self_type> const_closure_type;
         typedef matrix_reference<self_type> closure_type;
-        typedef mapped_vector<T, typename A::value_type> vector_temporary_type;
+        typedef mapped_vector<T> vector_temporary_type;
         typedef self_type matrix_temporary_type;
         typedef typename A::value_type::second_type vector_data_value_type;
         typedef sparse_tag storage_category;
@@ -1530,7 +1530,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class C>          // Container assignment without temporary
         BOOST_UBLAS_INLINE
         mapped_vector_of_mapped_vector &operator = (const matrix_container<C> &m) {
-            resize (m ().size1 (), m ().size2 ());
+            resize (m ().size1 (), m ().size2 (), false);
             assign (m);
             return *this;
         }
diff --git a/include/boost/numeric/ublas/operation/num_columns.hpp b/include/boost/numeric/ublas/operation/num_columns.hpp
index 30fcc7b6..ee48eea6 100644
--- a/include/boost/numeric/ublas/operation/num_columns.hpp
+++ b/include/boost/numeric/ublas/operation/num_columns.hpp
@@ -20,6 +20,8 @@
 
 
 #include <boost/numeric/ublas/detail/config.hpp>
+#include <boost/numeric/ublas/expression_types.hpp>
+#include <boost/numeric/ublas/traits.hpp>
 
 
 namespace boost { namespace numeric { namespace ublas {
@@ -32,9 +34,9 @@ namespace boost { namespace numeric { namespace ublas {
      */
     template <typename MatrixExprT>
     BOOST_UBLAS_INLINE
-    typename MatrixExprT::size_type num_columns(MatrixExprT const& m)
+    typename matrix_traits<MatrixExprT>::size_type num_columns(matrix_expression<MatrixExprT> const& me)
     {
-        return m.size2();
+        return me().size2();
     }
 
 }}} // Namespace boost::numeric::ublas
diff --git a/include/boost/numeric/ublas/operation/num_rows.hpp b/include/boost/numeric/ublas/operation/num_rows.hpp
index 856256ff..246aed3c 100644
--- a/include/boost/numeric/ublas/operation/num_rows.hpp
+++ b/include/boost/numeric/ublas/operation/num_rows.hpp
@@ -5,7 +5,7 @@
  *
  * \brief The \c num_rows operation.
  *
- * Copyright (c) 2009, Marco Guazzone
+ * Copyright (c) 2009-2012, Marco Guazzone
  *
  * Distributed under the Boost Software License, Version 1.0. (See
  * accompanying file LICENSE_1_0.txt or copy at
@@ -19,6 +19,8 @@
 
 
 #include <boost/numeric/ublas/detail/config.hpp>
+#include <boost/numeric/ublas/expression_types.hpp>
+#include <boost/numeric/ublas/traits.hpp>
 
 
 namespace boost { namespace numeric { namespace ublas {
@@ -31,9 +33,9 @@ namespace boost { namespace numeric { namespace ublas {
      */
     template <typename MatrixExprT>
     BOOST_UBLAS_INLINE
-    typename MatrixExprT::size_type num_rows(MatrixExprT const& m)
+    typename matrix_traits<MatrixExprT>::size_type num_rows(matrix_expression<MatrixExprT> const& me)
     {
-        return m.size1();
+        return me().size1();
     }
 
 }}} // Namespace boost::numeric::ublas
diff --git a/include/boost/numeric/ublas/traits.hpp b/include/boost/numeric/ublas/traits.hpp
index 6964031c..ae020700 100644
--- a/include/boost/numeric/ublas/traits.hpp
+++ b/include/boost/numeric/ublas/traits.hpp
@@ -41,6 +41,14 @@ namespace {
     // we'll find either std::abs or else another version via ADL:
     return abs (t);
   }
+  // unsigned types are always non-negative
+  template<> unsigned int boost_numeric_ublas_abs (const unsigned int& t) {
+    return t;
+  }
+  // unsigned types are always non-negative
+  template<> unsigned long boost_numeric_ublas_abs (const unsigned long& t) {
+    return t;
+  }
 }
 
 namespace boost { namespace numeric { namespace ublas {
@@ -594,10 +602,10 @@ namespace boost { namespace numeric { namespace ublas {
         struct has_trivial_destructor : public boost::has_trivial_destructor<T> {};
 
         template<typename FLT>
-        struct has_trivial_constructor<std::complex<FLT> > : public boost::true_type {};
+        struct has_trivial_constructor<std::complex<FLT> > : public has_trivial_constructor<FLT> {};
         
         template<typename FLT>
-        struct has_trivial_destructor<std::complex<FLT> > : public boost::true_type {};
+        struct has_trivial_destructor<std::complex<FLT> > : public has_trivial_destructor<FLT> {};
 
     }
 
diff --git a/include/boost/numeric/ublas/triangular.hpp b/include/boost/numeric/ublas/triangular.hpp
index b537df96..cce143f2 100644
--- a/include/boost/numeric/ublas/triangular.hpp
+++ b/include/boost/numeric/ublas/triangular.hpp
@@ -1936,6 +1936,92 @@ namespace boost { namespace numeric { namespace ublas {
             }
         }
     }
+
+    // Dense (proxy) case
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
+                        lower_tag, row_major_tag, dense_proxy_tag) {
+        typedef typename E2::size_type size_type;
+        typedef typename E2::difference_type difference_type;
+        typedef typename E2::value_type value_type;
+
+        BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
+        BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
+        size_type size = e2 ().size ();
+        for (size_type n = 0; n < size; ++ n) {
+#ifndef BOOST_UBLAS_SINGULAR_CHECK
+            BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
+#else
+            if (e1 () (n, n) == value_type/*zero*/())
+                singular ().raise ();
+#endif
+            value_type t = e2 () (n) /= e1 () (n, n);
+            if (t != value_type/*zero*/()) {
+                for (size_type m = n + 1; m < size; ++ m)
+                    e2 () (m) -= e1 () (m, n) * t;
+            }
+        }
+    }
+    // Packed (proxy) case
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
+                        lower_tag, row_major_tag, packed_proxy_tag) {
+        typedef typename E2::size_type size_type;
+        typedef typename E2::difference_type difference_type;
+        typedef typename E2::value_type value_type;
+
+        BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
+        BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
+        size_type size = e2 ().size ();
+        for (size_type n = 0; n < size; ++ n) {
+#ifndef BOOST_UBLAS_SINGULAR_CHECK
+            BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
+#else
+            if (e1 () (n, n) == value_type/*zero*/())
+                singular ().raise ();
+#endif
+            value_type t = e2 () (n);
+            typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, 0));
+            typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, n));
+            while (it2e1 != it2e1_end) {
+              t -= *it2e1 * e2 () (it2e1.index2());
+              ++ it2e1;
+            }
+            e2() (n) = t / e1 () (n, n);
+        }
+    }
+    // Sparse (proxy) case
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
+                        lower_tag, row_major_tag, unknown_storage_tag) {
+        typedef typename E2::size_type size_type;
+        typedef typename E2::difference_type difference_type;
+        typedef typename E2::value_type value_type;
+
+        BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
+        BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
+        size_type size = e2 ().size ();
+        for (size_type n = 0; n < size; ++ n) {
+#ifndef BOOST_UBLAS_SINGULAR_CHECK
+            BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
+#else
+            if (e1 () (n, n) == value_type/*zero*/())
+                singular ().raise ();
+#endif
+            value_type t = e2 () (n);
+            typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, 0));
+            typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, n));
+            while (it2e1 != it2e1_end) {
+              t -= *it2e1 * e2 () (it2e1.index2());
+              ++ it2e1;
+            }
+            e2() (n) = t / e1 () (n, n);
+        }
+    }
+
     // Redirectors :-)
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
@@ -1950,8 +2036,8 @@ namespace boost { namespace numeric { namespace ublas {
     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
                         lower_tag, row_major_tag) {
         typedef typename E1::storage_category storage_category;
-        inplace_solve (e2, trans (e1),
-                       upper_tag (), row_major_tag (), storage_category ());
+        inplace_solve (e1, e2,
+                       lower_tag (), row_major_tag (), storage_category ());
     }
     // Dispatcher
     template<class E1, class E2>
@@ -2020,9 +2106,9 @@ namespace boost { namespace numeric { namespace ublas {
             if (t != value_type/*zero*/()) {
                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
-                difference_type m (it1e1_rend - it1e1);
-                while (-- m >= 0)
-                    e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
+                while (it1e1 != it1e1_rend) {
+                  e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
+                }
             }
         }
     }
@@ -2049,11 +2135,100 @@ namespace boost { namespace numeric { namespace ublas {
             if (t != value_type/*zero*/()) {
                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
-                while (it1e1 != it1e1_rend)
-                    e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
+                while (it1e1 != it1e1_rend) {
+                  e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
+                }
             }
         }
     }
+
+    // Dense (proxy) case
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
+                        upper_tag, row_major_tag, dense_proxy_tag) {
+        typedef typename E2::size_type size_type;
+        typedef typename E2::difference_type difference_type;
+        typedef typename E2::value_type value_type;
+
+        BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
+        BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
+        size_type size = e1 ().size1 ();
+        for (difference_type n = size-1; n >=0; -- n) {
+#ifndef BOOST_UBLAS_SINGULAR_CHECK
+            BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
+#else
+            if (e1 () (n, n) == value_type/*zero*/())
+                singular ().raise ();
+#endif
+            value_type t = e2 () (n);
+            for (difference_type m = n + 1; m < e1 ().size2(); ++ m) {
+              t -= e1 () (n, m)  * e2 () (m);
+            }
+            e2() (n) = t / e1 () (n, n);
+        }
+    }
+    // Packed (proxy) case
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
+                        upper_tag, row_major_tag, packed_proxy_tag) {
+        typedef typename E2::size_type size_type;
+        typedef typename E2::difference_type difference_type;
+        typedef typename E2::value_type value_type;
+
+        BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
+        BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
+        size_type size = e1 ().size1 ();
+        for (difference_type n = size-1; n >=0; -- n) {
+#ifndef BOOST_UBLAS_SINGULAR_CHECK
+            BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
+#else
+            if (e1 () (n, n) == value_type/*zero*/())
+                singular ().raise ();
+#endif
+            value_type t = e2 () (n);
+            typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, n+1));
+            typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, e1 ().size2 ()));
+            while (it2e1 != it2e1_end) {
+              t -= *it2e1 * e2 () (it2e1.index2());
+              ++ it2e1;
+            }
+            e2() (n) = t / e1 () (n, n);
+
+        }
+    }
+    // Sparse (proxy) case
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
+                        upper_tag, row_major_tag, unknown_storage_tag) {
+        typedef typename E2::size_type size_type;
+        typedef typename E2::difference_type difference_type;
+        typedef typename E2::value_type value_type;
+
+        BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
+        BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
+        size_type size = e1 ().size1 ();
+        for (difference_type n = size-1; n >=0; -- n) {
+#ifndef BOOST_UBLAS_SINGULAR_CHECK
+            BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
+#else
+            if (e1 () (n, n) == value_type/*zero*/())
+                singular ().raise ();
+#endif
+            value_type t = e2 () (n);
+            typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, n+1));
+            typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, e1 ().size2 ()));
+            while (it2e1 != it2e1_end) {
+              t -= *it2e1 * e2 () (it2e1.index2());
+              ++ it2e1;
+            }
+            e2() (n) = t / e1 () (n, n);
+
+        }
+    }
+
     // Redirectors :-)
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
@@ -2068,8 +2243,8 @@ namespace boost { namespace numeric { namespace ublas {
     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
                         upper_tag, row_major_tag) {
         typedef typename E1::storage_category storage_category;
-        inplace_solve (e2, trans (e1),
-                       lower_tag (), row_major_tag (), storage_category ());
+        inplace_solve (e1, e2,
+                       upper_tag (), row_major_tag (), storage_category ());
     }
     // Dispatcher
     template<class E1, class E2>
@@ -2100,103 +2275,21 @@ namespace boost { namespace numeric { namespace ublas {
         return r;
     }
 
-    // Dense (proxy) case
-    template<class E1, class E2>
-    BOOST_UBLAS_INLINE
-    void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
-                        lower_tag, row_major_tag, dense_proxy_tag) {
-        typedef typename E1::size_type size_type;
-        typedef typename E1::difference_type difference_type;
-        typedef typename E1::value_type value_type;
 
-        BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
-        BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
-        size_type size = e1 ().size ();
-        for (difference_type n = size - 1; n >= 0; -- n) {
-#ifndef BOOST_UBLAS_SINGULAR_CHECK
-            BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
-#else
-            if (e2 () (n, n) == value_type/*zero*/())
-                singular ().raise ();
-#endif
-            value_type t = e1 () (n) /= e2 () (n, n);
-            if (t != value_type/*zero*/()) {
-                for (difference_type m = n - 1; m >= 0; -- m)
-                    e1 () (m) -= t * e2 () (n, m);
-            }
-        }
-    }
-    // Packed (proxy) case
-    template<class E1, class E2>
-    BOOST_UBLAS_INLINE
-    void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
-                        lower_tag, row_major_tag, packed_proxy_tag) {
-        typedef typename E1::size_type size_type;
-        typedef typename E1::difference_type difference_type;
-        typedef typename E1::value_type value_type;
-
-        BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
-        BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
-        size_type size = e1 ().size ();
-        for (difference_type n = size - 1; n >= 0; -- n) {
-#ifndef BOOST_UBLAS_SINGULAR_CHECK
-            BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
-#else
-            if (e2 () (n, n) == value_type/*zero*/())
-                singular ().raise ();
-#endif
-            value_type t = e1 () (n) /= e2 () (n, n);
-            if (t != value_type/*zero*/()) {
-                typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));
-                typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));
-                difference_type m (it2e2_rend - it2e2);
-                while (-- m >= 0)
-                    e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
-            }
-        }
-    }
-    // Sparse (proxy) case
-    template<class E1, class E2>
-    BOOST_UBLAS_INLINE
-    void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
-                        lower_tag, row_major_tag, unknown_storage_tag) {
-        typedef typename E1::size_type size_type;
-        typedef typename E1::difference_type difference_type;
-        typedef typename E1::value_type value_type;
-
-        BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
-        BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
-        size_type size = e1 ().size ();
-        for (difference_type n = size - 1; n >= 0; -- n) {
-#ifndef BOOST_UBLAS_SINGULAR_CHECK
-            BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
-#else
-            if (e2 () (n, n) == value_type/*zero*/())
-                singular ().raise ();
-#endif
-            value_type t = e1 () (n) /= e2 () (n, n);
-            if (t != value_type/*zero*/()) {
-                typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));
-                typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));
-                while (it2e2 != it2e2_rend)
-                    e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
-            }
-        }
-    }
     // Redirectors :-)
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
                         lower_tag, row_major_tag) {
-        typedef typename E1::storage_category storage_category;
-        inplace_solve (e1, e2,
-                       lower_tag (), row_major_tag (), storage_category ());
+        typedef typename E2::storage_category storage_category;
+        inplace_solve (trans(e2), e1,
+                       upper_tag (), column_major_tag (), storage_category ());
     }
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
                         lower_tag, column_major_tag) {
-        typedef typename E1::storage_category storage_category;
+        typedef typename E2::storage_category storage_category;
         inplace_solve (trans (e2), e1,
                        upper_tag (), row_major_tag (), storage_category ());
     }
@@ -2218,103 +2311,21 @@ namespace boost { namespace numeric { namespace ublas {
                        unit_lower_tag (), orientation_category ());
     }
 
-    // Dense (proxy) case
-    template<class E1, class E2>
-    BOOST_UBLAS_INLINE
-    void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
-                        upper_tag, row_major_tag, dense_proxy_tag) {
-        typedef typename E1::size_type size_type;
-        typedef typename E1::difference_type difference_type;
-        typedef typename E1::value_type value_type;
 
-        BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
-        BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
-        size_type size = e1 ().size ();
-        for (size_type n = 0; n < size; ++ n) {
-#ifndef BOOST_UBLAS_SINGULAR_CHECK
-            BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
-#else
-            if (e2 () (n, n) == value_type/*zero*/())
-                singular ().raise ();
-#endif
-            value_type t = e1 () (n) /= e2 () (n, n);
-            if (t != value_type/*zero*/()) {
-                for (size_type m = n + 1; m < size; ++ m)
-                    e1 () (m) -= t * e2 () (n, m);
-            }
-        }
-    }
-    // Packed (proxy) case
-    template<class E1, class E2>
-    BOOST_UBLAS_INLINE
-    void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
-                        upper_tag, row_major_tag, packed_proxy_tag) {
-        typedef typename E1::size_type size_type;
-        typedef typename E1::difference_type difference_type;
-        typedef typename E1::value_type value_type;
-
-        BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
-        BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
-        size_type size = e1 ().size ();
-        for (size_type n = 0; n < size; ++ n) {
-#ifndef BOOST_UBLAS_SINGULAR_CHECK
-            BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
-#else
-            if (e2 () (n, n) == value_type/*zero*/())
-                singular ().raise ();
-#endif
-            value_type t = e1 () (n) /= e2 () (n, n);
-            if (t != value_type/*zero*/()) {
-                typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));
-                typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));
-                difference_type m (it2e2_end - it2e2);
-                while (-- m >= 0)
-                    e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
-            }
-        }
-    }
-    // Sparse (proxy) case
-    template<class E1, class E2>
-    BOOST_UBLAS_INLINE
-    void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
-                        upper_tag, row_major_tag, unknown_storage_tag) {
-        typedef typename E1::size_type size_type;
-        typedef typename E1::difference_type difference_type;
-        typedef typename E1::value_type value_type;
-
-        BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
-        BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
-        size_type size = e1 ().size ();
-        for (size_type n = 0; n < size; ++ n) {
-#ifndef BOOST_UBLAS_SINGULAR_CHECK
-            BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
-#else
-            if (e2 () (n, n) == value_type/*zero*/())
-                singular ().raise ();
-#endif
-            value_type t = e1 () (n) /= e2 () (n, n);
-            if (t != value_type/*zero*/()) {
-                typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));
-                typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));
-                while (it2e2 != it2e2_end)
-                    e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
-            }
-        }
-    }
     // Redirectors :-)
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
                         upper_tag, row_major_tag) {
-        typedef typename E1::storage_category storage_category;
-        inplace_solve (e1, e2,
-                       upper_tag (), row_major_tag (), storage_category ());
+        typedef typename E2::storage_category storage_category;
+        inplace_solve (trans(e2), e1,
+                       lower_tag (), column_major_tag (), storage_category ());
     }
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
                         upper_tag, column_major_tag) {
-        typedef typename E1::storage_category storage_category;
+        typedef typename E2::storage_category storage_category;
         inplace_solve (trans (e2), e1,
                        lower_tag (), row_major_tag (), storage_category ());
     }
diff --git a/include/boost/numeric/ublas/vector_expression.hpp b/include/boost/numeric/ublas/vector_expression.hpp
index a5cf3104..fe4d1c3f 100644
--- a/include/boost/numeric/ublas/vector_expression.hpp
+++ b/include/boost/numeric/ublas/vector_expression.hpp
@@ -1409,7 +1409,9 @@ namespace boost { namespace numeric { namespace ublas {
     // (v / t) [i] = v [i] / t
     template<class E1, class T2>
     BOOST_UBLAS_INLINE
+    typename enable_if< is_convertible<T2, typename E1::value_type >,    
     typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
+    >::type
     operator / (const vector_expression<E1> &e1,
                 const T2 &e2) {
         typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 388561b1..423d1999 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1,4 +1,4 @@
-# Copyright (c) 2004 Michael Stevens
+# Copyright (c) 2004-2011 Michael Stevens, David Bellot
 # Use, modification and distribution are subject to the
 # Boost Software License, Version 1.0. (See accompanying file
 # LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
@@ -31,12 +31,20 @@ UBLAS_TESTSET_SPARSE ?=
             # USE_RANGE USE_SLICE	 # Too complex for regression testing
             USE_UNBOUNDED_ARRAY
 			USE_MAP_ARRAY USE_STD_MAP
-            USE_MAPPED_VECTOR USE_COMPRESSED_VECTOR USE_COORDINATE_VECTOR
-            USE_MAPPED_MATRIX USE_COMPRESSED_MATRIX USE_COORDINATE_MATRIX 
+            USE_MAPPED_VECTOR USE_COMPRESSED_VECTOR 
+            USE_MAPPED_MATRIX USE_COMPRESSED_MATRIX 
 			;
 # Generalize VofV still failing
 #            USE_GENERALIZED_VECTOR_OF_VECTOR
 
+UBLAS_TESTSET_SPARSE_COO = [ modules.peek : UBLAS_TESTSET_SPARSE_COO ] ;
+UBLAS_TESTSET_SPARSE_COO ?= 
+            USE_DOUBLE USE_STD_COMPLEX
+            USE_UNBOUNDED_ARRAY
+            USE_COORDINATE_VECTOR
+            USE_COORDINATE_MATRIX 
+			;
+
 
 # Project settings
 project
@@ -69,6 +77,20 @@ test-suite numeric/uBLAS
             test33.cpp
         : : :
             <define>$(UBLAS_TESTSET_SPARSE)
+            <define>$(UBLAS_TESTSET_SPARSE_COO)
+      ]
+      [ run test3.cpp
+            test31.cpp
+            test32.cpp
+            test33.cpp
+        : : :
+            <define>USE_FLOAT
+            <define>USE_DOUBLE 
+            <define>USE_STD_COMPLEX 
+            <define>USE_STD_MAP
+            <define>USE_MAPPED_VECTOR_OF_MAPPED_VECTOR 
+        : test3_mvov 
+        :
       ]
       [ run test4.cpp
             test42.cpp
@@ -132,6 +154,31 @@ test-suite numeric/uBLAS
       ]
       [ run test_complex_norms.cpp
       ]
-	[ run test_assignment.cpp
-	]
+	  [ run test_assignment.cpp
+	  ]
+	  [ run test_triangular.cpp
+	  ]
+	  [ run test_ticket7296.cpp
+	  ]
+      [ run test_inplace_solve.cpp
+        : 
+        : 
+        : <define>$(UBLAS_TESTSET) 
+        : test_inplace_solve_basic
+        :
+      ]
+      [ run test_inplace_solve.cpp
+        : 
+        : 
+        : <define>$(UBLAS_TESTSET_SPARSE) <define>$(UBLAS_TESTSET_SPARSE_COO)
+        : test_inplace_solve_sparse
+        :
+      ]
+      [ run test_inplace_solve.cpp
+        : 
+        : 
+        : <define>USE_MAPPED_VECTOR_OF_MAPPED_VECTOR
+        : test_inplace_solve_mvov
+        :
+      ]
     ;
diff --git a/test/README b/test/README
index 2e3e7177..76112314 100644
--- a/test/README
+++ b/test/README
@@ -1,4 +1,4 @@
-Copyright (c) 2000-2004 Joerg Walter, Mathias Koch
+Copyright (c) 2000-2011 Joerg Walter, Mathias Koch, David Bellot
 
 Distributed under the Boost Software License, Version 1.0. (See
 accompanying file LICENSE_1_0.txt or copy at
diff --git a/test/begin_end.cpp b/test/begin_end.cpp
index db4d8246..ff924ed0 100644
--- a/test/begin_end.cpp
+++ b/test/begin_end.cpp
@@ -45,7 +45,7 @@ BOOST_UBLAS_TEST_DEF( test_vector_iteration )
             ++it
     ) {
         BOOST_UBLAS_DEBUG_TRACE( "*it = " << *it << " ==> " << v(ix) );
-        BOOST_UBLAS_TEST_CHECK( std::fabs(*it - v(ix)) <= TOL );
+        BOOST_UBLAS_TEST_CHECK( std::abs(*it - v(ix)) <= TOL );
         ++ix;
     }
 }
@@ -74,7 +74,7 @@ BOOST_UBLAS_TEST_DEF( test_vector_const_iteration )
             ++it
     ) {
         BOOST_UBLAS_DEBUG_TRACE( "*it = " << *it << " ==> " << v(ix) );
-        BOOST_UBLAS_TEST_CHECK( std::fabs(*it - v(ix)) <= TOL );
+        BOOST_UBLAS_TEST_CHECK( std::abs(*it - v(ix)) <= TOL );
         ++ix;
     }
 }
@@ -112,7 +112,7 @@ BOOST_UBLAS_TEST_DEF( test_row_major_matrix_iteration )
                 ++inner_it
         ) {
             BOOST_UBLAS_DEBUG_TRACE( "*it = " << *inner_it << " ==> " << A(row,col) );
-            BOOST_UBLAS_TEST_CHECK( std::fabs(*inner_it - A(row,col)) <= TOL );
+            BOOST_UBLAS_TEST_CHECK( std::abs(*inner_it - A(row,col)) <= TOL );
 
             ++col;
         }
@@ -154,7 +154,7 @@ BOOST_UBLAS_TEST_DEF( test_col_major_matrix_iteration )
                 ++inner_it
         ) {
             BOOST_UBLAS_DEBUG_TRACE( "*it = " << *inner_it << " ==> " << A(row,col) );
-            BOOST_UBLAS_TEST_CHECK( std::fabs(*inner_it - A(row,col)) <= TOL );
+            BOOST_UBLAS_TEST_CHECK( std::abs(*inner_it - A(row,col)) <= TOL );
 
             ++row;
         }
diff --git a/test/test32.cpp b/test/test32.cpp
index 2ffc5021..eaf49297 100644
--- a/test/test32.cpp
+++ b/test/test32.cpp
@@ -321,4 +321,34 @@ void test_matrix_vector () {
 #endif
 #endif
 #endif
+
+#ifdef USE_MAPPED_VECTOR_OF_MAPPED_VECTOR
+#ifdef USE_STD_MAP
+#ifdef USE_FLOAT
+    std::cout << "float mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix_vector<ublas::mapped_vector<float, std::map<std::size_t, float> >,
+                          ublas::mapped_vector_of_mapped_vector<float>, 3 > () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "double mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix_vector<ublas::mapped_vector<double, std::map<std::size_t, double> >,
+                          ublas::mapped_vector_of_mapped_vector<double>, 3 > () ();
+#endif
+
+#ifdef USE_STD_COMPLEX
+#ifdef USE_FLOAT
+    std::cout << "std::complex<float> mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix_vector<ublas::mapped_vector<std::complex<float>, std::map<std::size_t, std::complex<float> > >,
+                          ublas::mapped_vector_of_mapped_vector<std::complex<float> >, 3 > () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "std::complex<double> mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix_vector<ublas::mapped_vector<std::complex<double>, std::map<std::size_t, std::complex<double> > >,
+                          ublas::mapped_vector_of_mapped_vector<std::complex<double> >, 3 > () ();
+#endif
+#endif
+#endif
+#endif
 }
diff --git a/test/test33.cpp b/test/test33.cpp
index bf7a3d28..260d8b9a 100644
--- a/test/test33.cpp
+++ b/test/test33.cpp
@@ -344,4 +344,28 @@ void test_matrix () {
 #endif
 #endif
 #endif
+
+#ifdef USE_MAPPED_VECTOR_OF_MAPPED_VECTOR
+#ifdef USE_FLOAT
+    std::cout << "float mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix<ublas::mapped_vector_of_mapped_vector<float>, 3 > () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "double mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix<ublas::mapped_vector_of_mapped_vector<double>, 3 > () ();
+#endif
+
+#ifdef USE_STD_COMPLEX
+#ifdef USE_FLOAT
+    std::cout << "std::complex<float> mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix<ublas::mapped_vector_of_mapped_vector<std::complex<float> >, 3 > () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "std::complex<double> mapped_vector_of_mapped_vector" << std::endl;
+    test_my_matrix<ublas::mapped_vector_of_mapped_vector<std::complex<double> >, 3 > () ();
+#endif
+#endif
+#endif
 }
diff --git a/test/test_assignment.cpp b/test/test_assignment.cpp
index cbfe0e56..4a3897c7 100644
--- a/test/test_assignment.cpp
+++ b/test/test_assignment.cpp
@@ -26,7 +26,7 @@ typename AE::value_type mean_square(const matrix_expression<AE> &me) {
     typename AE::size_type i, j;
     for (i=0; i!= me().size1(); i++) {
         for (j=0; j!= me().size2(); j++) {
-            s+=std::fabs(me()(i,j));
+          s+= scalar_traits<typename AE::value_type>::type_abs(me()(i,j));
         }
     }
     return s/me().size1()*me().size2();
@@ -39,7 +39,7 @@ typename AE::value_type mean_square(const vector_expression<AE> &ve) {
     typename AE::value_type s(0);
     typename AE::size_type i;
     for (i=0; i!= ve().size(); i++) {
-            s+=std::fabs(ve()(i));
+            s+=scalar_traits<typename AE::value_type>::type_abs(ve()(i));
         }
     return s/ve().size();
 }
diff --git a/test/test_ticket7296.cpp b/test/test_ticket7296.cpp
new file mode 100644
index 00000000..05c042b6
--- /dev/null
+++ b/test/test_ticket7296.cpp
@@ -0,0 +1,292 @@
+/**
+ * \file libs/numeric/ublas/test/test_utils.hpp
+ *
+ * \brief Test suite for utils.hpp.
+ *
+ * Copyright (c) 2012, Marco Guazzone
+ * 
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * \author Marco Guazzone (marco.guazzone@gmail.com)
+ */
+
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <complex>
+#include <cstddef>
+#include <libs/numeric/ublas/test/utils.hpp>
+
+
+namespace ublas = boost::numeric::ublas;
+
+
+static const double tol(1e-6);
+static const double mul(tol*10);
+
+
+BOOST_UBLAS_TEST_DEF( check )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check'" );
+
+	BOOST_UBLAS_TEST_CHECK( true );
+}
+
+BOOST_UBLAS_TEST_DEF( check_eq )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_eq'" );
+
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_EQ( short(1), short(1) );
+	BOOST_UBLAS_TEST_CHECK_EQ( int(1), int(1) );
+	BOOST_UBLAS_TEST_CHECK_EQ( long(1), long(1) );
+	BOOST_UBLAS_TEST_CHECK_EQ( unsigned(1), unsigned(1) );
+
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_EQ( short(1), int(1) );
+	BOOST_UBLAS_TEST_CHECK_EQ( short(1), int(1) );
+	BOOST_UBLAS_TEST_CHECK_EQ( int(1), long(1) );
+	BOOST_UBLAS_TEST_CHECK_EQ( long(1), int(1) );
+
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test aliases." );
+	BOOST_UBLAS_TEST_CHECK_EQUAL( int(1), int(1) );
+}
+
+BOOST_UBLAS_TEST_DEF( check_close )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_close'" );
+
+	const double c1(1*mul);
+	const double c2(2*mul);
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( float(c1), float(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( double(c1), double(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( ::std::complex<float>(c1,c2), ::std::complex<float>(c1,c2), tol );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( ::std::complex<double>(c1,c2), ::std::complex<double>(c1,c2), tol );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( float(c1), double(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( double(c1), float(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( ::std::complex<float>(c1,c2), ::std::complex<double>(c1,c2), tol );
+	BOOST_UBLAS_TEST_CHECK_CLOSE( ::std::complex<double>(c1,c2), ::std::complex<float>(c1,c2), tol );
+
+	// Check alias
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test aliases." );
+	BOOST_UBLAS_TEST_CHECK_PRECISION( float(c1), float(c1), tol );
+}
+
+BOOST_UBLAS_TEST_DEF( check_rel_close )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_rel_close'" );
+
+	const double c1(1*mul);
+	const double c2(2*mul);
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( float(c1), float(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( double(c1), double(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( ::std::complex<float>(c1,c2), ::std::complex<float>(c1,c2), tol );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( ::std::complex<double>(c1,c2), ::std::complex<double>(c1,c2), tol );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( float(c1), double(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( double(c1), float(c1), tol );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( ::std::complex<float>(c1,c2), ::std::complex<double>(c1,c2), tol );
+	BOOST_UBLAS_TEST_CHECK_REL_CLOSE( ::std::complex<double>(c1,c2), ::std::complex<float>(c1,c2), tol );
+
+	// Check alias
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test aliases." );
+	BOOST_UBLAS_TEST_CHECK_REL_PRECISION( float(c1), float(c1), tol );
+}
+
+BOOST_UBLAS_TEST_DEF( check_vector_eq )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_vector_eq'" );
+
+	const ::std::size_t n(5);
+
+	ublas::vector<short> sv = ublas::scalar_vector<short>(n, 1);
+	ublas::vector<int> iv = ublas::scalar_vector<int>(n, 1);
+	ublas::vector<long> lv = ublas::scalar_vector<long>(n, 1L);
+	ublas::vector<unsigned> uv = ublas::scalar_vector<unsigned>(n, 1u);
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( sv, sv, n );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( iv, iv, n );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( lv, lv, n );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( uv, uv, n );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( sv, iv, n );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( iv, sv, n );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( iv, lv, n );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_EQ( lv, iv, n );
+}
+
+BOOST_UBLAS_TEST_DEF( check_vector_close )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_vector_close'" );
+
+	const ::std::size_t n(5);
+
+	ublas::vector<float> fv = ublas::scalar_vector<float>(n, 1);
+	ublas::vector<float> dv = ublas::scalar_vector<float>(n, 1);
+	ublas::vector< ::std::complex<float> > cfv = ublas::scalar_vector< ::std::complex<float> >(n, ::std::complex<float>(1,2));
+	ublas::vector< ::std::complex<double> > cdv = ublas::scalar_vector< ::std::complex<double> >(n, ::std::complex<double>(1,2));
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( fv, fv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( dv, dv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( cfv, cfv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( cdv, cdv, n, tol );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( fv, dv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( dv, fv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( cfv, cdv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE( cdv, cfv, n, tol );
+}
+
+BOOST_UBLAS_TEST_DEF( check_vector_rel_close )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_vector_rel_close'" );
+
+	const ::std::size_t n(5);
+	const double c1(1*mul);
+	const double c2(2*mul);
+
+	ublas::vector<float> fv = ublas::scalar_vector<float>(n, c1);
+	ublas::vector<double> dv = ublas::scalar_vector<double>(n, c1);
+	ublas::vector< ::std::complex<float> > cfv = ublas::scalar_vector< ::std::complex<float> >(n, ::std::complex<float>(c1,c2));
+	ublas::vector< ::std::complex<double> > cdv = ublas::scalar_vector< ::std::complex<double> >(n, ::std::complex<double>(c1,c2));
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( fv, fv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( dv, dv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( cfv, cfv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( cdv, cdv, n, tol );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( fv, dv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( dv, fv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( cfv, cdv, n, tol );
+	BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE( cdv, cfv, n, tol );
+}
+
+BOOST_UBLAS_TEST_DEF( check_matrix_eq )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_matrix_eq'" );
+
+	const ::std::size_t nr(3);
+	const ::std::size_t nc(4);
+
+	ublas::matrix<short> sv = ublas::scalar_matrix<short>(nr, nc, 1);
+	ublas::matrix<int> iv = ublas::scalar_matrix<int>(nr, nc, 1);
+	ublas::matrix<long> lv = ublas::scalar_matrix<long>(nr, nc, 1L);
+	ublas::matrix<unsigned> uv = ublas::scalar_matrix<unsigned>(nr, nc, 1u);
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( sv, sv, nr, nc );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( iv, iv, nr, nc );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( lv, lv, nr, nc );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( uv, uv, nr, nc );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( sv, iv, nr, nc );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( iv, sv, nr, nc );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( iv, lv, nr, nc );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_EQ( lv, iv, nr, nc );
+}
+
+BOOST_UBLAS_TEST_DEF( check_matrix_close )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_matrix_close'" );
+
+	const ::std::size_t nr(3);
+	const ::std::size_t nc(4);
+	const double c1(1*mul);
+	const double c2(2*mul);
+
+	ublas::matrix<float> fA = ublas::scalar_matrix<float>(nr, nc, c1);
+	ublas::matrix<double> dA = ublas::scalar_matrix<double>(nr, nc, c1);
+	ublas::matrix< ::std::complex<float> > cfA = ublas::scalar_matrix< ::std::complex<float> >(nr, nc, ::std::complex<float>(c1,c2));
+	ublas::matrix< ::std::complex<double> > cdA = ublas::scalar_matrix< ::std::complex<double> >(nr, nc, ::std::complex<double>(c1,c2));
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( fA, fA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( dA, dA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( cfA, cfA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( cdA, cdA, nr, nc, tol );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( fA, dA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( dA, fA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( cfA, cdA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE( cdA, cfA, nr, nc, tol );
+}
+
+
+BOOST_UBLAS_TEST_DEF( check_matrix_rel_close )
+{
+	BOOST_UBLAS_TEST_TRACE( "Test case: 'check_matrix_rel_close'" );
+
+	const ::std::size_t nr(3);
+	const ::std::size_t nc(4);
+	const double c1(1*mul);
+	const double c2(2*mul);
+
+	ublas::matrix<float> fA = ublas::scalar_matrix<float>(nr, nc, c1);
+	ublas::matrix<double> dA = ublas::scalar_matrix<double>(nr, nc, c1);
+	ublas::matrix< ::std::complex<float> > cfA = ublas::scalar_matrix< ::std::complex<float> >(nr, nc, ::std::complex<float>(c1,c2));
+	ublas::matrix< ::std::complex<double> > cdA = ublas::scalar_matrix< ::std::complex<double> >(nr, nc, ::std::complex<double>(c1,c2));
+
+	// Check T vs. T
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against same types." );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( fA, fA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( dA, dA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( cfA, cfA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( cdA, cdA, nr, nc, tol );
+
+	// Check T1 vs. T2
+	BOOST_UBLAS_DEBUG_TRACE( "-- Test against different types." );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( fA, dA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( dA, fA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( cfA, cdA, nr, nc, tol );
+	BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE( cdA, cfA, nr, nc, tol );
+}
+
+
+int main()
+{
+	BOOST_UBLAS_TEST_SUITE( "Test 'utils.hpp' functionalities" );
+
+	BOOST_UBLAS_TEST_BEGIN();
+		BOOST_UBLAS_TEST_DO( check );
+		BOOST_UBLAS_TEST_DO( check_eq );
+		BOOST_UBLAS_TEST_DO( check_close );
+		BOOST_UBLAS_TEST_DO( check_rel_close );
+		BOOST_UBLAS_TEST_DO( check_vector_eq );
+		BOOST_UBLAS_TEST_DO( check_vector_close );
+		BOOST_UBLAS_TEST_DO( check_vector_rel_close );
+		BOOST_UBLAS_TEST_DO( check_matrix_eq );
+		BOOST_UBLAS_TEST_DO( check_matrix_close );
+		BOOST_UBLAS_TEST_DO( check_matrix_rel_close );
+	BOOST_UBLAS_TEST_END();
+}
diff --git a/test/test_triangular.cpp b/test/test_triangular.cpp
new file mode 100644
index 00000000..56b1d344
--- /dev/null
+++ b/test/test_triangular.cpp
@@ -0,0 +1,139 @@
+#include <iostream>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/numeric/ublas/io.hpp>
+
+
+#include <sys/times.h>
+#include <sys/time.h>
+#include <sys/stat.h>
+
+
+namespace ublas  = boost::numeric::ublas;
+
+template<class mat, class vec>
+double diff(const mat& A, const vec& x, const vec& b) {
+  vec temp(prod(A, x) - b);
+  double result = 0;
+  for (typename vec::size_type i=0; i<temp.size(); ++i) {
+    result += temp(i)*temp(i);
+  }
+  return sqrt(result);
+}
+
+template<class mat, class vec>
+double diff(const vec& x, const mat& A, const vec& b) {
+  return diff(trans(A), x, b);
+}
+
+namespace ublas  = boost::numeric::ublas;
+
+
+int main() {
+  const int n=10000;
+#if 1
+  ublas::compressed_matrix<double, ublas::row_major>     mat_row_upp(n, n);
+  ublas::compressed_matrix<double, ublas::column_major>  mat_col_upp(n, n);
+  ublas::compressed_matrix<double, ublas::row_major>     mat_row_low(n, n);
+  ublas::compressed_matrix<double, ublas::column_major>  mat_col_low(n, n);
+#else
+  ublas::matrix<double, ublas::row_major>     mat_row_upp(n, n, 0);
+  ublas::matrix<double, ublas::column_major>  mat_col_upp(n, n, 0);
+  ublas::matrix<double, ublas::row_major>     mat_row_low(n, n, 0);
+  ublas::matrix<double, ublas::column_major>  mat_col_low(n, n, 0);
+#endif
+  ublas::vector<double>  b(n, 1);
+
+  std::cerr << "Constructing..." << std::endl;
+  for (int i=0; i<n; ++i) {
+    b(i) = rand() % 10;
+    double main = -10 + rand() % 20 ;
+    if (main == 0) main+=1;
+    double side = -10 + rand() % 20 ;
+    if (i-1>=0) {
+      mat_row_low(i, i-1) = side;
+    }
+    mat_row_low(i, i) = main;
+
+    mat_col_low(i, i) = main;
+    if (i+1<n) {
+      mat_col_low(i+1, i) = side;
+    }
+
+    mat_row_upp(i, i) = main;
+    if (i+1<n) {
+      mat_row_upp(i, i+1) = side;
+    }
+
+    if (i-1>=0) {
+      mat_col_upp(i-1, i) = side;
+    }
+    mat_col_upp(i, i) = main;
+  }
+
+  std::cerr << "Starting..." << std::endl;
+  {
+    clock_t start = clock();
+    ublas::vector<double>  x(b);
+    ublas::inplace_solve(mat_col_low,  x, ublas::lower_tag());
+    std::cerr << "Col_low x: " << clock()-start << std::endl;
+    std::cerr << "delta: " << diff(mat_col_low, x, b) << "\n";
+  }
+  {
+    clock_t start = clock();
+    ublas::vector<double>  x(b);
+    ublas::inplace_solve(mat_row_low, x, ublas::lower_tag());
+    std::cerr << "Row_low x: " << clock()-start << std::endl;
+    std::cerr << "delta: " << diff(mat_row_low, x, b) << "\n";
+  }
+
+  {
+    clock_t start = clock();
+    ublas::vector<double>  x(b);
+    ublas::inplace_solve(mat_col_upp,  x, ublas::upper_tag());
+    std::cerr << "col_upp x: " << clock()-start << std::endl;
+    std::cerr << "delta: " << diff(mat_col_upp, x, b) << "\n";
+  }
+  {
+    clock_t start = clock();
+    ublas::vector<double>  x(b);
+    ublas::inplace_solve(mat_row_upp, x, ublas::upper_tag());
+    std::cerr << "row_upp x: " << clock()-start << std::endl;
+    std::cerr << "delta: " << diff(mat_row_upp, x, b) << "\n";
+  }
+
+//   {
+//     clock_t start = clock();
+//     ublas::vector<double>  x(b);
+//     ublas::inplace_solve(x, mat_col_low, ublas::lower_tag());
+//     std::cerr << "x col_low: " << clock()-start << std::endl;
+//     std::cerr << "delta: " << diff(x, mat_col_low, b) << "\n";
+//   }
+  {
+    clock_t start = clock();
+    ublas::vector<double>  x(b);
+    ublas::inplace_solve(x, mat_row_low, ublas::lower_tag());
+    std::cerr << "x row_low: " << clock()-start << std::endl;
+    std::cerr << "delta: " << diff(x, mat_row_low, b) << "\n";
+  }
+
+//   {
+//     clock_t start = clock();
+//     ublas::vector<double>  x(b);
+//     ublas::inplace_solve(x, mat_col_upp, ublas::upper_tag());
+//     std::cerr << "x col_upp: " << clock()-start << std::endl;
+//     std::cerr << "delta: " << diff(x, mat_col_upp, b) << "\n";
+//   }
+  {
+    clock_t start = clock();
+    ublas::vector<double>  x(b);
+    ublas::inplace_solve(x, mat_row_upp, ublas::upper_tag());
+    std::cerr << "x row_upp: " << clock()-start << std::endl;
+    std::cerr << "delta: " << diff(x, mat_row_upp, b) << "\n";
+  }
+
+
+}
diff --git a/test/utils.hpp b/test/utils.hpp
index c0526ce8..ba14e860 100644
--- a/test/utils.hpp
+++ b/test/utils.hpp
@@ -1,34 +1,383 @@
-/** -*- c++ -*- \file utils.hpp \brief Test utilities. */
+/**
+ * \file util.hpp
+ *
+ * \brief Utility macros/functions for testing and debugging purpose.
+ *
+ * Basic usage:
+ * <pre>
+ * BOOST_UBLAS_TEST_DEF( test_case_1 )
+ * {
+ *   // do your test stuff
+ * }
+ *
+ * BOOST_UBLAS_TEST_DEF( test_case_2 )
+ * {
+ *   // do your test stuff
+ * }
+ *
+ * // ...
+ *
+ * BOOST_UBLAS_TEST_DEF( test_case_n )
+ * {
+ *   // do your test stuff
+ * }
+ *
+ * int main()
+ * {
+ *   BOOST_UBLAS_TEST_SUITE( "My Test Suite" ); // optional
+ *
+ *   BOOST_UBLAS_TEST_BEGIN();
+ *     BOOST_UBLAS_TEST_DO( test_case_1 );
+ *     BOOST_UBLAS_TEST_DO( test_case_2 );
+ *     // ...
+ *     BOOST_UBLAS_TEST_DO( test_case_n );
+ *   BOOST_UBLAS_TEST_END();
+ * }
+ * </pre>
+ * Inside each <em>test_case_<code>k</code></em> you can use the various
+ * \c BOOST_UBLAS_TEST_CHECK* macros.
+ *
+ * <hr/>
+ *
+ *  Copyright (c) 2009-2012, Marco Guazzone
+ *
+ *  Distributed under the Boost Software License, Version 1.0. (See
+ *  accompanying file LICENSE_1_0.txt or copy at
+ *  http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * \author Marco Guazzone, marco.guazzone@gmail.com
+ */
 
-#ifndef TEST_UTILS_HPP
-#define TEST_UTILS_HPP
+#ifndef BOOST_NUMERIC_UBLAS_TEST_UTILS_HPP
+#define BOOST_NUMERIC_UBLAS_TEST_UTILS_HPP
 
 
+#include <boost/numeric/ublas/detail/config.hpp>
+#include <boost/numeric/ublas/traits.hpp>
+#include <cmath>
+#include <complex>
+#include <cstddef>
 #include <iostream>
+#include <limits>
+#include <stdexcept>
+
+namespace boost { namespace numeric { namespace ublas { namespace test { namespace detail { namespace /*<unnamed>*/ {
+
+/// Check if the given complex number is a NaN.
+template <typename T>
+BOOST_UBLAS_INLINE
+bool isnan(::std::complex<T> const& z)
+{
+	// According to IEEE, NaN is different even by itself
+	return (z != z) || ::std::isnan(z.real()) || ::std::isnan(z.imag());
+}
+
+/// Check if two (real) numbers are close each other (wrt a given tolerance).
+template <typename T1, typename T2, typename T3>
+BOOST_UBLAS_INLINE
+bool close_to(T1 x, T2 y, T3 tol)
+{
+	typedef typename promote_traits<typename promote_traits<T1,T2>::promote_type,
+									T3>::promote_type real_type;
+
+    if (::std::isnan(x) || ::std::isnan(y))
+    {
+        // According to IEEE, NaN is different even by itself
+        return false;
+    }
+    return ::std::abs(x-y) <= (::std::max(static_cast<real_type>(::std::abs(x)), static_cast<real_type>(::std::abs(y)))*tol);
+}
+
+/// Check if two complex numbers are close each other (wrt a given tolerance).
+template <typename T1, typename T2, typename T3>
+BOOST_UBLAS_INLINE
+bool close_to(::std::complex<T1> const& x, ::std::complex<T2> const& y, T3 tol)
+{
+	typedef typename promote_traits<typename promote_traits<T1,T2>::promote_type,
+									T3>::promote_type real_type;
+
+    if (isnan(x) || isnan(y))
+    {
+        // According to IEEE, NaN is different even by itself
+        return false;
+    }
+	::std::complex<real_type> xx(x);
+	::std::complex<real_type> yy(y);
+    return ::std::abs(xx-yy) <= (::std::max(::std::abs(xx), ::std::abs(yy))*tol);
+}
+
+/// Check if two (real) numbers are close each other (wrt a given tolerance).
+template <typename T1, typename T2, typename T3>
+BOOST_UBLAS_INLINE
+bool rel_close_to(T1 x, T2 y, T3 tol)
+{
+	typedef typename promote_traits<typename promote_traits<T1,T2>::promote_type,
+									T3>::promote_type real_type;
+
+    if (::std::isnan(x) || ::std::isnan(y))
+    {
+        // According to IEEE, NaN is different even by itself
+        return false;
+    }
+    return ::std::abs(x-y)/::std::abs(y) <= tol;
+}
+
+/// Check if two complex numbers are close each other (wrt a given tolerance).
+template <typename T1, typename T2, typename T3>
+BOOST_UBLAS_INLINE
+bool rel_close_to(::std::complex<T1> const& x, ::std::complex<T2> const& y, T3 tol)
+{
+	typedef typename promote_traits<typename promote_traits<T1,T2>::promote_type,
+									T3>::promote_type real_type;
+
+    if (isnan(x) || isnan(y))
+    {
+        // According to IEEE, NaN is different even by itself
+        return false;
+    }
+	::std::complex<real_type> xx(x);
+	::std::complex<real_type> yy(y);
+    return ::std::abs(xx-yy)/::std::abs(yy) <= tol;
+}
+
+}}}}}} // Namespace boost::numeric::ublas::test::detail::<unnamed>
 
 
-#define EXPAND_(x) x
+/// Expand its argument \a x.
+#define BOOST_UBLAS_TEST_EXPAND_(x) x
 
-#define STRINGIFY_(x) #x
 
-#define JOIN_(x,y) x ## y
+/// Expand its argument \a x inside parenthesis.
+#define BOOST_UBLAS_TEST_EXPANDP_(x) (x)
 
+
+/// Transform its argument \a x into a string.
+#define BOOST_UBLAS_TEST_STRINGIFY_(x) #x
+
+
+/// Concatenate its two \e string arguments \a x and \a y.
+#define BOOST_UBLAS_TEST_JOIN_(x,y) x ## y
+
+
+/// Output the message \a x if in debug-mode; otherwise output nothing.
+/// Note: we don't use macro expansion inside parenthesis to let \a m be an
+///  expression of the form <code>a &lt;&lt; b</code>.
 #ifndef NDEBUG
-#       define BOOST_UBLAS_DEBUG_TRACE(x) std::cerr << "[Debug>> " << EXPAND_(x) << std::endl
+# 	define BOOST_UBLAS_DEBUG_TRACE(x) ::std::cerr << "[Debug>> " << BOOST_UBLAS_TEST_EXPAND_(x) << ::std::endl
 #else
-#       define BOOST_UBLAS_DEBUG_TRACE(x) /**/
+# 	define BOOST_UBLAS_DEBUG_TRACE(x) /**/
 #endif // NDEBUG
 
-#define BOOST_UBLAS_TEST_BEGIN() unsigned int test_fails_(0)
 
-#define BOOST_UBLAS_TEST_DEF(x) void EXPAND_(x)(unsigned int& test_fails_)
-
-#define BOOST_UBLAS_TEST_DO(x) EXPAND_(x)(test_fails_)
-
-#define BOOST_UBLAS_TEST_END() if (test_fails_ > 0) { std::cerr << "Number of failed tests: " << test_fails_ << std::endl; return -1; \
-} else { std::cerr << "No failed test" << std::endl; return 0; }
-
-#define BOOST_UBLAS_TEST_CHECK(x) if (!(x)) { std::cerr << "Failed assertion: " << STRINGIFY_(x) << std::endl; ++test_fails_; }
+/// Define the name \a m of the entire test suite.
+#define BOOST_UBLAS_TEST_SUITE(m) ::std::cerr << "--- Test Suite: " << BOOST_UBLAS_TEST_EXPAND_(m) << " ---" << ::std::endl;
 
 
-#endif // TEST_UTILS_HPP
+/// Define the beginning of a test suite.
+#define BOOST_UBLAS_TEST_BEGIN() 	/* [BOOST_UBLAS_TEST_BEGIN] */ \
+									{ \
+										/* Begin of Test Suite */ \
+										::std::size_t test_fails__(0) \
+									/* [/BOOST_UBLAS_TEST_BEGIN] */
+
+
+/// Define a test case \a x inside the current test suite.
+#define BOOST_UBLAS_TEST_DEF(x) static void BOOST_UBLAS_TEST_EXPAND_(x)(::std::size_t& test_fails__)
+
+
+/// Call the test case \a x.
+#define BOOST_UBLAS_TEST_DO(x) 	/* [BOOST_UBLAS_TEST_DO] */ \
+								try \
+								{ \
+									BOOST_UBLAS_TEST_EXPAND_(x)(test_fails__); \
+								} \
+								catch (::std::exception& e) \
+								{ \
+									++test_fails__; \
+									BOOST_UBLAS_TEST_ERROR( e.what() ); \
+								} \
+								catch (...) \
+								{ \
+									++test_fails__; \
+								} \
+								/* [/BOOST_UBLAS_TEST_DO] */
+
+
+/// Define the end of a test suite.
+#define BOOST_UBLAS_TEST_END() 	/* [BOOST_UBLAS_TEST_END] */ \
+								if (test_fails__ > 0) \
+								{ \
+									::std::cerr << "Number of failed tests: " << test_fails__ << ::std::endl; \
+								} \
+								else \
+								{ \
+									::std::cerr << "No failed test" << ::std::endl; \
+								} \
+								} /* End of test suite */ \
+								/* [/BOOST_UBLAS_TEST_END] */
+
+
+/// Output the message \a m.
+/// Note: we don't use macro expansion inside parenthesis to let \a m be an
+///  expression of the form <code>a &lt;&lt; b</code>.
+#define BOOST_UBLAS_TEST_TRACE(m) ::std::cerr << "[Info>> " << BOOST_UBLAS_TEST_EXPAND_(m) << ::std::endl
+
+
+/// Check the truth of assertion \a x.
+#define BOOST_UBLAS_TEST_CHECK(x)	/* [BOOST_UBLAS_TEST_CHECK] */ \
+									if (!BOOST_UBLAS_TEST_EXPANDP_(x)) \
+									{ \
+										BOOST_UBLAS_TEST_ERROR( "Failed assertion: " << BOOST_UBLAS_TEST_STRINGIFY_(x) ); \
+										++test_fails__; \
+									} \
+									/* [/BOOST_UBLAS_TEST_CHECK] */
+
+
+/// Check for the equality of \a x against \a y.
+#define BOOST_UBLAS_TEST_CHECK_EQ(x,y)	/* [BOOST_UBLAS_TEST_CHECK_EQUAL] */ \
+										if (!(BOOST_UBLAS_TEST_EXPANDP_(x) == BOOST_UBLAS_TEST_EXPANDP_(y))) \
+										{ \
+											BOOST_UBLAS_TEST_ERROR( "Failed assertion: (" << BOOST_UBLAS_TEST_STRINGIFY_(x) << " == " << BOOST_UBLAS_TEST_STRINGIFY_(y) << ")" ); \
+											++test_fails__; \
+										} \
+										/* [/BOOST_UBLAS_TEST_CHECK_EQUAL] */
+
+
+/// Alias for macro \c BOOST_UBLAS_TEST_CHECK_EQ (for backward compatibility).
+#define BOOST_UBLAS_TEST_CHECK_EQUAL(x,y) BOOST_UBLAS_TEST_CHECK_EQ(x,y)
+
+
+/// Check that \a x and \a y are close with respect to a given precision \a e.
+#define BOOST_UBLAS_TEST_CHECK_CLOSE(x,y,e)	/* [BOOST_UBLAS_TEST_CHECK_CLOSE] */ \
+											if (!::boost::numeric::ublas::test::detail::close_to(BOOST_UBLAS_TEST_EXPAND_(x), BOOST_UBLAS_TEST_EXPAND_(y), BOOST_UBLAS_TEST_EXPAND_(e))) \
+											{ \
+												BOOST_UBLAS_TEST_ERROR( "Failed assertion: abs(" << BOOST_UBLAS_TEST_STRINGIFY_(x) << "-" << BOOST_UBLAS_TEST_STRINGIFY_(y) << ") <= " << BOOST_UBLAS_TEST_STRINGIFY_(e) << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y) << " and " << BOOST_UBLAS_TEST_STRINGIFY_(e) << " == " << BOOST_UBLAS_TEST_EXPANDP_(e) << "]" ); \
+												++test_fails__; \
+											} \
+											/* [/BOOST_UBLAS_TEST_CHECK_CLOSE] */
+
+
+/// Alias for macro \c BOOST_UBLAS_TEST_CHECK_CLOSE (for backward compatibility),
+#define BOOST_UBLAS_TEST_CHECK_PRECISION(x,y,e) BOOST_UBLAS_TEST_CHECK_CLOSE(x,y,e)
+
+
+/// Check that \a x is close to \a y with respect to a given relative precision \a e.
+#define BOOST_UBLAS_TEST_CHECK_REL_CLOSE(x,y,e)	/* [BOOST_UBLAS_TEST_CHECK_REL_CLOSE] */ \
+												if (!::boost::numeric::ublas::test::detail::rel_close_to(BOOST_UBLAS_TEST_EXPAND_(x), BOOST_UBLAS_TEST_EXPAND_(y), BOOST_UBLAS_TEST_EXPAND_(e))) \
+												{ \
+													BOOST_UBLAS_TEST_ERROR( "Failed assertion: abs((" << BOOST_UBLAS_TEST_STRINGIFY_(x) << "-" << BOOST_UBLAS_TEST_STRINGIFY_(y) << ")/" << BOOST_UBLAS_TEST_STRINGIFY_(y) << ") <= " << BOOST_UBLAS_TEST_STRINGIFY_(e)  << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y) << " and " << BOOST_UBLAS_TEST_STRINGIFY_(e) << " == " << BOOST_UBLAS_TEST_EXPANDP_(e) << "]" ); \
+													++test_fails__; \
+												} \
+												/* [/BOOST_UBLAS_TEST_CHECK_REL_CLOSE] */
+
+
+/// Alias for macro \c BOOST_UBLAS_TEST_CHECK_REL_CLOSE (for backward compatibility),
+#define BOOST_UBLAS_TEST_CHECK_REL_PRECISION(x,y,e) BOOST_UBLAS_TEST_CHECK_REL_CLOSE(x,y,e)
+
+
+/// Check that elements of \a x and \a y are equal.
+#define BOOST_UBLAS_TEST_CHECK_VECTOR_EQ(x,y,n)	/* [BOOST_UBLAS_TEST_CHECK_VECTOR_EQ] */ \
+												if (BOOST_UBLAS_TEST_EXPANDP_(n) > 0) \
+												{ \
+													::std::size_t n__ = BOOST_UBLAS_TEST_EXPAND_(n); \
+													for (::std::size_t i__ = n__; i__ > 0; --i__) \
+													{ \
+														if (!(BOOST_UBLAS_TEST_EXPANDP_(x)[n__-i__]==BOOST_UBLAS_TEST_EXPANDP_(y)[n__-i__])) \
+														{ \
+															BOOST_UBLAS_TEST_ERROR( "Failed assertion: (" << BOOST_UBLAS_TEST_STRINGIFY_(x[i__]) << "==" << BOOST_UBLAS_TEST_STRINGIFY_(y[i__]) << ")" << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x[i__]) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x)[n__-i__] << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y[i__]) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y)[n__-i__] << ", " << BOOST_UBLAS_TEST_STRINGIFY_(i__) << " == " << i__ << " and " << BOOST_UBLAS_TEST_STRINGIFY_(n) << " == " << n__ << "]" ); \
+															++test_fails__; \
+														} \
+													} \
+												} \
+												/* [/BOOST_UBLAS_TEST_CHECK_VECTOR_EQ] */
+
+
+/// Check that elements of \a x and \a y are close with respect to a given precision \a e.
+#define BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE(x,y,n,e)	/* [BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE] */ \
+														if (BOOST_UBLAS_TEST_EXPANDP_(n) > 0) \
+														{ \
+															::std::size_t n__ = BOOST_UBLAS_TEST_EXPAND_(n); \
+															for (::std::size_t i__ = n__; i__ > 0; --i__) \
+															{ \
+																if (!::boost::numeric::ublas::test::detail::close_to(BOOST_UBLAS_TEST_EXPANDP_(x)[n__-i__], BOOST_UBLAS_TEST_EXPANDP_(y)[n__-i__], BOOST_UBLAS_TEST_EXPANDP_(e))) \
+																{ \
+																	BOOST_UBLAS_TEST_ERROR( "Failed assertion: abs((" << BOOST_UBLAS_TEST_STRINGIFY_(x[i__]) << "-" << BOOST_UBLAS_TEST_STRINGIFY_(y[i__]) << ") <= " << BOOST_UBLAS_TEST_STRINGIFY_(e)  << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x[i__]) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x)[n__-i__] << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y[i__]) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y)[n__-i__] << ", " << BOOST_UBLAS_TEST_STRINGIFY_(i__) << " == " << i__ << " and " << BOOST_UBLAS_TEST_STRINGIFY_(n) << " == " << n__ << "]" ); \
+																	++test_fails__; \
+																} \
+															} \
+														} \
+														/* [/BOOST_UBLAS_TEST_CHECK_VECTOR_CLOSE] */
+
+
+/// Check that elements of \a x and \a y are close with respect to a given relative precision \a e.
+#define BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE(x,y,n,e)	/* [BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE] */ \
+														if (BOOST_UBLAS_TEST_EXPANDP_(n) > 0) \
+														{ \
+															::std::size_t n__ = BOOST_UBLAS_TEST_EXPAND_(n); \
+															for (::std::size_t i__ = n__; i__ > 0; --i__) \
+															{ \
+																if (!::boost::numeric::ublas::test::detail::rel_close_to(BOOST_UBLAS_TEST_EXPANDP_(x)[n__-i__], BOOST_UBLAS_TEST_EXPANDP_(y)[n__-i__], BOOST_UBLAS_TEST_EXPANDP_(e))) \
+																{ \
+																	BOOST_UBLAS_TEST_ERROR( "Failed assertion: abs((" << BOOST_UBLAS_TEST_STRINGIFY_(x[i__]) << "-" << BOOST_UBLAS_TEST_STRINGIFY_(y[i__]) << ") <= " << BOOST_UBLAS_TEST_STRINGIFY_(e)  << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x[i__]) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x)[n__-i__] << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y[i__]) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y)[n__-i__] << ", " << BOOST_UBLAS_TEST_STRINGIFY_(i__) << " == " << i__ << " and " << BOOST_UBLAS_TEST_STRINGIFY_(n) << " == " << n__ << "]" ); \
+																	++test_fails__; \
+																} \
+															} \
+														} \
+														/* [/BOOST_UBLAS_TEST_CHECK_VECTOR_REL_CLOSE] */
+
+
+/// Check that elements of matrices \a x and \a y are equal.
+#define BOOST_UBLAS_TEST_CHECK_MATRIX_EQ(x,y,nr,nc)	/* [BOOST_UBLAS_TEST_CHECK_MATRIX_EQ] */ \
+													for (::std::size_t i__ = 0; i__ < BOOST_UBLAS_TEST_EXPANDP_(nr); ++i__) \
+													{ \
+														for (::std::size_t j__ = 0; j__ < BOOST_UBLAS_TEST_EXPANDP_(nc); ++j__) \
+														{ \
+															if (!(BOOST_UBLAS_TEST_EXPANDP_(x)(i__,j__)==BOOST_UBLAS_TEST_EXPANDP_(y)(i__,j__))) \
+															{ \
+																BOOST_UBLAS_TEST_ERROR( "Failed assertion: (" << BOOST_UBLAS_TEST_STRINGIFY_(x(i__,j__)) << " == " << BOOST_UBLAS_TEST_STRINGIFY_(y(i__,j__)) << ") [with " << BOOST_UBLAS_TEST_STRINGIFY_(x(i__,j__)) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x)(i__,j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y(i__,j__)) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y)(i__,j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(i__) << " == " << i__ << ", " << BOOST_UBLAS_TEST_STRINGIFY_(j__) << " == " << BOOST_UBLAS_TEST_EXPANDP_(j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(nr) << " == " << BOOST_UBLAS_TEST_EXPANDP_(nr) << " and " << BOOST_UBLAS_TEST_STRINGIFY_(nc) << " == " << BOOST_UBLAS_TEST_EXPANDP_(nc) << "]" ); \
+																++test_fails__; \
+															} \
+														} \
+													} \
+													/* [/BOOST_UBLAS_TEST_CHECK_MATRIX_EQ] */
+
+
+/// Check that elements of matrices \a x and \a y are close with respect to a given precision \a e.
+#define BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE(x,y,nr,nc,e)	/* [BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE] */ \
+															for (::std::size_t i__ = 0; i__ < BOOST_UBLAS_TEST_EXPANDP_(nr); ++i__) \
+															{ \
+																for (::std::size_t j__ = 0; j__ < BOOST_UBLAS_TEST_EXPANDP_(nc); ++j__) \
+																{ \
+																	if (!::boost::numeric::ublas::test::detail::close_to(BOOST_UBLAS_TEST_EXPANDP_(x)(i__,j__), BOOST_UBLAS_TEST_EXPANDP_(y)(i__,j__), BOOST_UBLAS_TEST_EXPANDP_(e))) \
+																	{ \
+																		BOOST_UBLAS_TEST_ERROR( "Failed assertion: abs((" << BOOST_UBLAS_TEST_STRINGIFY_(x(i__,j__)) << "-" << BOOST_UBLAS_TEST_STRINGIFY_(y(i__,j__)) << ") <= " << BOOST_UBLAS_TEST_STRINGIFY_(e)  << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x(i__,j__)) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x)(i__,j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y(i__,j__)) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y)(i__,j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(i__) << " == " << i__ << ", " << BOOST_UBLAS_TEST_STRINGIFY_(j__) << " == " << BOOST_UBLAS_TEST_EXPANDP_(j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(nr) << " == " << BOOST_UBLAS_TEST_EXPANDP_(nr) << " and " << BOOST_UBLAS_TEST_STRINGIFY_(nc) << " == " << BOOST_UBLAS_TEST_EXPANDP_(nc) << "]" ); \
+																		++test_fails__; \
+																	} \
+																} \
+															} \
+															/* [/BOOST_UBLAS_TEST_CHECK_MATRIX_CLOSE] */
+
+
+/// Check that elements of matrices \a x and \a y are close with respect to a given relative precision \a e.
+#define BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE(x,y,nr,nc,e)	/* [BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE] */ \
+																for (::std::size_t i__ = 0; i__ < BOOST_UBLAS_TEST_EXPANDP_(nr); ++i__) \
+																{ \
+																	for (::std::size_t j__ = 0; j__ < BOOST_UBLAS_TEST_EXPANDP_(nc); ++j__) \
+																	{ \
+																		if (!::boost::numeric::ublas::test::detail::rel_close_to(BOOST_UBLAS_TEST_EXPANDP_(x)(i__,j__), BOOST_UBLAS_TEST_EXPANDP_(y)(i__,j__), BOOST_UBLAS_TEST_EXPANDP_(e))) \
+																		{ \
+																			BOOST_UBLAS_TEST_ERROR( "Failed assertion: abs((" << BOOST_UBLAS_TEST_STRINGIFY_(x(i__,j__)) << "-" << BOOST_UBLAS_TEST_STRINGIFY_(y(i__,j__)) << ") <= " << BOOST_UBLAS_TEST_STRINGIFY_(e)  << " [with " << BOOST_UBLAS_TEST_STRINGIFY_(x(i__,j__)) << " == " << BOOST_UBLAS_TEST_EXPANDP_(x)(i__,j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(y(i__,j__)) << " == " << BOOST_UBLAS_TEST_EXPANDP_(y)(i__,j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(i__) << " == " << i__ << ", " << BOOST_UBLAS_TEST_STRINGIFY_(j__) << " == " << BOOST_UBLAS_TEST_EXPANDP_(j__) << ", " << BOOST_UBLAS_TEST_STRINGIFY_(nr) << " == " << BOOST_UBLAS_TEST_EXPANDP_(nr) << " and " << BOOST_UBLAS_TEST_STRINGIFY_(nc) << " == " << BOOST_UBLAS_TEST_EXPANDP_(nc) << "]" ); \
+																			++test_fails__; \
+																		} \
+																	} \
+																} \
+																/* [/BOOST_UBLAS_TEST_CHECK_MATRIX_REL_CLOSE] */
+
+///< Output the error message \a x.
+#ifdef _MSC_VER
+# define BOOST_UBLAS_TEST_ERROR(x) ::std::cerr << "[Error (" << __FILE__ << ":" << __FUNCTION__ << ":" << __LINE__ << ")>> " << BOOST_UBLAS_TEST_EXPAND_(x) << ::std::endl
+#else
+# define BOOST_UBLAS_TEST_ERROR(x) ::std::cerr << "[Error (" << __FILE__ << ":" << __func__ << ":" << __LINE__ << ")>> " << BOOST_UBLAS_TEST_EXPAND_(x) << ::std::endl
+#endif
+
+#endif // BOOST_NUMERIC_UBLAS_TEST_UTILS_HPP