diff --git a/bench5/assignment_bench.cpp b/bench5/assignment_bench.cpp
new file mode 100644
index 00000000..532f3791
--- /dev/null
+++ b/bench5/assignment_bench.cpp
@@ -0,0 +1,141 @@
+//
+//  Copyright (c) 2010 Athanasios Iliopoulos
+//
+//  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)
+//
+
+#include <boost/numeric/ublas/assignment.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/timer.hpp>
+
+using namespace boost::numeric::ublas;
+
+int main() {
+
+    boost::timer timer;
+
+    unsigned int iterations = 1000000000;
+    double elapsed_exp, elapsed_assigner;
+
+    std::cout << "Ublas vector<double> Benchmarks------------------------ " << "\n";
+
+    {
+    std::cout << "Size 2 vector: " << "\n";
+    vector<double> a(2);
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++) {
+        a(0)=0; a(1)=1;
+    }
+    elapsed_exp = timer.elapsed();
+    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++)
+        a <<= 0, 1;
+    elapsed_assigner = timer.elapsed();
+    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
+    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
+    }
+
+    {
+    std::cout << "Size 3 vector: " << "\n";
+    vector<double> a(3);
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++) {
+        a(0)=0; a(1)=1; a(2)=2;
+    }
+    elapsed_exp = timer.elapsed();
+    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++)
+        a <<= 0, 1, 2;
+    elapsed_assigner = timer.elapsed();
+    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
+    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
+    }
+
+    iterations = 100000000;
+
+    {
+    std::cout << "Size 8 vector: " << "\n";
+    vector<double> a(8);
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++) {
+        a(0)=0; a(1)=1; a(2)=2; a(3)=3; a(4)=4; a(5)=5; a(6)=6; a(7)=7;
+    }
+    elapsed_exp = timer.elapsed();
+    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++)
+        a <<= 0, 1, 2, 3, 4, 5, 6, 7;
+    elapsed_assigner = timer.elapsed();
+    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
+    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
+    }
+
+
+    std::cout << "Ublas matrix<double> Benchmarks------------------------ " << "\n";
+
+    iterations = 200000000;
+    {
+    std::cout << "Size 3x3 matrix: " << "\n";
+    matrix<double> a(3,3);
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++) {
+        a(0,0)=0; a(0,1)=1; a(0,2)=2;
+        a(1,0)=3; a(1,1)=4; a(1,2)=5;
+        a(2,0)=6; a(2,1)=7; a(2,2)=8;
+    }
+    elapsed_exp = timer.elapsed();
+    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++)
+        a <<= 0, 1, 2, 3, 4, 5, 6, 7, 8;
+    elapsed_assigner = timer.elapsed();
+    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
+    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
+    }
+
+    std::cout << "Size 2x2 matrix: " << "\n";
+    iterations = 500000000;
+    {
+    matrix<double> a(2,2);
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++) {
+        a(0,0)=0; a(0,1)=1;
+        a(1,0)=3; a(1,1)=4;
+    }
+    elapsed_exp = timer.elapsed();
+    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++)
+        a <<= 0, 1, 3, 4;
+    elapsed_assigner = timer.elapsed();
+    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
+
+    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
+
+    timer.restart();
+    for(unsigned int i=0; i!=iterations; i++)
+        a <<= traverse_policy::by_row_no_wrap(), 0, 1, next_row(), 3, 4;
+    elapsed_assigner = timer.elapsed();
+    std::cout << "Assigner time no_wrap: " << elapsed_assigner << " secs" << "\n";
+    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
+    }
+
+    return 0;
+}
+
diff --git a/doc/samples/assignment_examples.cpp b/doc/samples/assignment_examples.cpp
new file mode 100644
index 00000000..bfad1f54
--- /dev/null
+++ b/doc/samples/assignment_examples.cpp
@@ -0,0 +1,319 @@
+//
+//  Copyright (c) 2010 Athanasios Iliopoulos
+//
+//  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)
+//
+
+#include <boost/numeric/ublas/assignment.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector_sparse.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+
+using namespace boost::numeric::ublas;
+
+int main() {
+        // Simple vector fill
+    vector<double> a(3);
+    a <<= 0, 1, 2;
+    std::cout << a << std::endl;
+    // [ 0 1 2]
+
+    // Vector from vector
+    vector<double> b(7);
+    b <<= a, 10, a;
+    std::cout << b << std::endl;
+    // [ 0 1 2 10 0 1 2]
+
+    // Simple matrix fill
+    matrix<double> A(3,3);
+    A <<= 0, 1, 2,
+         3, 4, 5,
+         6, 7, 8;
+    std::cout << A << std::endl;
+    // [ 0 1 2 ]
+    // [ 3 4 5 ]
+    // [ 6 7 8 ]
+
+    // Matrix from vector
+    A <<= 0, 1, 2,
+         3, 4, 5,
+         a;
+    std::cout << A << std::endl;
+    // [ 0 1 2 ]
+    // [ 3 4 5 ]
+    // [ 0 1 2 ]
+
+    // Matrix from vector - column assignment
+    A <<= move(0,2), traverse_policy::by_column(),
+         a;
+    std::cout << A << std::endl;
+    // [ 0 1 0 ]
+    // [ 3 4 1 ]
+    // [ 0 1 2 ]
+
+    // Another matrix from vector example (watch the wraping);
+    vector<double> c(9); c <<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    A <<= c;
+    std::cout << A << std::endl;
+    // [ 1 2 3 ]
+    // [ 4 5 6 ]
+    // [ 7 8 9 ]
+
+    // If for performance(Benchmarks are not definite about that) or consistency reasons you need to disable wraping:
+    static next_row_manip endr; //This can be defined globally
+    A <<= traverse_policy::by_row_no_wrap(),
+            1, 2, 3, endr,
+            4, 5, 6, endr,
+            7, 8, 9, endr;
+    // [ 1 2 3 ]
+    // [ 4 5 6 ]
+    // [ 7 8 9 ]
+    // If by default you need to disable wraping define
+    // BOOST_UBLAS_DEFAULT_NO_WRAP_POLICY, in the compilation options,
+    // so that you avoid typing the "traverse_policy::by_row_no_wrap()".
+
+    //  Plus and minus assign:
+    A <<= fill_policy::index_plus_assign(),
+         3,2,1;
+    std::cout << A << std::endl;
+    // [ 4 4 4 ]
+    // [ 4 5 6 ]
+    // [ 7 8 9 ]
+
+    // Matrix from proxy
+    A <<= 0, 1, 2,
+         project(b, range(3,6)),
+         a;
+    std::cout << A << std::endl;
+    // [ 0 1 2 ]
+    // [10 0 1 ]
+    // [ 6 7 8 ]
+
+    // Matrix from matrix
+    matrix<double> B(6,6);
+    B <<= A, A,
+         A, A;
+    std::cout << B << std::endl;
+    // [ A A ]
+    // [ A A ]
+
+    // Matrix range (vector is similar)
+    B = zero_matrix<double>(6,6);
+    matrix_range<matrix<double> > mrB (B, range (1, 4), range (1, 4));
+    mrB <<= 1,2,3,4,5,6,7,8,9;
+    std::cout << B << std::endl;
+    // [ 0 0 0 0 0 0]
+    // [ 0 1 2 3 0 0]
+    // [ 0 4 5 6 0 0]
+    // [ 0 0 0 0 0 0]
+    // [ 0 0 0 0 0 0]
+    // [ 0 0 0 0 0 0]
+
+    // Horizontal concatenation can be achieved using this trick:
+    matrix<double> BH(3,9);
+    BH <<= A, A, A;
+    std::cout << BH << std::endl;
+    // [ A A A]
+
+    // Vertical concatenation can be achieved using this trick:
+    matrix<double> BV(9,3);
+    BV <<= A,
+          A,
+          A;
+    std::cout << BV << std::endl;
+    // [ A ]
+    // [ A ]
+    // [ A ]
+
+    // Watch the difference when assigning matrices for different traverse policies:
+    matrix<double> BR(9,9, 0);
+    BR <<= traverse_policy::by_row(), // This is the default, so this might as well be omitted.
+          A, A, A;
+    std::cout << BR << std::endl;
+    // [ A A A]
+    // [ 0 0 0]
+    // [ 0 0 0]
+
+    matrix<double> BC(9,9, 0);
+    BC <<= traverse_policy::by_column(),
+          A, A, A;
+    std::cout << BC << std::endl;
+    // [ A 0 0]
+    // [ A 0 0]
+    // [ A 0 0]
+
+    // The following will throw a run-time exception in debug mode (matrix mid-assignment wrap is not allowed) :
+    // matrix<double> C(7,7);
+    // C <<= A, A, A;
+
+    // Matrix from matrix with index manipulators
+    matrix<double> C(6,6,0);
+    C <<= A, move(3,0), A;
+    // [ A 0 ]
+    // [ 0 A ]
+
+    // A faster way for to construct this dense matrix.
+    matrix<double> D(6,6);
+    D <<= A, zero_matrix<double>(3,3),
+         zero_matrix<double>(3,3), A;
+    // [ A 0 ]
+    // [ 0 A ]
+
+    // The next_row and next_column index manipulators:
+    // note: next_row and next_column functions return
+    // a next_row_manip and and next_column_manip object.
+    // This is the manipulator we used earlier when we disabled
+    // wrapping.
+    matrix<double> E(2,4,0);
+    E <<= 1, 2, next_row(),
+         3, 4, next_column(),5;
+    std::cout << E << std::endl;
+    // [ 1 2 0 5 ]
+    // [ 3 4 0 0 ]
+
+    // The begin1 (moves to the begining of the column) index manipulator, begin2 does the same for the row:
+    matrix<double> F(2,4,0);
+    F <<= 1, 2, next_row(),
+         3, 4, begin1(),5;
+    std::cout << F << std::endl;
+    // [ 1 2 5 0 ]
+    // [ 3 4 0 0 ]
+
+    // The move (relative) and move_to(absolute) index manipulators (probably the most useful manipulators):
+    matrix<double> G(2,4,0);
+    G <<= 1, 2, move(0,1), 3,
+         move_to(1,3), 4;
+    std::cout << G << std::endl;
+    // [ 1 2 0 3 ]
+    // [ 0 0 0 4 ]
+
+    // Static equivallents (faster) when sizes are known at compile time:
+    matrix<double> Gs(2,4,0);
+    Gs <<= 1, 2, move<0,1>(), 3,
+         move_to<1,3>(), 4;
+    std::cout << Gs << std::endl;
+    // [ 1 2 0 3 ]
+    // [ 0 0 0 4 ]
+
+    // Choice of traverse policy (default is "row by row" traverse):
+
+    matrix<double> H(2,4,0);
+    H <<= 1, 2, 3, 4,
+         5, 6, 7, 8;
+    std::cout << H << std::endl;
+    // [ 1 2 3 4 ]
+    // [ 5 6 7 8 ]
+
+    H <<= traverse_policy::by_column(),
+        1, 2, 3, 4,
+        5, 6, 7, 8;
+    std::cout << H << std::endl;
+    // [ 1 3 5 7 ]
+    // [ 2 4 6 8 ]
+
+    // traverse policy can be changed mid assignment if desired.
+     matrix<double> H1(4,4,0);
+     H1 <<= 1, 2, 3, traverse_policy::by_column(), 1, 2, 3;
+
+    std::cout << H << std::endl;
+    // [1 2 3 1]
+    // [0 0 0 2]
+    // [0 0 0 3]
+    // [0 0 0 0]
+
+    // note: fill_policy and traverse_policy are namespaces, so you can use them
+    // by a using statement.
+
+    // For compressed and coordinate matrix types a push_back or insert fill policy can be chosen for faster assginment:
+    compressed_matrix<double> I(2, 2);
+    I <<=    fill_policy::sparse_push_back(),
+            0, 1, 2, 3;
+    std::cout << I << std::endl;
+    // [ 0 1 ]
+    // [ 2 3 ]
+
+    coordinate_matrix<double> J(2,2);
+    J<<=fill_policy::sparse_insert(),
+        1, 2, 3, 4;
+    std::cout << J << std::endl;
+    // [ 1 2 ]
+    // [ 3 4 ]
+
+    // A sparse matrix from another matrix works as with other types.
+    coordinate_matrix<double> K(3,3);
+    K<<=fill_policy::sparse_insert(),
+        J;
+    std::cout << K << std::endl;
+    // [ 1 2 0 ]
+    // [ 3 4 0 ]
+    // [ 0 0 0 ]
+
+    // Be careful this will not work:
+    //compressed_matrix<double> J2(4,4);
+    //J2<<=fill_policy::sparse_push_back(),
+     //   J,J;
+    // That's because the second J2's elements
+    // are attempted to be assigned at positions
+    // that come before the elements already pushed.
+    // Unfortunatelly that's the only thing you can do in this case
+    // (or of course make a custom agorithm):
+    compressed_matrix<double> J2(4,4);
+    J2<<=fill_policy::sparse_push_back(),
+        J, fill_policy::sparse_insert(),
+        J;
+
+    std::cout << J2 << std::endl;
+    // [  J   J  ]
+    // [ 0 0 0 0 ]
+    // [ 0 0 0 0 ]
+
+    // A different traverse policy doesn't change the result, only they order it is been assigned.
+    coordinate_matrix<double> L(3,3);
+    L<<=fill_policy::sparse_insert(), traverse_policy::by_column(),
+        J;
+    std::cout << L << std::endl;
+    // (same as previous)
+    // [ 1 2 0 ]
+    // [ 3 4 0 ]
+    // [ 0 0 0 ]
+
+    typedef coordinate_matrix<double>::size_type cmst;
+    const cmst size = 30;
+    //typedef fill_policy::sparse_push_back spb;
+    // Although the above could have been used the following is may be faster if
+    //  you use the policy often and for relatively small containers.
+    static fill_policy::sparse_push_back spb;
+
+    // A block diagonal sparse using a loop:
+    compressed_matrix<double> M(size, size, 4*15);
+    for (cmst i=0; i!=size; i+=J.size1())
+        M <<= spb, move_to(i,i), J;
+
+
+    // If typedef was used above the last expression should start
+    // with M <<= spb()...
+
+    // Displaying so that blocks can be easily seen:
+    for (unsigned int i=0; i!=M.size1(); i++) {
+        std::cout << M(i,0);
+        for (unsigned int j=1; j!=M.size2(); j++) std::cout << ", " << M(i,j);
+        std::cout << "\n";
+    }
+    // [ J 0 0 0 ... 0]
+    // [ 0 J 0 0 ... 0]
+    // [ 0 . . . ... 0]
+    // [ 0 0 ... 0 0 J]
+
+
+    // A "repeat" trasverser may by provided so that this becomes faster and an on-liner like:
+    // M <<= spb, repeat(0, size, J.size1(), 0, size, J.size1()), J;
+    // An alternate would be to create a :repeater" matrix and vector expression that can be used in other places as well. The latter is probably better,
+    return 0;
+}
+
diff --git a/test/test_assignment.cpp b/test/test_assignment.cpp
new file mode 100644
index 00000000..cbfe0e56
--- /dev/null
+++ b/test/test_assignment.cpp
@@ -0,0 +1,811 @@
+//
+//  Copyright (c) 2010 Athanasios Iliopoulos
+//
+//  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)
+//
+
+#include <boost/numeric/ublas/assignment.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <boost/numeric/ublas/vector_sparse.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include "libs/numeric/ublas/test/utils.hpp"
+#include <boost/timer.hpp>
+#include <ctime>
+
+using namespace boost::numeric::ublas;
+
+namespace tans {
+template <class AE>
+typename AE::value_type mean_square(const matrix_expression<AE> &me) {
+    typename AE::value_type s(0);
+    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));
+        }
+    }
+    return s/me().size1()*me().size2();
+}
+
+
+template <class AE>
+typename AE::value_type mean_square(const vector_expression<AE> &ve) {
+    // We could have use norm2 here, but ublas' ABS does not support unsigned types.
+    typename AE::value_type s(0);
+    typename AE::size_type i;
+    for (i=0; i!= ve().size(); i++) {
+            s+=std::fabs(ve()(i));
+        }
+    return s/ve().size();
+}
+const double TOL=0.0;
+
+}
+
+template <class V>
+bool test_vector() {
+    bool pass = true;
+    using namespace tans;
+
+    V a(3), ra(3);
+    a <<=  1, 2, 3;
+    ra(0) = 1; ra(1) = 2; ra(2) = 3;
+    pass &= (mean_square(a-ra)<=TOL);
+
+    V b(7), rb(7);
+    b<<= a, 10, a;
+    rb(0) = 1; rb(1) = 2; rb(2) = 3; rb(3)=10, rb(4)= 1; rb(5)=2; rb(6)=3;
+    pass &= (mean_square(b-rb)<=TOL);
+
+    V c(6), rc(6);
+    c <<= 1, move(2), 3 ,4, 5, move(-5), 10, 10;
+    rc(0) = 1; rc(1) = 10; rc(2) = 10; rc(3) = 3; rc(4) = 4; rc(5) = 5;
+    pass &= (mean_square(c-rc)<=TOL);
+
+    V d(6), rd(6);
+    d <<= 1, move_to(3), 3 ,4, 5, move_to(1), 10, 10;
+    rd(0) = 1; rd(1) = 10; rd(2) = 10; rd(3) = 3; rd(4) = 4; rd(5) = 5;
+    pass &= (mean_square(d-rd)<=TOL);
+
+    {
+    V c(6), rc(6);
+    c <<= 1, move<2>(), 3 ,4, 5, move<-5>(), 10, 10;
+    rc(0) = 1; rc(1) = 10; rc(2) = 10; rc(3) = 3; rc(4) = 4; rc(5) = 5;
+    pass &= (mean_square(c-rc)<=TOL);
+
+    V d(6), rd(6);
+    d <<= 1, move_to<3>(), 3 ,4, 5, move_to<1>(), 10, 10;
+    rd(0) = 1; rd(1) = 10; rd(2) = 10; rd(3) = 3; rd(4) = 4; rd(5) = 5;
+    pass &= (mean_square(d-rd)<=TOL);
+    }
+
+
+    {
+    V f(6), rf(6);
+    f <<= 5, 5, 5, 5, 5, 5;
+    V fa(3); fa<<= 1, 2, 3;
+    f <<= fill_policy::index_plus_assign(), fa;
+    rf <<= 6,7,8, 5, 5, 5;
+    pass &= (mean_square(f-rf)<=TOL);
+    }
+
+    {
+    V f(6), rf(6);
+    f <<= 5, 5, 5, 5, 5, 5;
+    V fa(3); fa<<= 1, 2, 3;
+    f <<= fill_policy::index_minus_assign(), fa;
+    rf <<= 4,3,2, 5, 5, 5;
+    pass &= (mean_square(f-rf)<=TOL);
+    }
+
+    return pass;
+}
+
+template <class V>
+bool test_vector_sparse_push_back() {
+    bool pass = true;
+    using namespace tans;
+
+    V a(3), ra(3);
+    a <<= fill_policy::sparse_push_back(), 1, 2, 3;
+    ra(0) = 1; ra(1) = 2; ra(2) = 3;
+    pass &= (mean_square(a-ra)<=TOL);
+
+    V b(7), rb(7);
+    b<<= fill_policy::sparse_push_back(), a, 10, a;
+    rb(0) = 1; rb(1) = 2; rb(2) = 3; rb(3)=10, rb(4)= 1; rb(5)=2; rb(6)=3;
+    pass &= (mean_square(b-rb)<=TOL);
+
+    V c(6), rc(6);
+    c <<= fill_policy::sparse_push_back(), 1, move(2), 3 ,4, 5; // Move back (i.e. negative is dangerous for push_back)
+    rc(0) = 1; rc(1) = 0; rc(2) = 0; rc(3) = 3; rc(4) = 4; rc(5) = 5;
+    pass &= (mean_square(c-rc)<=TOL);
+
+    V d(6), rd(6);
+    d <<= fill_policy::sparse_push_back(), 1, move_to(3), 3 ,4, 5; // Move back (i.e. before current index is dangerous for push_back)
+    rd(0) = 1; rd(1) = 0; rd(2) = 0; rd(3) = 3; rd(4) = 4; rd(5) = 5;
+    pass &= (mean_square(d-rd)<=TOL);
+
+    V e(6), re(6);
+    e <<= fill_policy::sparse_push_back(), 1, move_to(3), 3 ,4, 5, fill_policy::sparse_insert(), move_to(1), 10, 10; // If you want to move back, use this
+    re(0) = 1; re(1) = 10; re(2) = 10; re(3) = 3; re(4) = 4; re(5) = 5;
+    pass &= (mean_square(e-re)<=TOL);
+
+    return pass;
+}
+
+
+template <class V>
+bool test_vector_sparse_insert() {
+    bool pass = true;
+    using namespace tans;
+
+    V a(3), ra(3);
+    a <<= fill_policy::sparse_insert(), 1, 2, 3;
+    ra(0) = 1; ra(1) = 2; ra(2) = 3;
+    pass &= (mean_square(a-ra)<=TOL);
+
+    V b(7), rb(7);
+    b<<= fill_policy::sparse_insert(), a, 10, a;
+    rb(0) = 1; rb(1) = 2; rb(2) = 3; rb(3)=10, rb(4)= 1; rb(5)=2; rb(6)=3;
+    pass &= (mean_square(b-rb)<=TOL);
+
+    V c(6), rc(6);
+    c <<= fill_policy::sparse_insert(), 1, move(2), 3 ,4, 5, move(-5), 10, 10; // Move back (i.e. negative is dangerous for sparse)
+    rc(0) = 1; rc(1) = 10; rc(2) = 10; rc(3) = 3; rc(4) = 4; rc(5) = 5;
+    pass &= (mean_square(c-rc)<=TOL);
+
+
+    V d(6), rd(6);
+    d <<= fill_policy::sparse_insert(), 1, move_to(3), 3 ,4, 5, move_to(1), 10, 10; // Move back (i.e.before is dangerous for sparse)
+    rd(0) = 1; rd(1) = 10; rd(2) = 10; rd(3) = 3; rd(4) = 4; rd(5) = 5;
+    pass &= (mean_square(d-rd)<=TOL);
+
+
+    return pass;
+}
+
+
+template <class V>
+bool test_matrix() {
+    bool pass = true;
+    using namespace tans;
+
+    V A(3,3), RA(3,3);
+    A <<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    RA(0,0)= 1; RA(0,1)=2; RA(0,2)=3;
+    RA(1,0)= 4; RA(1,1)=5; RA(1,2)=6;
+    RA(2,0)= 7; RA(2,1)=8; RA(2,2)=9;
+    pass &= (mean_square(A-RA)<=TOL);
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(3);
+    b<<= 4,5,6;
+    B<<= 1, 2, 3, b, 7, project(b, range(1,3));
+    RB<<=1, 2, 3, 4, 5, 6, 7, 5, 6; // If the first worked we can now probably use it.
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(3);
+    b<<= 4,5,6;
+    B<<= move(1,0), b, move_to(0,0), 1, 2, 3, move(1,0), 7, project(b, range(1,3));
+    RB<<=1, 2, 3, 4, 5, 6, 7, 5, 6;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(9);
+    b<<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    B<<=b;
+    RB<<=1, 2, 3, 4, 5, 6, 7, 8, 9;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<=    2, 3,
+            4, 5;
+    B<<= C,C,
+        C,C;
+    RB <<=   2,3,2,3,
+            4,5,4,5,
+            2,3,2,3,
+            4,5,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<= C, zero_matrix<typename V::value_type>(2,2),
+        zero_matrix<typename V::value_type>(2,2), C;
+    RB<<=    2,3,0,0,
+            4,5,0,0,
+            0,0,2,3,
+            0,0,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<= C, zero_matrix<typename V::value_type>(2,2),
+        zero_matrix<typename V::value_type>(2,2), C;
+    RB<<=    2,3,0,0,
+            4,5,0,0,
+            0,0,2,3,
+            0,0,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4); // We need that because of the non-zero instatiation of dense types.
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<= move(1,1), C;
+    RB<<=    0,0,0,0,
+            0,2,3,0,
+            0,4,5,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<= move_to(0,1), 2, 3, next_row(), 1, 2, next_row(), 4, 5;
+    RB<<=    0,2,3,0,
+            1,2,0,0,
+            4,5,0,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=traverse_policy::by_column(), move_to(0,1), 2, 3, 6, next_column(), 4, 5;
+    RB<<=    0,2,4,0,
+            0,3,5,0,
+            0,6,0,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=traverse_policy::by_column(), move_to(0,1), 2, 3, next_row(), traverse_policy::by_row(), 4, 5;
+    RB<<=    0,2,0,0,
+            0,3,0,0,
+            0,0,0,0,
+            4,5,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=traverse_policy::by_column(), move_to(0,1), 2, 3, begin2(), traverse_policy::by_row(), 4, 5, 6, 7, 8;
+    RB<<=    0,2,0,0,
+            0,3,0,0,
+            4,5,6,7,
+            8,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=traverse_policy::by_column(), move_to(0,1), 2, 3, begin2(), traverse_policy::by_row(), 4, 5, 6, 7, 8,9, begin1(), 1, 2;
+    RB<<=    0,2,1,2,
+            0,3,0,0,
+            4,5,6,7,
+            8,9,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = scalar_matrix<typename V::value_type>(4,4,1);
+    V C(2,2);
+    C <<= 1, 2, 3, 4;
+    B<<= fill_policy::index_plus_assign(), move(1,1), C;
+    RB<<=    1,1,1,1,
+            1,2,3,1,
+            1,4,5,1,
+            1,1,1,1;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = scalar_matrix<typename V::value_type>(4,4,5);
+    V C(2,2);
+    C <<= 1, 2, 3, 4;
+    B<<= fill_policy::index_minus_assign(), move(1,1), C;
+    RB<<=    5,5,5,5,
+            5,4,3,5,
+            5,2,1,5,
+            5,5,5,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+
+    return pass;
+}
+
+template <class V>
+bool test_matrix_sparse_push_back() {
+    bool pass = true;
+    using namespace tans;
+
+    V A(3,3), RA(3,3);
+    A <<= fill_policy::sparse_push_back(), 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    RA(0,0)= 1; RA(0,1)=2; RA(0,2)=3;
+    RA(1,0)= 4; RA(1,1)=5; RA(1,2)=6;
+    RA(2,0)= 7; RA(2,1)=8; RA(2,2)=9;
+    pass &= (mean_square(A-RA)<=TOL);
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(3);
+    b<<= 4,5,6;
+    B<<=fill_policy::sparse_push_back(), 1, 2, 3, b, 7, project(b, range(1,3));
+    RB<<= 1, 2, 3, 4, 5, 6, 7, 5, 6; // If the first worked we can now probably use it.
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(3);
+    b<<= 4,5,6;
+    B<<=fill_policy::sparse_push_back(), move(1,0), b, fill_policy::sparse_insert(), move_to(0,0), 1, 2, 3, move(1,0), 7, project(b, range(1,3));
+    RB<<=1, 2, 3, 4, 5, 6, 7, 5, 6;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(9);
+    b<<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    B<<=b;
+    RB<<=1, 2, 3, 4, 5, 6, 7, 8, 9;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<=    2, 3,
+            4, 5;
+    // It might get complicated for sparse push_back, this must go into the tutorial. (This way is not convient nor fast)
+    B<<=fill_policy::sparse_push_back(), C, move_to(2,2), C, fill_policy::sparse_insert(), move_to(0,2), C, C;
+    RB <<=   2,3,2,3,
+            4,5,4,5,
+            2,3,2,3,
+            4,5,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<=fill_policy::sparse_push_back(), C, move_to(2,2), C;
+    RB<<=    2,3,0,0,
+            4,5,0,0,
+            0,0,2,3,
+            0,0,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<=fill_policy::sparse_push_back(), move(1,1), C;
+    RB<<=    0,0,0,0,
+            0,2,3,0,
+            0,4,5,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_push_back(), move_to(0,1), 2, 3, next_row(), 1, 2, next_row(), 4, 5;
+    RB<<=    0,2,3,0,
+            1,2,0,0,
+            4,5,0,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+    // The next will not work with sparse push_back because elements that are prior to the ones already in are attempted to be added
+/*
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_push_back(),traverse_policy::by_column(), move_to(0,1), 2, 3, 6, next_column(), 4, 5;
+    RB<<=    0,2,4,0,
+            0,3,5,0,
+            0,6,0,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+*/
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_push_back(),traverse_policy::by_column(), move_to(0,1), 2, 3, next_row(), traverse_policy::by_row(), 4, 5;
+    RB<<=    0,2,0,0,
+            0,3,0,0,
+            0,0,0,0,
+            4,5,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_push_back(),traverse_policy::by_column(), move_to(0,1), 2, 3, begin2(), traverse_policy::by_row(), 4, 5, 6, 7, 8;
+    RB<<=    0,2,0,0,
+            0,3,0,0,
+            4,5,6,7,
+            8,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    // The next will not work with sparse push_back because elements that are prior to the ones already in are attempted to be added
+/*
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_push_back(),traverse_policy::by_column(), move_to(0,1), 2, 3, begin2(), traverse_policy::by_row(), 4, 5, 6, 7, 8,9, begin1(), 1, 2;
+    RB<<=    0,2,1,2,
+            0,3,0,0,
+            4,5,6,7,
+            8,9,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+*/
+    return pass;
+}
+
+template <class V>
+bool test_matrix_sparse_insert() {
+    bool pass = true;
+    using namespace tans;
+
+    V A(3,3), RA(3,3);
+    A <<= fill_policy::sparse_insert(), 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    RA(0,0)= 1; RA(0,1)=2; RA(0,2)=3;
+    RA(1,0)= 4; RA(1,1)=5; RA(1,2)=6;
+    RA(2,0)= 7; RA(2,1)=8; RA(2,2)=9;
+    pass &= (mean_square(A-RA)<=TOL);
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(3);
+    b<<= 4,5,6;
+    B<<=fill_policy::sparse_insert(), 1, 2, 3, b, 7, project(b, range(1,3));
+    RB<<=1, 2, 3, 4, 5, 6, 7, 5, 6; // If the first worked we can now probably use it.
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(3);
+    b<<= 4,5,6;
+    B<<=fill_policy::sparse_insert(), move(1,0), b, fill_policy::sparse_insert(), move_to(0,0), 1, 2, 3, move(1,0), 7, project(b, range(1,3));
+    RB<<=1, 2, 3, 4, 5, 6, 7, 5, 6;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(3,3), RB(3,3);
+    vector<typename V::value_type>  b(9);
+    b<<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
+    B<<=b;
+    RB<<=1, 2, 3, 4, 5, 6, 7, 8, 9;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<=    2, 3,
+            4, 5;
+    B<<=fill_policy::sparse_insert(), C, C, C, C;
+    RB <<=   2,3,2,3,
+            4,5,4,5,
+            2,3,2,3,
+            4,5,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<=fill_policy::sparse_insert(), C, move_to(2,2), C;
+    RB<<=    2,3,0,0,
+            4,5,0,0,
+            0,0,2,3,
+            0,0,4,5;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    V C(2,2);
+    C <<= 2, 3, 4, 5;
+    B<<=fill_policy::sparse_insert(), move(1,1), C;
+    RB<<=    0,0,0,0,
+            0,2,3,0,
+            0,4,5,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_insert(), move_to(0,1), 2, 3, next_row(), 1, 2, next_row(), 4, 5;
+    RB<<=    0,2,3,0,
+            1,2,0,0,
+            4,5,0,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_insert(),traverse_policy::by_column(), move_to(0,1), 2, 3, 6, next_column(), 4, 5;
+    RB<<=    0,2,4,0,
+            0,3,5,0,
+            0,6,0,0,
+            0,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_insert(),traverse_policy::by_column(), move_to(0,1), 2, 3, next_row(), traverse_policy::by_row(), 4, 5;
+    RB<<=    0,2,0,0,
+            0,3,0,0,
+            0,0,0,0,
+            4,5,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_insert(),traverse_policy::by_column(), move_to(0,1), 2, 3, begin2(), traverse_policy::by_row(), 4, 5, 6, 7, 8;
+    RB<<=    0,2,0,0,
+            0,3,0,0,
+            4,5,6,7,
+            8,0,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    {
+    V B(4,4), RB(4,4);
+    B = zero_matrix<typename V::value_type>(4,4);
+    B<<=fill_policy::sparse_insert(),traverse_policy::by_column(), move_to(0,1), 2, 3, begin2(), traverse_policy::by_row(), 4, 5, 6, 7, 8,9, begin1(), 1, 2;
+    RB<<=    0,2,1,2,
+            0,3,0,0,
+            4,5,6,7,
+            8,9,0,0;
+    pass &= (mean_square(B-RB)<=TOL);
+    }
+
+    return pass;
+}
+
+
+BOOST_UBLAS_TEST_DEF (test_vector) {
+
+    BOOST_UBLAS_DEBUG_TRACE( "Starting operator \"<<= \" vector assignment tests" );
+
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<double,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<float,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<long,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<unsigned long,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<int,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<unsigned int,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<std::size_t,7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_vector<bounded_vector<char,7> >()));
+
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<unsigned int> >())
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<std::size_t> >())
+    BOOST_UBLAS_TEST_CHECK(test_vector<mapped_vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<compressed_vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<long> >())
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<unsigned long> >())
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector<coordinate_vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<compressed_vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_push_back<coordinate_vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<compressed_vector<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_vector_sparse_insert<coordinate_vector<char> >());
+}
+
+BOOST_UBLAS_TEST_DEF (test_matrix) {
+
+    BOOST_UBLAS_DEBUG_TRACE( "Starting operator \"<<= \" matrix assignment tests" );
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<matrix<char> >());
+
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<double,7, 7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<float,7, 7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<long,7, 7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<unsigned long,7, 7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<int,7,7 > >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<unsigned int,7, 7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<char,7, 7> >()));
+    BOOST_UBLAS_TEST_CHECK((test_matrix<bounded_matrix<std::size_t,7, 7> >()));
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<unsigned int> >())
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<std::size_t> >())
+    BOOST_UBLAS_TEST_CHECK(test_matrix<mapped_matrix<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<compressed_matrix<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<long> >())
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<unsigned long> >())
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix<coordinate_matrix<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<compressed_matrix<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_push_back<coordinate_matrix<char> >());
+
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<compressed_matrix<char> >());
+
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<double> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<float> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<unsigned long> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<unsigned int> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<std::size_t> >());
+    BOOST_UBLAS_TEST_CHECK(test_matrix_sparse_insert<coordinate_matrix<char> >());
+}
+
+
+int main () {
+    BOOST_UBLAS_TEST_BEGIN();
+
+    BOOST_UBLAS_TEST_DO( test_vector );
+    BOOST_UBLAS_TEST_DO( test_matrix );
+
+    BOOST_UBLAS_TEST_END();
+
+    return EXIT_SUCCESS;;
+}