diff --git a/include/boost/numeric/ublas/functional.hpp b/include/boost/numeric/ublas/functional.hpp
index 7b071272..08757fe0 100644
--- a/include/boost/numeric/ublas/functional.hpp
+++ b/include/boost/numeric/ublas/functional.hpp
@@ -222,7 +222,7 @@ namespace boost { namespace numeric { namespace ublas {
 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
         static const bool computed ;
 #else
-      static const bool computed = true ;
+        static const bool computed = true ;
 #endif
 
         static BOOST_UBLAS_INLINE
@@ -331,35 +331,32 @@ namespace boost { namespace numeric { namespace ublas {
     // Vector functors
 
     // Unary returning scalar
-    template<class T>
+    template<class V>
     struct vector_scalar_unary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef T value_type;
-        typedef T result_type;
+        typedef typename V::value_type value_type;
+        typedef typename V::value_type result_type;
     };
 
-    template<class T>
+    template<class V>
     struct vector_sum: 
-        public vector_scalar_unary_functor<T> {
-        typedef typename vector_scalar_unary_functor<T>::size_type size_type;
-        typedef typename vector_scalar_unary_functor<T>::difference_type difference_type;
-        typedef typename vector_scalar_unary_functor<T>::value_type value_type;
-        typedef typename vector_scalar_unary_functor<T>::result_type result_type;
+        public vector_scalar_unary_functor<V> {
+        typedef typename vector_scalar_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_unary_functor<V>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const vector_expression<E> &e) { 
             result_type t = result_type (0);
-            size_type size (e ().size ());
-            for (size_type i = 0; i < size; ++ i)
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i)
                 t += e () (i);
             return t;
         }
         // Dense case
-        template<class I>
+        template<class D, class I>
         static BOOST_UBLAS_INLINE
-        result_type apply (difference_type size, I it) { 
+        result_type apply (D size, I it) { 
             result_type t = result_type (0);
             while (-- size >= 0)
                 t += *it, ++ it;
@@ -377,39 +374,36 @@ namespace boost { namespace numeric { namespace ublas {
     };
 
     // Unary returning real scalar 
-    template<class T>
+    template<class V>
     struct vector_scalar_real_unary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef T value_type;
-        typedef typename type_traits<T>::real_type real_type;
+        typedef typename V::value_type value_type;
+        typedef typename type_traits<value_type>::real_type real_type;
         typedef real_type result_type;
     };
 
-    template<class T>
+    template<class V>
     struct vector_norm_1:
-        public vector_scalar_real_unary_functor<T> {
-        typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
-        typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
-        typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
-        typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
-        typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
+        public vector_scalar_real_unary_functor<V> {
+        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
+        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const vector_expression<E> &e) {
             real_type t = real_type ();
-            size_type size (e ().size ());
-            for (size_type i = 0; i < size; ++ i) {
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) {
                 real_type u (type_traits<value_type>::type_abs (e () (i)));
                 t += u;
             }
             return t;
         }
         // Dense case
-        template<class I>
+        template<class D, class I>
         static BOOST_UBLAS_INLINE
-        result_type apply (difference_type size, I it) {
+        result_type apply (D size, I it) {
             real_type t = real_type ();
             while (-- size >= 0) {
                 real_type u (type_traits<value_type>::norm_1 (*it));
@@ -431,22 +425,21 @@ namespace boost { namespace numeric { namespace ublas {
             return t;
         }
     };
-    template<class T>
+    template<class V>
     struct vector_norm_2:
-        public vector_scalar_real_unary_functor<T> {
-        typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
-        typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
-        typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
-        typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
-        typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
+        public vector_scalar_real_unary_functor<V> {
+        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
+        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const vector_expression<E> &e) {
 #ifndef BOOST_UBLAS_SCALED_NORM
             real_type t = real_type ();
-            size_type size (e ().size ());
-            for (size_type i = 0; i < size; ++ i) {
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) {
                 real_type u (type_traits<value_type>::norm_2 (e () (i)));
                 t +=  u * u;
             }
@@ -470,9 +463,9 @@ namespace boost { namespace numeric { namespace ublas {
 #endif
         }
         // Dense case
-        template<class I>
+        template<class D, class I>
         static BOOST_UBLAS_INLINE
-        result_type apply (difference_type size, I it) {
+        result_type apply (D size, I it) {
 #ifndef BOOST_UBLAS_SCALED_NORM
             real_type t = real_type ();
             while (-- size >= 0) {
@@ -530,21 +523,20 @@ namespace boost { namespace numeric { namespace ublas {
 #endif
         }
     };
-    template<class T>
+    template<class V>
     struct vector_norm_inf:
-        public vector_scalar_real_unary_functor<T> {
-        typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
-        typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
-        typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
-        typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
-        typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
+        public vector_scalar_real_unary_functor<V> {
+        typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
+        typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const vector_expression<E> &e) {
             real_type t = real_type ();
-            size_type size (e ().size ());
-            for (size_type i = 0; i < size; ++ i) {
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) {
                 real_type u (type_traits<value_type>::norm_inf (e () (i)));
                 if (u > t)
                     t = u;
@@ -552,9 +544,9 @@ namespace boost { namespace numeric { namespace ublas {
             return t;
         }
         // Dense case
-        template<class I>
+        template<class D, class I>
         static BOOST_UBLAS_INLINE
-        result_type apply (difference_type size, I it) {
+        result_type apply (D size, I it) {
             real_type t = real_type ();
             while (-- size >= 0) {
                 real_type u (type_traits<value_type>::norm_inf (*it));
@@ -580,23 +572,19 @@ namespace boost { namespace numeric { namespace ublas {
     };
 
     // Unary returning index
-    template<class T>
+    template<class V>
     struct vector_scalar_index_unary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef T value_type;
-        typedef typename type_traits<T>::real_type real_type;
-        typedef size_type result_type;
+        typedef typename V::value_type value_type;
+        typedef typename type_traits<value_type>::real_type real_type;
+        typedef typename V::size_type result_type;
     };
 
-    template<class T>
+    template<class V>
     struct vector_index_norm_inf:
-        public vector_scalar_index_unary_functor<T> {
-        typedef typename vector_scalar_index_unary_functor<T>::size_type size_type;
-        typedef typename vector_scalar_index_unary_functor<T>::difference_type difference_type;
-        typedef typename vector_scalar_index_unary_functor<T>::value_type value_type;
-        typedef typename vector_scalar_index_unary_functor<T>::real_type real_type;
-        typedef typename vector_scalar_index_unary_functor<T>::result_type result_type;
+        public vector_scalar_index_unary_functor<V> {
+        typedef typename vector_scalar_index_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_index_unary_functor<V>::real_type real_type;
+        typedef typename vector_scalar_index_unary_functor<V>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
@@ -604,8 +592,9 @@ namespace boost { namespace numeric { namespace ublas {
             // ISSUE For CBLAS compatibility return 0 index in empty case
             result_type i_norm_inf (0);
             real_type t = real_type ();
-            size_type size (e ().size ());
-            for (size_type i = 0; i < size; ++ i) {
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) {
                 real_type u (type_traits<value_type>::norm_inf (e () (i)));
                 if (u > t) {
                     i_norm_inf = i;
@@ -615,9 +604,9 @@ namespace boost { namespace numeric { namespace ublas {
             return i_norm_inf;
         }
         // Dense case
-        template<class I>
+        template<class D, class I>
         static BOOST_UBLAS_INLINE
-        result_type apply (difference_type size, I it) {
+        result_type apply (D size, I it) {
             // ISSUE For CBLAS compatibility return 0 index in empty case
             result_type i_norm_inf (0);
             real_type t = real_type ();
@@ -651,21 +640,17 @@ namespace boost { namespace numeric { namespace ublas {
     };
 
     // Binary returning scalar
-    template<class T1, class T2, class TR>
+    template<class V1, class V2, class TV>
     struct vector_scalar_binary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef TR value_type;
-        typedef TR result_type;
+        typedef TV value_type;
+        typedef TV result_type;
     };
 
-    template<class T1, class T2, class TR>
+    template<class V1, class V2, class TV>
     struct vector_inner_prod:
-        public vector_scalar_binary_functor<T1, T2, TR> {
-        typedef typename vector_scalar_binary_functor<T1, T2, TR>::size_type size_type ;
-        typedef typename vector_scalar_binary_functor<T1, T2, TR>::difference_type difference_type;
-        typedef typename vector_scalar_binary_functor<T1, T2, TR>::value_type value_type;
-        typedef typename vector_scalar_binary_functor<T1, T2, TR>::result_type result_type;
+        public vector_scalar_binary_functor<V1, V2, TV> {
+        typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type;
+        typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type;
 
         template<class C1, class C2>
         static BOOST_UBLAS_INLINE
@@ -673,23 +658,24 @@ namespace boost { namespace numeric { namespace ublas {
                            const vector_container<C2> &c2) {
 #ifdef BOOST_UBLAS_USE_SIMD
             using namespace raw;
-            size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
-            const T1 *data1 = data_const (c1 ());
-            const T2 *data2 = data_const (c2 ());
-            size_type s1 = stride (c1 ());
-            size_type s2 = stride (c2 ());
+            typedef typename C1::size_type vector_size_type;
+            vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
+            const typename V1::value_type *data1 = data_const (c1 ());
+            const typename V1::value_type *data2 = data_const (c2 ());
+            vector_size_type s1 = stride (c1 ());
+            vector_size_type s2 = stride (c2 ());
             result_type t = result_type (0);
             if (s1 == 1 && s2 == 1) {
-                for (size_type i = 0; i < size; ++ i)
+                for (vector_size_type i = 0; i < size; ++ i)
                     t += data1 [i] * data2 [i];
             } else if (s2 == 1) {
-                for (size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
+                for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
                     t += data1 [i1] * data2 [i];
             } else if (s1 == 1) {
-                for (size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
+                for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
                     t += data1 [i] * data2 [i2];
             } else {
-                for (size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
+                for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
                     t += data1 [i1] * data2 [i2];
             }
             return t;
@@ -703,21 +689,22 @@ namespace boost { namespace numeric { namespace ublas {
         static BOOST_UBLAS_INLINE
         result_type apply (const vector_expression<E1> &e1,
                            const vector_expression<E2> &e2) {
-            size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
+            typedef typename E1::size_type vector_size_type;
+            vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
             result_type t = result_type (0);
 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
-            for (size_type i = 0; i < size; ++ i)
+            for (vector_size_type i = 0; i < size; ++ i)
                 t += e1 () (i) * e2 () (i);
 #else
-            size_type i (0);
+            vector_size_type i (0);
             DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
 #endif
             return t;
         }
         // Dense case
-        template<class I1, class I2>
+        template<class D, class I1, class I2>
         static BOOST_UBLAS_INLINE
-        result_type apply (difference_type size, I1 it1, I2 it2) {
+        result_type apply (D size, I1 it1, I2 it2) {
             result_type t = result_type (0);
 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
             while (-- size >= 0)
@@ -732,13 +719,14 @@ namespace boost { namespace numeric { namespace ublas {
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
             result_type t = result_type (0);
-            difference_type it1_size (it1_end - it1);
-            difference_type it2_size (it2_end - it2);
-            difference_type diff (0);
+            typedef typename I1::difference_type vector_difference_type;
+            vector_difference_type it1_size (it1_end - it1);
+            vector_difference_type it2_size (it2_end - it2);
+            vector_difference_type diff (0);
             if (it1_size > 0 && it2_size > 0)
                 diff = it2.index () - it1.index ();
             if (diff != 0) {
-                difference_type size = (std::min) (diff, it1_size);
+                vector_difference_type size = (std::min) (diff, it1_size);
                 if (size > 0) {
                     it1 += size;
                     it1_size -= size;
@@ -751,7 +739,7 @@ namespace boost { namespace numeric { namespace ublas {
                     diff += size;
                 }
             }
-            difference_type size ((std::min) (it1_size, it2_size));
+            vector_difference_type size ((std::min) (it1_size, it2_size));
             while (-- size >= 0)
                 t += *it1 * *it2, ++ it1, ++ it2;
             return t;
@@ -762,27 +750,18 @@ namespace boost { namespace numeric { namespace ublas {
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
             result_type t = result_type (0);
             if (it1 != it1_end && it2 != it2_end) {
-                size_type it1_index = it1.index (), it2_index = it2.index ();
                 while (true) {
-                    difference_type compare = it1_index - it2_index;
-                    if (compare == 0) {
+                    if (it1.index () == it2.index ()) {
                         t += *it1 * *it2, ++ it1, ++ it2;
-                        if (it1 != it1_end && it2 != it2_end) {
-                            it1_index = it1.index ();
-                            it2_index = it2.index ();
-                        } else
+                        if (it1 == it1_end || it2 == it2_end)
                             break;
-                    } else if (compare < 0) {
-                        increment (it1, it1_end, - compare);
-                        if (it1 != it1_end)
-                            it1_index = it1.index ();
-                        else
+                    } else if (it1.index () < it2.index ()) {
+                        increment (it1, it1_end, it2.index () - it1.index ());
+                        if (it1 == it1_end)
                             break;
-                    } else if (compare > 0) {
-                        increment (it2, it2_end, compare);
-                        if (it2 != it2_end)
-                            it2_index = it2.index ();
-                        else
+                    } else if (it1.index () > it2.index ()) {
+                        increment (it2, it2_end, it1.index () - it2.index ());
+                        if (it2 == it2_end)
                             break;
                     }
                 }
@@ -794,21 +773,21 @@ namespace boost { namespace numeric { namespace ublas {
     // Matrix functors
 
     // Binary returning vector
-    template<class T1, class T2, class TR>
+    template<class M1, class M2, class TV>
     struct matrix_vector_binary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef TR value_type;
-        typedef TR result_type;
+        typedef typename M1::size_type size_type;
+        typedef typename M1::difference_type difference_type;
+        typedef TV value_type;
+        typedef TV result_type;
     };
 
-    template<class T1, class T2, class TR>
+    template<class M1, class M2, class TV>
     struct matrix_vector_prod1:
-        public matrix_vector_binary_functor<T1, T2, TR> {
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;
+        public matrix_vector_binary_functor<M1, M2, TV> {
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
 
         template<class C1, class C2>
         static BOOST_UBLAS_INLINE
@@ -818,8 +797,8 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_USE_SIMD
             using namespace raw;
             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
-            const T1 *data1 = data_const (c1 ()) + i * stride1 (c1 ());
-            const T2 *data2 = data_const (c2 ());
+            const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
+            const typename M2::value_type *data2 = data_const (c2 ());
             size_type s1 = stride2 (c1 ());
             size_type s2 = stride (c2 ());
             result_type t = result_type (0);
@@ -846,7 +825,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class E1, class E2>
         static BOOST_UBLAS_INLINE
         result_type apply (const matrix_expression<E1> &e1,
-                                 const vector_expression<E2> &e2,
+                           const vector_expression<E2> &e2,
                            size_type i) {
             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
             result_type t = result_type (0);
@@ -905,7 +884,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
-                                 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
+                           sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
             result_type t = result_type (0);
             if (it1 != it1_end && it2 != it2_end) {
                 size_type it1_index = it1.index2 (), it2_index = it2.index ();
@@ -939,7 +918,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
-                                 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
+                           sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
             result_type t = result_type (0);
             while (it1 != it1_end) {
                 t += *it1 * it2 () (it1.index2 ());
@@ -951,7 +930,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
-                                 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
+                           packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
             result_type t = result_type (0);
             while (it2 != it2_end) {
                 t += it1 () (it1.index1 (), it2.index ()) * *it2;
@@ -963,20 +942,20 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
-                                 sparse_bidirectional_iterator_tag) {
+                           sparse_bidirectional_iterator_tag) {
             typedef typename I1::iterator_category iterator1_category;
             typedef typename I2::iterator_category iterator2_category;
             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
         }
     };
 
-    template<class T1, class T2, class TR>
+    template<class M1, class M2, class TV>
     struct matrix_vector_prod2:
-        public matrix_vector_binary_functor<T1, T2, TR> {
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
-        typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;
+        public matrix_vector_binary_functor<M1, M2, TV> {
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
+        typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
 
         template<class C1, class C2>
         static BOOST_UBLAS_INLINE
@@ -986,8 +965,8 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_USE_SIMD
             using namespace raw;
             size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
-            const T1 *data1 = data_const (c1 ());
-            const T2 *data2 = data_const (c2 ()) + i * stride2 (c2 ());
+            const typename M1::value_type *data1 = data_const (c1 ());
+            const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ());
             size_type s1 = stride (c1 ());
             size_type s2 = stride1 (c2 ());
             result_type t = result_type (0);
@@ -1014,7 +993,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class E1, class E2>
         static BOOST_UBLAS_INLINE
         result_type apply (const vector_expression<E1> &e1,
-                                 const matrix_expression<E2> &e2,
+                           const matrix_expression<E2> &e2,
                            size_type i) {
             size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
             result_type t = result_type (0);
@@ -1073,7 +1052,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
-                                 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
+                           sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
             result_type t = result_type (0);
             if (it1 != it1_end && it2 != it2_end) {
                 size_type it1_index = it1.index (), it2_index = it2.index1 ();
@@ -1107,7 +1086,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
-                                 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
+                           packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
             result_type t = result_type (0);
             while (it2 != it2_end) {
                 t += it1 () (it2.index1 ()) * *it2;
@@ -1119,7 +1098,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
-                                 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
+                           sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
             result_type t = result_type (0);
             while (it1 != it1_end) {
                 t += *it1 * it2 () (it1.index (), it2.index2 ());
@@ -1131,7 +1110,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class I1, class I2>
         static BOOST_UBLAS_INLINE
         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
-                                 sparse_bidirectional_iterator_tag) {
+                           sparse_bidirectional_iterator_tag) {
             typedef typename I1::iterator_category iterator1_category;
             typedef typename I2::iterator_category iterator2_category;
             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
@@ -1139,21 +1118,21 @@ namespace boost { namespace numeric { namespace ublas {
     };
 
     // Binary returning matrix
-    template<class T1, class T2, class TR>
+    template<class M1, class M2, class TV>
     struct matrix_matrix_binary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef TR value_type;
-        typedef TR result_type;
+        typedef typename M1::size_type size_type;
+        typedef typename M1::difference_type difference_type;
+        typedef TV value_type;
+        typedef TV result_type;
     };
 
-    template<class T1, class T2, class TR>
+    template<class M1, class M2, class TV>
     struct matrix_matrix_prod:
-        public matrix_matrix_binary_functor<T1, T2, TR> {
-        typedef typename matrix_matrix_binary_functor<T1, T2, TR>::size_type size_type;
-        typedef typename matrix_matrix_binary_functor<T1, T2, TR>::difference_type difference_type;
-        typedef typename matrix_matrix_binary_functor<T1, T2, TR>::value_type value_type;
-        typedef typename matrix_matrix_binary_functor<T1, T2, TR>::result_type result_type;
+        public matrix_matrix_binary_functor<M1, M2, TV> {
+        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type;
+        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type;
+        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type;
+        typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type;
 
         template<class C1, class C2>
         static BOOST_UBLAS_INLINE
@@ -1163,8 +1142,8 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_USE_SIMD
             using namespace raw;
             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
-            const T1 *data1 = data_const (c1 ()) + i * stride1 (c1 ());
-            const T2 *data2 = data_const (c2 ()) + j * stride2 (c2 ());
+            const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
+            const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ());
             size_type s1 = stride2 (c1 ());
             size_type s2 = stride1 (c2 ());
             result_type t = result_type (0);
@@ -1191,7 +1170,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class E1, class E2>
         static BOOST_UBLAS_INLINE
         result_type apply (const matrix_expression<E1> &e1,
-                                 const matrix_expression<E2> &e2,
+                           const matrix_expression<E2> &e2,
                            size_type i, size_type j) {
             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
             result_type t = result_type (0);
@@ -1254,7 +1233,7 @@ namespace boost { namespace numeric { namespace ublas {
             if (it1 != it1_end && it2 != it2_end) {
                 size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
                 while (true) {
-                    difference_type compare = it1_index - it2_index;
+                    difference_type compare = difference_type (it1_index - it2_index);
                     if (compare == 0) {
                         t += *it1 * *it2, ++ it1, ++ it2;
                         if (it1 != it1_end && it2 != it2_end) {
@@ -1282,33 +1261,30 @@ namespace boost { namespace numeric { namespace ublas {
     };
 
     // Unary returning scalar norm
-    template<class T>
+    template<class M>
     struct matrix_scalar_real_unary_functor {
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
-        typedef T value_type;
-        typedef typename type_traits<T>::real_type real_type;
+        typedef typename M::value_type value_type;
+        typedef typename type_traits<value_type>::real_type real_type;
         typedef real_type result_type;
     };
 
-    template<class T>
+    template<class M>
     struct matrix_norm_1:
-        public matrix_scalar_real_unary_functor<T> {
-        typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
+        public matrix_scalar_real_unary_functor<M> {
+        typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
+        typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
+        typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const matrix_expression<E> &e) {
             real_type t = real_type ();
-            size_type size2 (e ().size2 ());
-            for (size_type j = 0; j < size2; ++ j) {
+            typedef typename E::size_type matrix_size_type;
+            matrix_size_type size2 (e ().size2 ());
+            for (matrix_size_type j = 0; j < size2; ++ j) {
                 real_type u = real_type ();
-                size_type size1 (e ().size1 ());
-                for (size_type i = 0; i < size1; ++ i) {
+                matrix_size_type size1 (e ().size1 ());
+                for (matrix_size_type i = 0; i < size1; ++ i) {
                     real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
                     u += v;
                 }
@@ -1318,23 +1294,23 @@ namespace boost { namespace numeric { namespace ublas {
             return t; 
         }
     };
-    template<class T>
+
+    template<class M>
     struct matrix_norm_frobenius:
-        public matrix_scalar_real_unary_functor<T> {
-        typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
+        public matrix_scalar_real_unary_functor<M> {
+        typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
+        typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
+        typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const matrix_expression<E> &e) { 
             real_type t = real_type ();
-            size_type size1 (e ().size1 ());
-            for (size_type i = 0; i < size1; ++ i) {
-                size_type size2 (e ().size2 ());
-                for (size_type j = 0; j < size2; ++ j) {
+            typedef typename E::size_type matrix_size_type;
+            matrix_size_type size1 (e ().size1 ());
+            for (matrix_size_type i = 0; i < size1; ++ i) {
+                matrix_size_type size2 (e ().size2 ());
+                for (matrix_size_type j = 0; j < size2; ++ j) {
                     real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
                     t +=  u * u;
                 }
@@ -1342,24 +1318,24 @@ namespace boost { namespace numeric { namespace ublas {
             return type_traits<real_type>::type_sqrt (t); 
         }
     };
-    template<class T>
+
+    template<class M>
     struct matrix_norm_inf: 
-        public matrix_scalar_real_unary_functor<T> {
-        typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
-        typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
+        public matrix_scalar_real_unary_functor<M> {
+        typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
+        typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
+        typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
 
         template<class E>
         static BOOST_UBLAS_INLINE
         result_type apply (const matrix_expression<E> &e) {
             real_type t = real_type ();
-            size_type size1 (e ().size1 ());
-            for (size_type i = 0; i < size1; ++ i) {
+            typedef typename E::size_type matrix_size_type;
+            matrix_size_type size1 (e ().size1 ());
+            for (matrix_size_type i = 0; i < size1; ++ i) {
                 real_type u = real_type ();
-                size_type size2 (e ().size2 ());
-                for (size_type j = 0; j < size2; ++ j) {
+                matrix_size_type size2 (e ().size2 ());
+                for (matrix_size_type j = 0; j < size2; ++ j) {
                     real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
                     u += v;
                 }
diff --git a/include/boost/numeric/ublas/fwd.hpp b/include/boost/numeric/ublas/fwd.hpp
index 57f912c6..d7a9eb22 100644
--- a/include/boost/numeric/ublas/fwd.hpp
+++ b/include/boost/numeric/ublas/fwd.hpp
@@ -93,11 +93,11 @@ namespace boost { namespace numeric { namespace ublas {
     template<class T, std::size_t N>
     class bounded_vector;
 
-    template<class T = int>
+    template<class T = int, class ALLOC = std::allocator<T> >
     class unit_vector;
-    template<class T = int>
+    template<class T = int, class ALLOC = std::allocator<T> >
     class zero_vector;
-    template<class T = int>
+    template<class T = int, class ALLOC = std::allocator<T> >
     class scalar_vector;
 
     template<class T, std::size_t N>
@@ -130,11 +130,11 @@ namespace boost { namespace numeric { namespace ublas {
     template<class T, std::size_t M, std::size_t N, class L = row_major>
     class bounded_matrix;
 
-    template<class T = int>
+    template<class T = int, class ALLOC = std::allocator<T> >
     class identity_matrix;
-    template<class T = int>
+    template<class T = int, class ALLOC = std::allocator<T> >
     class zero_matrix;
-    template<class T = int>
+    template<class T = int, class ALLOC = std::allocator<T> >
     class scalar_matrix;
 
     template<class T, std::size_t M, std::size_t N>
diff --git a/include/boost/numeric/ublas/matrix.hpp b/include/boost/numeric/ublas/matrix.hpp
index 6e9d2e8f..4112feab 100644
--- a/include/boost/numeric/ublas/matrix.hpp
+++ b/include/boost/numeric/ublas/matrix.hpp
@@ -1968,18 +1968,18 @@ namespace boost { namespace numeric { namespace ublas {
 
 
     // Zero matrix class
-    template<class T>
+    template<class T, class ALLOC>
     class zero_matrix:
-        public matrix_container<zero_matrix<T> > {
+        public matrix_container<zero_matrix<T, ALLOC> > {
 
         typedef const T *const_pointer;
-        typedef zero_matrix<T> self_type;
+        typedef zero_matrix<T, ALLOC> self_type;
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using matrix_container<self_type>::operator ();
 #endif
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
+        typedef typename ALLOC::size_type size_type;
+        typedef typename ALLOC::difference_type difference_type;
         typedef T value_type;
         typedef const T &const_reference;
         typedef T &reference;
@@ -2323,23 +2323,23 @@ namespace boost { namespace numeric { namespace ublas {
         static const value_type zero_;
     };
 
-    template<class T>
-    const typename zero_matrix<T>::value_type zero_matrix<T>::zero_ (0);
+    template<class T, class ALLOC>
+    const typename zero_matrix<T, ALLOC>::value_type zero_matrix<T, ALLOC>::zero_ (0);
 
 
     // Identity matrix class
-    template<class T>
+    template<class T, class ALLOC>
     class identity_matrix:
-        public matrix_container<identity_matrix<T> > {
+        public matrix_container<identity_matrix<T, ALLOC> > {
 
         typedef const T *const_pointer;
-        typedef identity_matrix<T> self_type;
+        typedef identity_matrix<T, ALLOC> self_type;
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using matrix_container<self_type>::operator ();
 #endif
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
+        typedef typename ALLOC::size_type size_type;
+        typedef typename ALLOC::difference_type difference_type;
         typedef T value_type;
         typedef const T &const_reference;
         typedef T &reference;
@@ -2705,19 +2705,19 @@ namespace boost { namespace numeric { namespace ublas {
         static const value_type one_;
     };
 
-    template<class T>
-    const typename identity_matrix<T>::value_type identity_matrix<T>::zero_ (0);
-    template<class T>
-    const typename identity_matrix<T>::value_type identity_matrix<T>::one_ (1);
+    template<class T, class ALLOC>
+    const typename identity_matrix<T, ALLOC>::value_type identity_matrix<T, ALLOC>::zero_ (0);
+    template<class T, class ALLOC>
+    const typename identity_matrix<T, ALLOC>::value_type identity_matrix<T, ALLOC>::one_ (1);
 
 
     // Scalar matrix class
-    template<class T>
+    template<class T, class ALLOC>
     class scalar_matrix:
-        public matrix_container<scalar_matrix<T> > {
+        public matrix_container<scalar_matrix<T, ALLOC> > {
 
         typedef const T *const_pointer;
-        typedef scalar_matrix<T> self_type;
+        typedef scalar_matrix<T, ALLOC> self_type;
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using matrix_container<self_type>::operator ();
diff --git a/include/boost/numeric/ublas/matrix_expression.hpp b/include/boost/numeric/ublas/matrix_expression.hpp
index ac86f658..6da4e457 100644
--- a/include/boost/numeric/ublas/matrix_expression.hpp
+++ b/include/boost/numeric/ublas/matrix_expression.hpp
@@ -2482,7 +2482,7 @@ namespace boost { namespace numeric { namespace ublas {
     operator + (const matrix_expression<E1> &e1,
                 const matrix_expression<E2> &e2) {
         typedef typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
-                                                                              typename E2::value_type> >::expression_type expression_type;
+                                                                  typename E2::value_type> >::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -2494,7 +2494,7 @@ namespace boost { namespace numeric { namespace ublas {
     operator - (const matrix_expression<E1> &e1,
                 const matrix_expression<E2> &e2) {
         typedef typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
-                                                                               typename E2::value_type> >::expression_type expression_type;
+                                                                   typename E2::value_type> >::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -2506,7 +2506,7 @@ namespace boost { namespace numeric { namespace ublas {
     element_prod (const matrix_expression<E1> &e1,
                   const matrix_expression<E2> &e2) {
         typedef typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
-                                                                                    typename E2::value_type> >::expression_type expression_type;
+                                                                        typename E2::value_type> >::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -2514,11 +2514,11 @@ namespace boost { namespace numeric { namespace ublas {
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
     typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
-                                                       typename E2::value_type> >::result_type
+                                                         typename E2::value_type> >::result_type
     element_div (const matrix_expression<E1> &e1,
                  const matrix_expression<E2> &e2) {
         typedef typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
-                                                                                 typename E2::value_type> >::expression_type expression_type;
+                                                                     typename E2::value_type> >::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -3689,7 +3689,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef unknown_storage_tag storage_category;
         typedef row_major_tag orientation_category;
         typedef typename promote_traits<T1, T2>::promote_type promote_type;
-        typedef matrix_vector_binary1<E1, E2, matrix_vector_prod1<T1, T2, promote_type> > expression_type;
+        typedef matrix_vector_binary1<E1, E2, matrix_vector_prod1<E1, E2, promote_type> > expression_type;
 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
         typedef expression_type result_type;
 #else
@@ -3706,7 +3706,7 @@ namespace boost { namespace numeric { namespace ublas {
           unknown_storage_tag,
           row_major_tag) {
         typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
-                                                                  typename E2::value_type, E2>::expression_type expression_type;
+                                                      typename E2::value_type, E2>::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -3719,9 +3719,9 @@ namespace boost { namespace numeric { namespace ublas {
           const vector_expression<E2> &e2) {
         BOOST_STATIC_ASSERT (E2::complexity == 0);
         typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
-                                                                  typename E2::value_type, E2>::storage_category storage_category;
+                                                      typename E2::value_type, E2>::storage_category storage_category;
         typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
-                                                                  typename E2::value_type, E2>::orientation_category orientation_category;
+                                                      typename E2::value_type, E2>::orientation_category orientation_category;
         return prod (e1, e2, storage_category (), orientation_category ());
     }
 
@@ -3734,7 +3734,7 @@ namespace boost { namespace numeric { namespace ublas {
                unknown_storage_tag,
                row_major_tag) {
         typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                  typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
+                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -3747,9 +3747,9 @@ namespace boost { namespace numeric { namespace ublas {
                const vector_expression<E2> &e2) {
         BOOST_STATIC_ASSERT (E2::complexity == 0);
         typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                  typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
+                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
         typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                  typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
+                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
         return prec_prod (e1, e2, storage_category (), orientation_category ());
     }
 
@@ -4079,7 +4079,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef unknown_storage_tag storage_category;
         typedef column_major_tag orientation_category;
         typedef typename promote_traits<T1, T2>::promote_type promote_type;
-        typedef matrix_vector_binary2<E1, E2, matrix_vector_prod2<T1, T2, promote_type> > expression_type;
+        typedef matrix_vector_binary2<E1, E2, matrix_vector_prod2<E1, E2, promote_type> > expression_type;
 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
         typedef expression_type result_type;
 #else
@@ -4096,7 +4096,7 @@ namespace boost { namespace numeric { namespace ublas {
           unknown_storage_tag,
           column_major_tag) {
         typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
-                                                                  typename E2::value_type, E2>::expression_type expression_type;
+                                                      typename E2::value_type, E2>::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -4109,9 +4109,9 @@ namespace boost { namespace numeric { namespace ublas {
           const matrix_expression<E2> &e2) {
         BOOST_STATIC_ASSERT (E1::complexity == 0);
         typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
-                                                                  typename E2::value_type, E2>::storage_category storage_category;
+                                                      typename E2::value_type, E2>::storage_category storage_category;
         typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
-                                                                  typename E2::value_type, E2>::orientation_category orientation_category;
+                                                      typename E2::value_type, E2>::orientation_category orientation_category;
         return prod (e1, e2, storage_category (), orientation_category ());
     }
 
@@ -4124,7 +4124,7 @@ namespace boost { namespace numeric { namespace ublas {
                unknown_storage_tag,
                column_major_tag) {
         typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                  typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
+                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -4137,9 +4137,9 @@ namespace boost { namespace numeric { namespace ublas {
                const matrix_expression<E2> &e2) {
         BOOST_STATIC_ASSERT (E1::complexity == 0);
         typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                  typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
+                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
         typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                  typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
+                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
         return prec_prod (e1, e2, storage_category (), orientation_category ());
     }
 
@@ -4784,7 +4784,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef unknown_storage_tag storage_category;
         typedef unknown_orientation_tag orientation_category;
         typedef typename promote_traits<T1, T2>::promote_type promote_type;
-        typedef matrix_matrix_binary<E1, E2, matrix_matrix_prod<T1, T2, promote_type> > expression_type;
+        typedef matrix_matrix_binary<E1, E2, matrix_matrix_prod<E1, E2, promote_type> > expression_type;
 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
         typedef expression_type result_type;
 #else
@@ -4801,7 +4801,7 @@ namespace boost { namespace numeric { namespace ublas {
           unknown_storage_tag,
           unknown_orientation_tag) {
         typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
-                                                                 typename E2::value_type, E2>::expression_type expression_type;
+                                                     typename E2::value_type, E2>::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -4814,9 +4814,9 @@ namespace boost { namespace numeric { namespace ublas {
           const matrix_expression<E2> &e2) {
         BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
         typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
-                                                                 typename E2::value_type, E2>::storage_category storage_category;
+                                                     typename E2::value_type, E2>::storage_category storage_category;
         typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
-                                                                 typename E2::value_type, E2>::orientation_category orientation_category;
+                                                     typename E2::value_type, E2>::orientation_category orientation_category;
         return prod (e1, e2, storage_category (), orientation_category ());
     }
 
@@ -4829,7 +4829,7 @@ namespace boost { namespace numeric { namespace ublas {
                unknown_storage_tag,
                unknown_orientation_tag) {
         typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                 typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
+                                                     typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
@@ -4842,9 +4842,9 @@ namespace boost { namespace numeric { namespace ublas {
                const matrix_expression<E2> &e2) {
         BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
         typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                 typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
+                                                     typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
         typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
-                                                                 typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
+                                                     typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
         return prec_prod (e1, e2, storage_category (), orientation_category ());
     }
 
@@ -4925,25 +4925,25 @@ namespace boost { namespace numeric { namespace ublas {
 
     template<class E>
     BOOST_UBLAS_INLINE
-    typename matrix_scalar_unary_traits<E, matrix_norm_1<typename E::value_type> >::result_type
+    typename matrix_scalar_unary_traits<E, matrix_norm_1<E> >::result_type
     norm_1 (const matrix_expression<E> &e) {
-        typedef typename matrix_scalar_unary_traits<E, matrix_norm_1<typename E::value_type> >::expression_type expression_type;
+        typedef typename matrix_scalar_unary_traits<E, matrix_norm_1<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
     template<class E>
     BOOST_UBLAS_INLINE
-    typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<typename E::value_type> >::result_type
+    typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<E> >::result_type
     norm_frobenius (const matrix_expression<E> &e) {
-        typedef typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<typename E::value_type> >::expression_type expression_type;
+        typedef typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
     template<class E>
     BOOST_UBLAS_INLINE
-    typename matrix_scalar_unary_traits<E, matrix_norm_inf<typename E::value_type> >::result_type
+    typename matrix_scalar_unary_traits<E, matrix_norm_inf<E> >::result_type
     norm_inf (const matrix_expression<E> &e) {
-        typedef typename matrix_scalar_unary_traits<E, matrix_norm_inf<typename E::value_type> >::expression_type expression_type;
+        typedef typename matrix_scalar_unary_traits<E, matrix_norm_inf<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
diff --git a/include/boost/numeric/ublas/matrix_sparse.hpp b/include/boost/numeric/ublas/matrix_sparse.hpp
index 3b5fb893..70f6648f 100644
--- a/include/boost/numeric/ublas/matrix_sparse.hpp
+++ b/include/boost/numeric/ublas/matrix_sparse.hpp
@@ -2512,7 +2512,7 @@ namespace boost { namespace numeric { namespace ublas {
         // is_convertable (IA::size_type, TA::size_type)
         typedef typename IA::value_type size_type;
         // size_type for the data arrays.
-        typedef typename IA::size_type  array_size_type;
+        typedef typename IA::size_type array_size_type;
         // FIXME difference type for sparse storage iterators should it be in the container?
         typedef typename IA::difference_type difference_type;
         typedef T value_type;
@@ -3869,13 +3869,12 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using matrix_container<self_type>::operator ();
 #endif
-        // ISSUE require type consistency check
-        // is_convertable (IA::size_type, TA::size_type)
+        // ISSUE require type consistency check, is_convertable (IA::size_type, TA::size_type)
         typedef typename IA::value_type size_type;
+        // ISSUE difference_type cannot be deduced for sparse indices, we only know the value_type
+        typedef typename std::ptrdiff_t difference_type;
         // size_type for the data arrays.
-        typedef typename IA::size_type  array_size_type;
-        // FIXME difference type for sprase storage iterators should it be in the container?
-        typedef typename IA::difference_type difference_type;
+        typedef typename IA::size_type array_size_type;
         typedef T value_type;
         typedef const T &const_reference;
 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
@@ -3902,7 +3901,7 @@ namespace boost { namespace numeric { namespace ublas {
             storage_invariants ();
         }
         BOOST_UBLAS_INLINE
-        coordinate_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
+        coordinate_matrix (size_type size1, size_type size2, array_size_type non_zeros = 0):
             matrix_container<self_type> (),
             size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
@@ -3919,7 +3918,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
         template<class AE>
         BOOST_UBLAS_INLINE
-        coordinate_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
+        coordinate_matrix (const matrix_expression<AE> &ae, array_size_type non_zeros = 0):
             matrix_container<self_type> (),
             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
@@ -3993,9 +3992,9 @@ namespace boost { namespace numeric { namespace ublas {
         // Resizing
     private:
         BOOST_UBLAS_INLINE
-        size_type restrict_capacity (size_type non_zeros) const {
+        array_size_type restrict_capacity (array_size_type non_zeros) const {
             // minimum non_zeros
-            non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_));
+            non_zeros = (std::max) (non_zeros, array_size_type((std::min) (size1_, size2_)));
             // ISSUE no maximum as coordinate may contain inserted duplicates
             return non_zeros;
         }
@@ -4018,7 +4017,7 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Reserving
         BOOST_UBLAS_INLINE
-        void reserve (size_type non_zeros, bool preserve = true) {
+        void reserve (array_size_type non_zeros, bool preserve = true) {
             sort ();    // remove duplicate elements
             capacity_ = restrict_capacity (non_zeros);
             if (preserve) {
@@ -4368,7 +4367,7 @@ namespace boost { namespace numeric { namespace ublas {
                     }
                 } else /* if (direction < 0)  */ {
                     if (layout_type::fast_i ()) {
-                        if (it == index2_data_.begin () + zero_based (*itv))
+                        if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                             return const_iterator1 (*this, rank, i, j, itv, it);
                         i = zero_based (*(it - 1));
                     } else {
@@ -4409,7 +4408,7 @@ namespace boost { namespace numeric { namespace ublas {
                     }
                 } else /* if (direction < 0)  */ {
                     if (layout_type::fast_i ()) {
-                        if (it == index2_data_.begin () + zero_based (*itv))
+                        if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                             return iterator1 (*this, rank, i, j, itv, it);
                         i = zero_based (*(it - 1));
                     } else {
@@ -4450,7 +4449,7 @@ namespace boost { namespace numeric { namespace ublas {
                     }
                 } else /* if (direction < 0)  */ {
                     if (layout_type::fast_j ()) {
-                        if (it == index2_data_.begin () + zero_based (*itv))
+                        if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                             return const_iterator2 (*this, rank, i, j, itv, it);
                         j = zero_based (*(it - 1));
                     } else {
@@ -4491,7 +4490,7 @@ namespace boost { namespace numeric { namespace ublas {
                     }
                 } else /* if (direction < 0)  */ {
                     if (layout_type::fast_j ()) {
-                        if (it == index2_data_.begin () + zero_based (*itv))
+                        if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
                             return iterator2 (*this, rank, i, j, itv, it);
                         j = zero_based (*(it - 1));
                     } else {
diff --git a/include/boost/numeric/ublas/vector.hpp b/include/boost/numeric/ublas/vector.hpp
index 9e967d74..d77f97dc 100644
--- a/include/boost/numeric/ublas/vector.hpp
+++ b/include/boost/numeric/ublas/vector.hpp
@@ -569,18 +569,18 @@ namespace boost { namespace numeric { namespace ublas {
 
 
     // Zero vector class
-    template<class T>
+    template<class T, class ALLOC>
     class zero_vector:
-        public vector_container<zero_vector<T> > {
+        public vector_container<zero_vector<T, ALLOC> > {
 
         typedef const T *const_pointer;
-        typedef zero_vector<T> self_type;
+        typedef zero_vector<T, ALLOC> self_type;
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using vector_container<self_type>::operator ();
 #endif
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
+        typedef typename ALLOC::size_type size_type;
+        typedef typename ALLOC::difference_type difference_type;
         typedef T value_type;
         typedef const T &const_reference;
         typedef T &reference;
@@ -754,23 +754,23 @@ namespace boost { namespace numeric { namespace ublas {
         static const_value_type zero_;
     };
 
-    template<class T>
-    typename zero_vector<T>::const_value_type zero_vector<T>::zero_ (0);
+    template<class T, class ALLOC>
+    typename zero_vector<T, ALLOC>::const_value_type zero_vector<T, ALLOC>::zero_ (0);
 
 
     // Unit vector class
-    template<class T>
+    template<class T, class ALLOC>
     class unit_vector:
-        public vector_container<unit_vector<T> > {
+        public vector_container<unit_vector<T, ALLOC> > {
 
         typedef const T *const_pointer;
-        typedef unit_vector<T> self_type;
+        typedef unit_vector<T, ALLOC> self_type;
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using vector_container<self_type>::operator ();
 #endif
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
+        typedef typename ALLOC::size_type size_type;
+        typedef typename ALLOC::difference_type difference_type;
         typedef T value_type;
         typedef const T &const_reference;
         typedef T &reference;
@@ -966,25 +966,25 @@ namespace boost { namespace numeric { namespace ublas {
         static const_value_type one_;
     };
 
-    template<class T>
-    typename unit_vector<T>::const_value_type unit_vector<T>::zero_ (0);
-    template<class T>
-    typename unit_vector<T>::const_value_type unit_vector<T>::one_ (1);
+    template<class T, class ALLOC>
+    typename unit_vector<T, ALLOC>::const_value_type unit_vector<T, ALLOC>::zero_ (0);
+    template<class T, class ALLOC>
+    typename unit_vector<T, ALLOC>::const_value_type unit_vector<T, ALLOC>::one_ (1);
 
 
     // Scalar vector class
-    template<class T>
+    template<class T, class ALLOC>
     class scalar_vector:
-        public vector_container<scalar_vector<T> > {
+        public vector_container<scalar_vector<T, ALLOC> > {
 
         typedef const T *const_pointer;
-        typedef scalar_vector<T> self_type;
+        typedef scalar_vector<T, ALLOC> self_type;
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
         using vector_container<self_type>::operator ();
 #endif
-        typedef std::size_t size_type;
-        typedef std::ptrdiff_t difference_type;
+        typedef typename ALLOC::size_type size_type;
+        typedef typename ALLOC::difference_type difference_type;
         typedef T value_type;
         typedef const T &const_reference;
         typedef T &reference;
diff --git a/include/boost/numeric/ublas/vector_expression.hpp b/include/boost/numeric/ublas/vector_expression.hpp
index f5588de4..a466ae68 100644
--- a/include/boost/numeric/ublas/vector_expression.hpp
+++ b/include/boost/numeric/ublas/vector_expression.hpp
@@ -1427,8 +1427,6 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename E::const_iterator::iterator_category iterator_category;
         typedef vector_scalar_unary<E, F> self_type;
     public:
-        typedef typename F::size_type size_type;
-        typedef typename F::difference_type difference_type;
         typedef typename F::result_type value_type;
         typedef const self_type const_closure_type;
         typedef const_closure_type closure_type;
@@ -1501,9 +1499,9 @@ namespace boost { namespace numeric { namespace ublas {
     // sum v = sum (v [i])
     template<class E>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_unary_traits<E, vector_sum<typename E::value_type> >::result_type
+    typename vector_scalar_unary_traits<E, vector_sum<E> >::result_type
     sum (const vector_expression<E> &e) {
-        typedef typename vector_scalar_unary_traits<E, vector_sum<typename E::value_type> >::expression_type expression_type;
+        typedef typename vector_scalar_unary_traits<E, vector_sum<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
@@ -1511,9 +1509,9 @@ namespace boost { namespace numeric { namespace ublas {
     // complex: norm_1 v = sum (abs (real (v [i])) + abs (imag (v [i])))
     template<class E>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_unary_traits<E, vector_norm_1<typename E::value_type> >::result_type
+    typename vector_scalar_unary_traits<E, vector_norm_1<E> >::result_type
     norm_1 (const vector_expression<E> &e) {
-        typedef typename vector_scalar_unary_traits<E, vector_norm_1<typename E::value_type> >::expression_type expression_type;
+        typedef typename vector_scalar_unary_traits<E, vector_norm_1<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
@@ -1521,9 +1519,9 @@ namespace boost { namespace numeric { namespace ublas {
     // complex: norm_2 v = sqrt (sum (v [i] * conj (v [i])))
     template<class E>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_unary_traits<E, vector_norm_2<typename E::value_type> >::result_type
+    typename vector_scalar_unary_traits<E, vector_norm_2<E> >::result_type
     norm_2 (const vector_expression<E> &e) {
-        typedef typename vector_scalar_unary_traits<E, vector_norm_2<typename E::value_type> >::expression_type expression_type;
+        typedef typename vector_scalar_unary_traits<E, vector_norm_2<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
@@ -1531,18 +1529,18 @@ namespace boost { namespace numeric { namespace ublas {
     // complex: norm_inf v = maximum (maximum (abs (real (v [i])), abs (imag (v [i]))))
     template<class E>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_unary_traits<E, vector_norm_inf<typename E::value_type> >::result_type
+    typename vector_scalar_unary_traits<E, vector_norm_inf<E> >::result_type
     norm_inf (const vector_expression<E> &e) {
-        typedef typename vector_scalar_unary_traits<E, vector_norm_inf<typename E::value_type> >::expression_type expression_type;
+        typedef typename vector_scalar_unary_traits<E, vector_norm_inf<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
     // real: index_norm_inf v = minimum (i: abs (v [i]) == maximum (abs (v [i])))
     template<class E>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_unary_traits<E, vector_index_norm_inf<typename E::value_type> >::result_type
+    typename vector_scalar_unary_traits<E, vector_index_norm_inf<E> >::result_type
     index_norm_inf (const vector_expression<E> &e) {
-        typedef typename vector_scalar_unary_traits<E, vector_index_norm_inf<typename E::value_type> >::expression_type expression_type;
+        typedef typename vector_scalar_unary_traits<E, vector_index_norm_inf<E> >::expression_type expression_type;
         return expression_type (e ());
     }
 
@@ -1560,8 +1558,6 @@ namespace boost { namespace numeric { namespace ublas {
         typedef vector_scalar_binary<E1, E2, F> self_type;
     public:
         static const unsigned complexity = 1;
-        typedef typename F::size_type size_type;
-        typedef typename F::difference_type difference_type;
         typedef typename F::result_type value_type;
         typedef const self_type const_closure_type;
         typedef const_closure_type closure_type;
@@ -1639,31 +1635,27 @@ namespace boost { namespace numeric { namespace ublas {
     // inner_prod (v1, v2) = sum (v1 [i] * v2 [i])
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
-                                                                   typename E2::value_type,
+    typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
                                                                    typename promote_traits<typename E1::value_type,
                                                                                            typename E2::value_type>::promote_type> >::result_type
     inner_prod (const vector_expression<E1> &e1,
                 const vector_expression<E2> &e2) {
-        typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
-                                                                                           typename E2::value_type,
-                                                                                           typename promote_traits<typename E1::value_type,
-                                                                                                                               typename E2::value_type>::promote_type> >::expression_type expression_type;
+        typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
+                                                                   typename promote_traits<typename E1::value_type,
+                                                                                           typename E2::value_type>::promote_type> >::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }
 
     template<class E1, class E2>
     BOOST_UBLAS_INLINE
-    typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
-                                                                   typename E2::value_type,
+    typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
                                                                    typename type_traits<typename promote_traits<typename E1::value_type,
                                                                                                                 typename E2::value_type>::promote_type>::precision_type> >::result_type
     prec_inner_prod (const vector_expression<E1> &e1,
                      const vector_expression<E2> &e2) {
-        typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<typename E1::value_type,
-                                                                                           typename E2::value_type,
-                                                                                           typename type_traits<typename promote_traits<typename E1::value_type,
-                                                                                                                                                                typename E2::value_type>::promote_type>::precision_type> >::expression_type expression_type;
+        typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
+                                                                   typename type_traits<typename promote_traits<typename E1::value_type,
+                                                                                                                typename E2::value_type>::promote_type>::precision_type> >::expression_type expression_type;
         return expression_type (e1 (), e2 ());
     }