From 4e887b9948dda5ab5f7c2ca06310d71539308b56 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B6rg=20Walter?= <jhr.walter@t-online.de>
Date: Tue, 19 Nov 2002 08:32:55 +0000
Subject: [PATCH] Bugfixes and compatibility hacks.

---
 doc/index.htm                                 |    8 +
 doc/matrix_sparse.htm                         | 1298 ++++++++++++++++-
 doc/overview.htm                              |   92 +-
 doc/vector.htm                                |    8 +-
 doc/vector_sparse.htm                         |  923 +++++++++++-
 include/boost/numeric/ublas/banded.hpp        |   19 +-
 include/boost/numeric/ublas/concepts.hpp      |   14 +-
 include/boost/numeric/ublas/config.hpp        |  229 +--
 include/boost/numeric/ublas/exception.hpp     |   49 +-
 include/boost/numeric/ublas/hermitian.hpp     |   21 +-
 include/boost/numeric/ublas/matrix_assign.hpp |   40 +
 .../boost/numeric/ublas/matrix_expression.hpp |  237 ++-
 include/boost/numeric/ublas/matrix_proxy.hpp  |  280 ++--
 include/boost/numeric/ublas/matrix_sparse.hpp |   92 +-
 .../boost/numeric/ublas/storage_sparse.hpp    |   22 +-
 include/boost/numeric/ublas/symmetric.hpp     |   19 +-
 include/boost/numeric/ublas/triangular.hpp    |   48 +-
 include/boost/numeric/ublas/vector_assign.hpp |   11 +
 .../boost/numeric/ublas/vector_expression.hpp |  160 +-
 include/boost/numeric/ublas/vector_proxy.hpp  |  139 +-
 include/boost/numeric/ublas/vector_sparse.hpp |   62 +-
 test1/test1.dsp                               |    4 -
 test2/test2.dsp                               |    4 -
 test2/test2.hpp                               |    2 +
 test3/Jamfile                                 |    2 +-
 test3/test3.dsp                               |    4 -
 test4/test4.dsp                               |    4 -
 test5/test5.dsp                               |    4 -
 test6/test6.dsp                               |    4 -
 29 files changed, 3123 insertions(+), 676 deletions(-)

diff --git a/doc/index.htm b/doc/index.htm
index 9a70499d..d6d3b28a 100644
--- a/doc/index.htm
+++ b/doc/index.htm
@@ -103,6 +103,10 @@ href="http://www.netlib.org/clapack">CLAPACK</a> test suite.</p>
         type="circle">
             <li><a href="vector_sparse.htm#sparse_vector">Sparse
                 Vector</a></li>
+            <li><a href="vector_sparse.htm#compressed_vector">Compressed
+                Vector</a></li>
+            <li><a href="vector_sparse.htm#coordinate_vector">Coordinate
+                Vector</a></li>
         </ul>
     </li>
     <li><a href="vector_proxy.htm">Vector Proxies</a><ul
@@ -166,6 +170,10 @@ href="http://www.netlib.org/clapack">CLAPACK</a> test suite.</p>
         type="circle">
             <li><a href="matrix_sparse.htm#sparse_matrix">Sparse
                 Matrix</a></li>
+            <li><a href="matrix_sparse.htm#compressed_matrix">Compressed
+                Matrix</a></li>
+            <li><a href="matrix_sparse.htm#coordinate_matrix">Coordinate
+                Matrix</a></li>
         </ul>
     </li>
     <li><a href="matrix_proxy.htm">Matrix Proxies</a><ul
diff --git a/doc/matrix_sparse.htm b/doc/matrix_sparse.htm
index 3fa5dca3..ca05ce26 100644
--- a/doc/matrix_sparse.htm
+++ b/doc/matrix_sparse.htm
@@ -19,10 +19,20 @@ width="277" height="86">Sparse Matrix</h1>
 <p>The templated class <code>sparse_matrix&lt;T, F, A&gt; </code>is
 the base container adaptor for sparse matrices. For a <em>(m x n</em>)-dimensional
 sparse matrix and <em>0 &lt;= i &lt; m</em>,<em> 0 &lt;= j &lt; n</em>
-every non-zero element <em>m</em><sub><em>i, j</em></sub> is
-mapped to the <em>(i x n + j)-</em>th element of the container
-for row major orientation or the <em>(i + j x m)-</em>th element
-of the container for column major orientation.</p>
+the non-zero elements <em>m</em><sub><em>i, j</em></sub> are
+mapped via <em>(i x n + j)</em> for row major orientation or via <em>(i
++ j x m)</em> for column major orientation to consecutive
+elements of the associative container, i.e. for elements <em>k</em>
+= <em>m</em><sub><em>i</em></sub><sub><sub><em>1</em></sub></sub><sub><em>,j</em></sub><sub><sub><em>1</em></sub></sub><sub><em>
+</em></sub>and <em>k + 1 = m</em><sub><em>i</em></sub><sub><sub><em>2</em></sub></sub><sub><em>,j</em></sub><sub><sub><em>2</em></sub></sub><sub><em>
+</em></sub>of the container holds <em>i</em><sub><em>1</em></sub><em>
+&lt; i</em><sub><em>2</em></sub> or <em>(i</em><sub><em>1</em></sub><em>
+= i</em><sub><em>2 </em></sub>and <em>j</em><sub><em>1</em></sub><em>
+&lt; j</em><sub><em>2</em></sub><em>)</em> with row major
+orientation or <em>j</em><sub><em>1</em></sub><em> &lt; j</em><sub><em>2</em></sub>
+or <em>(j</em><sub><em>1</em></sub><em> = j</em><sub><em>2 </em></sub>and
+<em>i</em><sub><em>1</em></sub><em> &lt; i</em><sub><em>2</em></sub><em>)</em>
+with column major orientation.</p>
 
 <h4>Example</h4>
 
@@ -54,12 +64,14 @@ of the container for column major orientation.</p>
     </tr>
     <tr>
         <td><code>F</code></td>
-        <td>Functor describing the storage organization. <a href="#sparse_matrix_1">[1]</a></td>
+        <td>Functor describing the storage organization. <a
+        href="#sparse_matrix_1">[1]</a></td>
         <td><code>row_major</code></td>
     </tr>
     <tr>
         <td><code>A</code></td>
-        <td>The type of the adapted array. <a href="#sparse_matrix_2">[2]</a></td>
+        <td>The type of the adapted array. <a
+        href="#sparse_matrix_2">[2]</a></td>
         <td><code>map_array&lt;std::size_t, T&gt;</code></td>
     </tr>
 </table>
@@ -309,12 +321,12 @@ href="container.htm#matrix">Matrix</a>.</p>
 
 <h4>Notes</h4>
 
-<p><a name="#sparse_matrix_1">[1]</a> Supported parameters for the storage organization are
-<code>row_major</code> and <code>column_major</code>.</p>
+<p><a name="#sparse_matrix_1">[1]</a> Supported parameters for
+the storage organization are <code>row_major</code> and <code>column_major</code>.</p>
 
-<p><a name="#sparse_matrix_2">[2]</a> Supported parameters for the adapted array are
-<code>map_array&lt;std::size_t, T&gt</code> and <code>std::map&lt;std::size_t, T&gt</code>.
-</p>
+<p><a name="#sparse_matrix_2">[2]</a> Supported parameters for
+the adapted array are <code>map_array&lt;std::size_t, T&gt;</code>
+and <code>std::map&lt;std::size_t, T&gt;</code>. </p>
 
 <h4>Interface</h4>
 
@@ -608,6 +620,1270 @@ href="container.htm#matrix">Matrix</a>.</p>
         reverse_iterator2 rend2 ();
     };</code></pre>
 
+<h2><a name="compressed_matrix"></a>Compressed Matrix</h2>
+
+<h4>Description</h4>
+
+<p>The templated class <code>compressed_matrix&lt;T, F, IB, IA,
+TA&gt; </code>is the base container adaptor for compressed
+matrices. For a <em>(m x n</em>)-dimensional compressed matrix
+and <em>0 &lt;= i &lt; m</em>,<em> 0 &lt;= j &lt; n</em> the
+non-zero elements <em>m</em><sub><em>i, j</em></sub> are mapped
+via <em>(i x n + j)</em> for row major orientation or via <em>(i
++ j x m)</em> for column major orientation to consecutive
+elements of the index and value containers, i.e. for elements <em>k</em>
+= <em>m</em><sub><em>i</em></sub><sub><sub><em>1</em></sub></sub><sub><em>,j</em></sub><sub><sub><em>1</em></sub></sub><sub><em>
+</em></sub>and <em>k + 1 = m</em><sub><em>i</em></sub><sub><sub><em>2</em></sub></sub><sub><em>,j</em></sub><sub><sub><em>2</em></sub></sub><sub><em>
+</em></sub>of the container holds <em>i</em><sub><em>1</em></sub><em>
+&lt; i</em><sub><em>2</em></sub> or <em>(i</em><sub><em>1</em></sub><em>
+= i</em><sub><em>2 </em></sub>and <em>j</em><sub><em>1</em></sub><em>
+&lt; j</em><sub><em>2</em></sub><em>)</em> with row major
+orientation or <em>j</em><sub><em>1</em></sub><em> &lt; j</em><sub><em>2</em></sub>
+or <em>(j</em><sub><em>1</em></sub><em> = j</em><sub><em>2 </em></sub>and
+<em>i</em><sub><em>1</em></sub><em> &lt; i</em><sub><em>2</em></sub><em>)</em>
+with column major orientation.</p>
+
+<h4>Example</h4>
+
+<pre>int main () {
+    using namespace boost::numeric::ublas;
+    compressed_matrix&lt;double&gt; m (3, 3, 3 * 3);
+    for (int i = 0; i &lt; m.size1 (); ++ i) 
+        for (int j = 0; j &lt; m.size2 (); ++ j) 
+            m (i, j) = 3 * i + j;
+    std::cout &lt;&lt; m &lt;&lt; std::endl;
+}</pre>
+
+<h4>Definition</h4>
+
+<p>Defined in the header matrix_sparse.hpp.</p>
+
+<h4>Template parameters</h4>
+
+<table border="1">
+    <tr>
+        <th>Parameter </th>
+        <th>Description </th>
+        <th>Default </th>
+    </tr>
+    <tr>
+        <td><code>T</code> </td>
+        <td>The type of object stored in the compressed matrix. </td>
+        <td>&nbsp;</td>
+    </tr>
+    <tr>
+        <td><code>F</code></td>
+        <td>Functor describing the storage organization. <a
+        href="#compressed_matrix_1">[1]</a></td>
+        <td><code>row_major</code></td>
+    </tr>
+    <tr>
+        <td><code>IB</code></td>
+        <td>The index base of the compressed vector. <a
+        href="#compressed_matrix_2">[2]</a></td>
+        <td><code>0</code></td>
+    </tr>
+    <tr>
+        <td><code>IA</code></td>
+        <td>The type of the adapted array for indices. <a
+        href="#compressed_matrix_3">[3]</a></td>
+        <td><code>unbounded_array&lt;std::size_t&gt;</code></td>
+    </tr>
+    <tr>
+        <td><code>TA</code></td>
+        <td>The type of the adapted array for values. <a
+        href="#compressed_matrix_3">[3]</a></td>
+        <td><code>unbounded_array&lt;T&gt;</code></td>
+    </tr>
+</table>
+
+<h4>Model of</h4>
+
+<p><a href="container.htm#matrix">Matrix</a>. </p>
+
+<h4>Type requirements</h4>
+
+<p>None, except for those imposed by the requirements of <a
+href="container.htm#matrix">Matrix</a>.</p>
+
+<h4>Public base classes</h4>
+
+<p><code>matrix_expression&lt;compressed_matrix&lt;T, F, IB, IA,
+TA&gt; &gt;</code> </p>
+
+<h4>Members</h4>
+
+<table border="1">
+    <tr>
+        <th>Member </th>
+        <th>Description </th>
+    </tr>
+    <tr>
+        <td><code>compressed_matrix ()</code> </td>
+        <td>Allocates a <code>compressed_matrix </code>that holds
+        at most zero rows of zero elements.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_matrix (size_type size1, size_type2,
+        size_type non_zeros)</code></td>
+        <td>Allocates a <code>compressed_matrix </code>that holds
+        at most <code>size1 </code>rows of <code>size2 </code>elements.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_matrix (const compressed_matrix
+        &amp;m)</code></td>
+        <td>The copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix (size_type non_zeros, const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>void resize (size_type size1, size_type size2,
+        size_type non_zeros)</code></td>
+        <td>Reallocates a <code>compressed_matrix </code>to hold
+        at most <code>size1 </code>rows of <code>size2 </code>elements.
+        The content of the <code>compressed_matrix </code>is
+        preserved.</td>
+    </tr>
+    <tr>
+        <td><code>size_type size1 () const</code></td>
+        <td>Returns the number of rows. </td>
+    </tr>
+    <tr>
+        <td><code>size_type size2 () const</code></td>
+        <td>Returns the number of columns. </td>
+    </tr>
+    <tr>
+        <td><code>const_reference operator () (size_type i,
+        size_type j) const</code></td>
+        <td>Returns the value of the <code>j</code>-th element in
+        the<code> i</code>-th row. </td>
+    </tr>
+    <tr>
+        <td><code>reference operator () (size_type i, size_type
+        j)</code></td>
+        <td>Returns a reference of the <code>j</code>-th element
+        in the<code> i</code>-th row. </td>
+    </tr>
+    <tr>
+        <td><code>compressed_matrix &amp;operator = (const
+        compressed_matrix &amp;m)</code></td>
+        <td>The assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_matrix &amp;assign_temporary
+        (compressed_matrix &amp;m)</code></td>
+        <td>Assigns a temporary. May change the compressed matrix
+        <code>m</code>.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix &amp;operator = (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix &amp;assign (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Assigns a matrix expression to the compressed matrix.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix &amp;operator += (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Adds the matrix
+        expression to the compressed matrix.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix &amp;plus_assign (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Adds a matrix expression to the compressed matrix.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix &amp;operator -= (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Subtracts the matrix
+        expression from the compressed matrix.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_matrix &amp;minus_assign (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Subtracts a matrix expression from the compressed
+        matrix. Left and right hand side of the assignment should
+        be independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        compressed_matrix &amp;operator *= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Multiplies the
+        compressed matrix with a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        compressed_matrix &amp;operator /= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Divides the
+        compressed matrix through a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>void swap (compressed_matrix &amp;m)</code></td>
+        <td>Swaps the contents of the compressed matrices. </td>
+    </tr>
+    <tr>
+        <td><code>void insert (size_type i, size_type j,
+        const_reference t)</code></td>
+        <td>Inserts the value <code>t</code> at the <code>j</code>-th
+        element of the <code>i</code>-th row.</td>
+    </tr>
+    <tr>
+        <td><code>void erase (size_type i, size_type j)</code></td>
+        <td>Erases the value at the <code>j</code>-th element of
+        the <code>i</code>-th row.</td>
+    </tr>
+    <tr>
+        <td><code>void clear ()</code></td>
+        <td>Clears the compressed matrix.</td>
+    </tr>
+    <tr>
+        <td><code>const_iterator1 begin1 () const</code></td>
+        <td>Returns a <code>const_iterator1</code> pointing to
+        the beginning of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator1 end1 () const</code></td>
+        <td>Returns a <code>const_iterator1</code> pointing to
+        the end of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator1 begin1 () </code></td>
+        <td>Returns a <code>iterator1</code> pointing to the
+        beginning of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator1 end1 () </code></td>
+        <td>Returns a <code>iterator1</code> pointing to the end
+        of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator2 begin2 () const</code></td>
+        <td>Returns a <code>const_iterator2</code> pointing to
+        the beginning of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator2 end2 () const</code></td>
+        <td>Returns a <code>const_iterator2</code> pointing to
+        the end of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator2 begin2 () </code></td>
+        <td>Returns a <code>iterator2</code> pointing to the
+        beginning of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator2 end2 () </code></td>
+        <td>Returns a <code>iterator2</code> pointing to the end
+        of the <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator1 rbegin1 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator1</code>
+        pointing to the beginning of the reversed <code>compressed_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator1 rend1 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator1</code>
+        pointing to the end of the reversed <code>compressed_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator1 rbegin1 () </code></td>
+        <td>Returns a <code>reverse_iterator1</code> pointing to
+        the beginning of the reversed <code>compressed_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator1 rend1 () </code></td>
+        <td>Returns a <code>reverse_iterator1</code> pointing to
+        the end of the reversed <code>compressed_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator2 rbegin2 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator2</code>
+        pointing to the beginning of the reversed <code>compressed_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator2 rend2 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator2</code>
+        pointing to the end of the reversed <code>compressed_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator2 rbegin2 () </code></td>
+        <td>Returns a <code>reverse_iterator2</code> pointing to
+        the beginning of the reversed <code>compressed_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator2 rend2 () </code></td>
+        <td>Returns a <code>reverse_iterator2</code> pointing to
+        the end of the reversed <code>compressed_matrix</code>. </td>
+    </tr>
+</table>
+
+<h4>Notes</h4>
+
+<p><a name="#compressed_matrix_1">[1]</a> Supported parameters
+for the storage organization are <code>row_major</code> and <code>column_major</code>.</p>
+
+<p><a name="#compressed_matrix_2">[2]</a> Supported parameters
+for the index base are <code>0 </code>and <code>1</code> at
+least. </p>
+
+<p><a name="#compressed_matrix_3">[3]</a> Supported parameters
+for the adapted array are <code>unbounded_array&lt;&gt;</code>, <code>bounded_array&lt;&gt;</code>
+and <code>std::vector&lt;&gt;</code>. </p>
+
+<h4>Interface</h4>
+
+<pre><code>    // Array based sparse matrix class 
+    template&lt;class T, class F, std::size_t IB, class IA, class TA&gt;
+    class compressed_matrix: 
+        public matrix_expression&lt;compressed_matrix&lt;T, F, IB, IA, TA&gt; &gt; {
+    public:      
+        typedef std::size_t size_type;
+        typedef std::ptrdiff_t difference_type;
+        typedef T value_type;
+        typedef const T &amp;const_reference;
+        typedef T &amp;reference;
+        typedef const T *const_pointer;
+        typedef T *pointer;
+        typedef F functor_type;
+        typedef IA index_array_type;
+        typedef TA value_array_type;
+        typedef const compressed_matrix&lt;T, F, IB, IA, TA&gt; const_self_type;
+        typedef compressed_matrix&lt;T, F, IB, IA, TA&gt; self_type;
+        typedef const matrix_const_reference&lt;const_self_type&gt; const_closure_type;
+        typedef matrix_reference&lt;self_type&gt; closure_type;
+        typedef typename A::const_iterator const_iterator_type;
+        typedef typename A::iterator iterator_type;
+        typedef sparse_tag storage_category;
+        typedef typename F::orientation_category orientation_category;
+
+        // Construction and destruction
+        compressed_matrix ();
+        compressed_matrix (size_type size1, size_type size2, size_type non_zeros = 0); 
+        compressed_matrix (const compressed_matrix &amp;m);
+        template&lt;class AE&gt;
+        compressed_matrix (const matrix_expression&lt;AE&gt; &amp;ae, size_type non_zeros = 0);
+
+        // Accessors
+        size_type size1 () const;
+        size_type size2 () const;
+        size_type non_zeros () const;
+        static size_type index_base ();
+        const index_array_type &amp;index1_data () const;
+        index_array_type &amp;index1_data ();
+        const index_array_type &amp;index2_data () const;
+        index_array_type &amp;index2_data ();
+        const value_array_type &amp;value_data () const;
+        value_array_type &amp;value_data ();
+
+        // Resizing
+        void resize (size_type size1, size_type size2, size_type non_zeros = 0);
+
+        // Element access
+        const_reference operator () (size_type i, size_type j) const;
+        reference operator () (size_type i, size_type j);
+
+        // Assignment
+        compressed_matrix &amp;operator = (const compressed_matrix &amp;m);
+        compressed_matrix &amp;assign_temporary (compressed_matrix &amp;m);
+        template&lt;class AE&gt;
+        compressed_matrix &amp;operator = (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_matrix &amp;reset (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_matrix &amp;assign (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_matrix&amp; operator += (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_matrix &amp;plus_assign (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_matrix&amp; operator -= (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_matrix &amp;minus_assign (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AT&gt;
+        compressed_matrix&amp; operator *= (const AT &amp;at);
+        template&lt;class AT&gt;
+        compressed_matrix&amp; operator /= (const AT &amp;at);
+
+        // Swapping
+        void swap (compressed_matrix &amp;m);
+        friend void swap (compressed_matrix &amp;m1, compressed_matrix &amp;m2);
+
+        // Element insertion and erasure
+        void insert (size_type i, size_type j, const_reference t);
+        void erase (size_type i, size_type j);
+        void clear ();
+
+        class const_iterator1;
+        class iterator1;
+        class const_iterator2;
+        class iterator2;
+        typedef reverse_iterator_base1&lt;const_iterator1&gt; const_reverse_iterator1;
+        typedef reverse_iterator_base1&lt;iterator1&gt; reverse_iterator1;
+        typedef reverse_iterator_base2&lt;const_iterator2&gt; const_reverse_iterator2;
+        typedef reverse_iterator_base2&lt;iterator2&gt; reverse_iterator2;
+
+        // Element lookup
+        const_iterator1 find1 (int rank, size_type i, size_type j) const;
+        iterator1 find1 (int rank, size_type i, size_type j);
+        const_iterator2 find2 (int rank, size_type i, size_type j) const;
+        iterator2 find2 (int rank, size_type i, size_type j);
+        const_iterator1 find_first1 (int rank, size_type i, size_type j) const;
+        iterator1 find_first1 (int rank, size_type i, size_type j);
+        const_iterator1 find_last1 (int rank, size_type i, size_type j) const;
+        iterator1 find_last1 (int rank, size_type i, size_type j);
+        const_iterator2 find_first2 (int rank, size_type i, size_type j) const;
+        iterator2 find_first2 (int rank, size_type i, size_type j);
+        const_iterator2 find_last2 (int rank, size_type i, size_type j) const;
+        iterator2 find_last2 (int rank, size_type i, size_type j);
+
+        // Iterators simply are pointers.
+
+        class const_iterator1:
+            public container_const_reference&lt;compressed_matrix&gt;,
+            public bidirectional_iterator_base&lt;const_iterator1, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename compressed_matrix::difference_type difference_type;
+            typedef typename compressed_matrix::value_type value_type;
+            typedef typename compressed_matrix::const_reference reference;
+            typedef typename compressed_matrix::const_pointer pointer;
+            typedef const_iterator2 dual_iterator_type;
+            typedef const_reverse_iterator2 dual_reverse_iterator_type;
+            typedef typename functor_type::functor1_type functor1_type;
+
+            // Construction and destruction
+            const_iterator1 ();
+            const_iterator1 (const compressed_matrix &amp;m, int rank, size_type i, size_type j, const const_iterator_type &amp;it);
+            const_iterator1 (const iterator1 &amp;it);
+
+            // Arithmetic
+            const_iterator1 &amp;operator ++ ();
+            const_iterator1 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            const_iterator2 begin () const;
+            const_iterator2 end () const;
+            const_reverse_iterator2 rbegin () const;
+            const_reverse_iterator2 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            const_iterator1 &amp;operator = (const const_iterator1 &amp;it);
+
+            // Comparison
+            bool operator == (const const_iterator1 &amp;it) const;
+        };
+
+        const_iterator1 begin1 () const;
+        const_iterator1 end1 () const;
+
+        class iterator1:
+            public container_reference&lt;compressed_matrix&gt;,
+            public bidirectional_iterator_base&lt;iterator1, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename compressed_matrix::difference_type difference_type;
+            typedef typename compressed_matrix::value_type value_type;
+            typedef typename compressed_matrix::reference reference;
+            typedef typename compressed_matrix::pointer pointer;
+            typedef iterator2 dual_iterator_type;
+            typedef reverse_iterator2 dual_reverse_iterator_type;
+            typedef typename functor_type::functor1_type functor1_type;
+
+            // Construction and destruction
+            iterator1 ();
+            iterator1 (compressed_matrix &amp;m, int rank, size_type i, size_type j, const iterator_type &amp;it);
+
+            // Arithmetic
+            iterator1 &amp;operator ++ ();
+            iterator1 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            iterator2 begin () const;
+            iterator2 end () const;
+            reverse_iterator2 rbegin () const;
+            reverse_iterator2 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            iterator1 &amp;operator = (const iterator1 &amp;it);
+
+            // Comparison
+            bool operator == (const iterator1 &amp;it) const;
+        };
+
+        iterator1 begin1 ();
+        iterator1 end1 ();
+
+        class const_iterator2:
+            public container_const_reference&lt;compressed_matrix&gt;,
+            public bidirectional_iterator_base&lt;const_iterator2, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename compressed_matrix::difference_type difference_type;
+            typedef typename compressed_matrix::value_type value_type;
+            typedef typename compressed_matrix::const_reference reference;
+            typedef typename compressed_matrix::const_pointer pointer;
+            typedef const_iterator1 dual_iterator_type;
+            typedef const_reverse_iterator1 dual_reverse_iterator_type;
+            typedef typename functor_type::functor2_type functor2_type;
+
+            // Construction and destruction
+            const_iterator2 ();
+            const_iterator2 (const compressed_matrix &amp;m, int rank, size_type i, size_type j, const const_iterator_type &amp;it);
+            const_iterator2 (const iterator2 &amp;it);
+
+            // Arithmetic
+            const_iterator2 &amp;operator ++ ();
+            const_iterator2 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            const_iterator1 begin () const;
+            const_iterator1 end () const;
+            const_reverse_iterator1 rbegin () const;
+            const_reverse_iterator1 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            const_iterator2 &amp;operator = (const const_iterator2 &amp;it);
+
+            // Comparison
+            bool operator == (const const_iterator2 &amp;it) const;
+        };
+
+        const_iterator2 begin2 () const;
+        const_iterator2 end2 () const;
+
+        class iterator2:
+            public container_reference&lt;compressed_matrix&gt;,
+            public bidirectional_iterator_base&lt;iterator2, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename compressed_matrix::difference_type difference_type;
+            typedef typename compressed_matrix::value_type value_type;
+            typedef typename compressed_matrix::reference reference;
+            typedef typename compressed_matrix::pointer pointer;
+            typedef iterator1 dual_iterator_type;
+            typedef reverse_iterator1 dual_reverse_iterator_type;
+            typedef typename functor_type::functor2_type functor2_type;
+
+            // Construction and destruction
+            iterator2 ();
+            iterator2 (compressed_matrix &amp;m, int rank, size_type i, size_type j, const iterator_type &amp;it);
+
+            // Arithmetic
+            iterator2 &amp;operator ++ ();
+            iterator2 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            iterator1 begin () const;
+            iterator1 end () const;
+            reverse_iterator1 rbegin () const;
+            reverse_iterator1 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            iterator2 &amp;operator = (const iterator2 &amp;it);
+
+            // Comparison
+            bool operator == (const iterator2 &amp;it) const;
+        };
+
+        iterator2 begin2 ();
+        iterator2 end2 ();
+
+        // Reverse iterators
+
+        const_reverse_iterator1 rbegin1 () const;
+        const_reverse_iterator1 rend1 () const;
+
+        reverse_iterator1 rbegin1 ();
+        reverse_iterator1 rend1 ();
+
+        const_reverse_iterator2 rbegin2 () const;
+        const_reverse_iterator2 rend2 () const;
+
+        reverse_iterator2 rbegin2 ();
+        reverse_iterator2 rend2 ();
+    };</code></pre>
+
+<h2><a name="coordinate_matrix"></a>Coordinate Matrix</h2>
+
+<h4>Description</h4>
+
+<p>The templated class <code>coordinate_matrix&lt;T, F, IB, IA,
+TA&gt; </code>is the base container adaptor for compressed
+matrices. For a <em>(m x n</em>)-dimensional sorted coordinate
+matrix and <em>0 &lt;= i &lt; m</em>,<em> 0 &lt;= j &lt; n</em>
+the non-zero elements <em>m</em><sub><em>i, j</em></sub> are
+mapped via <em>(i x n + j)</em> for row major orientation or via <em>(i
++ j x m)</em> for column major orientation to consecutive
+elements of the index and value containers, i.e. for elements <em>k</em>
+= <em>m</em><sub><em>i</em></sub><sub><sub><em>1</em></sub></sub><sub><em>,j</em></sub><sub><sub><em>1</em></sub></sub><sub><em>
+</em></sub>and <em>k + 1 = m</em><sub><em>i</em></sub><sub><sub><em>2</em></sub></sub><sub><em>,j</em></sub><sub><sub><em>2</em></sub></sub><sub><em>
+</em></sub>of the container holds <em>i</em><sub><em>1</em></sub><em>
+&lt; i</em><sub><em>2</em></sub> or <em>(i</em><sub><em>1</em></sub><em>
+= i</em><sub><em>2 </em></sub>and <em>j</em><sub><em>1</em></sub><em>
+&lt; j</em><sub><em>2</em></sub><em>)</em> with row major
+orientation or <em>j</em><sub><em>1</em></sub><em> &lt; j</em><sub><em>2</em></sub>
+or <em>(j</em><sub><em>1</em></sub><em> = j</em><sub><em>2 </em></sub>and
+<em>i</em><sub><em>1</em></sub><em> &lt; i</em><sub><em>2</em></sub><em>)</em>
+with column major orientation.</p>
+
+<h4>Example</h4>
+
+<pre>int main () {
+    using namespace boost::numeric::ublas;
+    coordinate_matrix&lt;double&gt; m (3, 3, 3 * 3);
+    for (int i = 0; i &lt; m.size1 (); ++ i) 
+        for (int j = 0; j &lt; m.size2 (); ++ j) 
+            m (i, j) = 3 * i + j;
+    std::cout &lt;&lt; m &lt;&lt; std::endl;
+}</pre>
+
+<h4>Definition</h4>
+
+<p>Defined in the header matrix_sparse.hpp.</p>
+
+<h4>Template parameters</h4>
+
+<table border="1">
+    <tr>
+        <th>Parameter </th>
+        <th>Description </th>
+        <th>Default </th>
+    </tr>
+    <tr>
+        <td><code>T</code> </td>
+        <td>The type of object stored in the coordinate matrix. </td>
+        <td>&nbsp;</td>
+    </tr>
+    <tr>
+        <td><code>F</code></td>
+        <td>Functor describing the storage organization. <a
+        href="#coordinate_matrix_1">[1]</a></td>
+        <td><code>row_major</code></td>
+    </tr>
+    <tr>
+        <td><code>IB</code></td>
+        <td>The index base of the coordinate vector. <a
+        href="#coordinate_matrix_2">[2]</a></td>
+        <td><code>0</code></td>
+    </tr>
+    <tr>
+        <td><code>IA</code></td>
+        <td>The type of the adapted array for indices. <a
+        href="#coordinate_matrix_3">[3]</a></td>
+        <td><code>unbounded_array&lt;std::size_t&gt;</code></td>
+    </tr>
+    <tr>
+        <td><code>TA</code></td>
+        <td>The type of the adapted array for values. <a
+        href="#coordinate_matrix_3">[3]</a></td>
+        <td><code>unbounded_array&lt;T&gt;</code></td>
+    </tr>
+</table>
+
+<h4>Model of</h4>
+
+<p><a href="container.htm#matrix">Matrix</a>. </p>
+
+<h4>Type requirements</h4>
+
+<p>None, except for those imposed by the requirements of <a
+href="container.htm#matrix">Matrix</a>.</p>
+
+<h4>Public base classes</h4>
+
+<p><code>matrix_expression&lt;coordinate_matrix&lt;T, F, IB, IA,
+TA&gt; &gt;</code> </p>
+
+<h4>Members</h4>
+
+<table border="1">
+    <tr>
+        <th>Member </th>
+        <th>Description </th>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix ()</code> </td>
+        <td>Allocates a <code>coordinate_matrix </code>that holds
+        at most zero rows of zero elements.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix (size_type size1, size_type2,
+        size_type non_zeros)</code></td>
+        <td>Allocates a <code>coordinate_matrix </code>that holds
+        at most <code>size1 </code>rows of <code>size2 </code>elements.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix (const coordinate_matrix
+        &amp;m)</code></td>
+        <td>The copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix (size_type non_zeros, const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>void resize (size_type size1, size_type size2,
+        size_type non_zeros)</code></td>
+        <td>Reallocates a <code>coordinate_matrix </code>to hold
+        at most <code>size1 </code>rows of <code>size2 </code>elements.
+        The content of the <code>coordinate_matrix </code>is
+        preserved.</td>
+    </tr>
+    <tr>
+        <td><code>size_type size1 () const</code></td>
+        <td>Returns the number of rows. </td>
+    </tr>
+    <tr>
+        <td><code>size_type size2 () const</code></td>
+        <td>Returns the number of columns. </td>
+    </tr>
+    <tr>
+        <td><code>const_reference operator () (size_type i,
+        size_type j) const</code></td>
+        <td>Returns the value of the <code>j</code>-th element in
+        the<code> i</code>-th row. </td>
+    </tr>
+    <tr>
+        <td><code>reference operator () (size_type i, size_type
+        j)</code></td>
+        <td>Returns a reference of the <code>j</code>-th element
+        in the<code> i</code>-th row. </td>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix &amp;operator = (const
+        coordinate_matrix &amp;m)</code></td>
+        <td>The assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix &amp;assign_temporary
+        (coordinate_matrix &amp;m)</code></td>
+        <td>Assigns a temporary. May change the coordinate matrix
+        <code>m</code>.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix &amp;operator = (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix &amp;assign (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Assigns a matrix expression to the coordinate matrix.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix &amp;operator += (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Adds the matrix
+        expression to the coordinate matrix.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix &amp;plus_assign (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Adds a matrix expression to the coordinate matrix.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix &amp;operator -= (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Subtracts the matrix
+        expression from the coordinate matrix.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_matrix &amp;minus_assign (const
+        matrix_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Subtracts a matrix expression from the coordinate
+        matrix. Left and right hand side of the assignment should
+        be independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        coordinate_matrix &amp;operator *= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Multiplies the
+        coordinate matrix with a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        coordinate_matrix &amp;operator /= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Divides the
+        coordinate matrix through a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>void swap (coordinate_matrix &amp;m)</code></td>
+        <td>Swaps the contents of the coordinate matrices. </td>
+    </tr>
+    <tr>
+        <td><code>void insert (size_type i, size_type j,
+        const_reference t)</code></td>
+        <td>Inserts the value <code>t</code> at the <code>j</code>-th
+        element of the <code>i</code>-th row.</td>
+    </tr>
+    <tr>
+        <td><code>void erase (size_type i, size_type j)</code></td>
+        <td>Erases the value at the <code>j</code>-th element of
+        the <code>i</code>-th row.</td>
+    </tr>
+    <tr>
+        <td><code>void clear ()</code></td>
+        <td>Clears the coordinate matrix.</td>
+    </tr>
+    <tr>
+        <td><code>const_iterator1 begin1 () const</code></td>
+        <td>Returns a <code>const_iterator1</code> pointing to
+        the beginning of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator1 end1 () const</code></td>
+        <td>Returns a <code>const_iterator1</code> pointing to
+        the end of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator1 begin1 () </code></td>
+        <td>Returns a <code>iterator1</code> pointing to the
+        beginning of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator1 end1 () </code></td>
+        <td>Returns a <code>iterator1</code> pointing to the end
+        of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator2 begin2 () const</code></td>
+        <td>Returns a <code>const_iterator2</code> pointing to
+        the beginning of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator2 end2 () const</code></td>
+        <td>Returns a <code>const_iterator2</code> pointing to
+        the end of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator2 begin2 () </code></td>
+        <td>Returns a <code>iterator2</code> pointing to the
+        beginning of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator2 end2 () </code></td>
+        <td>Returns a <code>iterator2</code> pointing to the end
+        of the <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator1 rbegin1 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator1</code>
+        pointing to the beginning of the reversed <code>coordinate_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator1 rend1 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator1</code>
+        pointing to the end of the reversed <code>coordinate_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator1 rbegin1 () </code></td>
+        <td>Returns a <code>reverse_iterator1</code> pointing to
+        the beginning of the reversed <code>coordinate_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator1 rend1 () </code></td>
+        <td>Returns a <code>reverse_iterator1</code> pointing to
+        the end of the reversed <code>coordinate_matrix</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator2 rbegin2 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator2</code>
+        pointing to the beginning of the reversed <code>coordinate_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator2 rend2 () const</code></td>
+        <td>Returns a <code>const_reverse_iterator2</code>
+        pointing to the end of the reversed <code>coordinate_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator2 rbegin2 () </code></td>
+        <td>Returns a <code>reverse_iterator2</code> pointing to
+        the beginning of the reversed <code>coordinate_matrix</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator2 rend2 () </code></td>
+        <td>Returns a <code>reverse_iterator2</code> pointing to
+        the end of the reversed <code>coordinate_matrix</code>. </td>
+    </tr>
+</table>
+
+<h4>Notes</h4>
+
+<p><a name="#coordinate_matrix_1">[1]</a> Supported parameters
+for the storage organization are <code>row_major</code> and <code>column_major</code>.</p>
+
+<p><a name="#coordinate_matrix_2">[2]</a> Supported parameters
+for the index base are <code>0 </code>and <code>1</code> at
+least. </p>
+
+<p><a name="#coordinate_matrix_3">[3]</a> Supported parameters
+for the adapted array are <code>unbounded_array&lt;&gt;</code>, <code>bounded_array&lt;&gt;</code>
+and <code>std::vector&lt;&gt;</code>. </p>
+
+<h4>Interface</h4>
+
+<pre><code>    // Array based sparse matrix class 
+    template&lt;class T, class F, std::size_t IB, class IA, class TA&gt;
+    class coordinate_matrix: 
+        public matrix_expression&lt;coordinate_matrix&lt;T, F, IB, IA, TA&gt; &gt; {
+    public:      
+        typedef std::size_t size_type;
+        typedef std::ptrdiff_t difference_type;
+        typedef T value_type;
+        typedef const T &amp;const_reference;
+        typedef T &amp;reference;
+        typedef const T *const_pointer;
+        typedef T *pointer;
+        typedef F functor_type;
+        typedef IA index_array_type;
+        typedef TA value_array_type;
+        typedef const coordinate_matrix&lt;T, F, IB, IA, TA&gt; const_self_type;
+        typedef coordinate_matrix&lt;T, F, IB, IA, TA&gt; self_type;
+        typedef const matrix_const_reference&lt;const_self_type&gt; const_closure_type;
+        typedef matrix_reference&lt;self_type&gt; closure_type;
+        typedef typename A::const_iterator const_iterator_type;
+        typedef typename A::iterator iterator_type;
+        typedef sparse_tag storage_category;
+        typedef typename F::orientation_category orientation_category;
+
+        // Construction and destruction
+        coordinate_matrix ();
+        coordinate_matrix (size_type size1, size_type size2, size_type non_zeros = 0); 
+        coordinate_matrix (const coordinate_matrix &amp;m);
+        template&lt;class AE&gt;
+        coordinate_matrix (const matrix_expression&lt;AE&gt; &amp;ae, size_type non_zeros = 0);
+
+        // Accessors
+        size_type size1 () const;
+        size_type size2 () const;
+        size_type non_zeros () const;
+        static size_type index_base ();
+        const index_array_type &amp;index1_data () const;
+        index_array_type &amp;index1_data ();
+        const index_array_type &amp;index2_data () const;
+        index_array_type &amp;index2_data ();
+        const value_array_type &amp;value_data () const;
+        value_array_type &amp;value_data ();
+
+        // Resizing
+        void resize (size_type size1, size_type size2, size_type non_zeros = 0);
+
+        // Element access
+        const_reference operator () (size_type i, size_type j) const;
+        reference operator () (size_type i, size_type j);
+
+        // Assignment
+        coordinate_matrix &amp;operator = (const coordinate_matrix &amp;m);
+        coordinate_matrix &amp;assign_temporary (coordinate_matrix &amp;m);
+        template&lt;class AE&gt;
+        coordinate_matrix &amp;operator = (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_matrix &amp;reset (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_matrix &amp;assign (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_matrix&amp; operator += (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_matrix &amp;plus_assign (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_matrix&amp; operator -= (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_matrix &amp;minus_assign (const matrix_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AT&gt;
+        coordinate_matrix&amp; operator *= (const AT &amp;at);
+        template&lt;class AT&gt;
+        coordinate_matrix&amp; operator /= (const AT &amp;at);
+
+        // Swapping
+        void swap (coordinate_matrix &amp;m);
+        friend void swap (coordinate_matrix &amp;m1, coordinate_matrix &amp;m2);
+
+        // Element insertion and erasure
+        void insert (size_type i, size_type j, const_reference t);
+        void erase (size_type i, size_type j);
+        void clear ();
+
+        class const_iterator1;
+        class iterator1;
+        class const_iterator2;
+        class iterator2;
+        typedef reverse_iterator_base1&lt;const_iterator1&gt; const_reverse_iterator1;
+        typedef reverse_iterator_base1&lt;iterator1&gt; reverse_iterator1;
+        typedef reverse_iterator_base2&lt;const_iterator2&gt; const_reverse_iterator2;
+        typedef reverse_iterator_base2&lt;iterator2&gt; reverse_iterator2;
+
+        // Element lookup
+        const_iterator1 find1 (int rank, size_type i, size_type j) const;
+        iterator1 find1 (int rank, size_type i, size_type j);
+        const_iterator2 find2 (int rank, size_type i, size_type j) const;
+        iterator2 find2 (int rank, size_type i, size_type j);
+        const_iterator1 find_first1 (int rank, size_type i, size_type j) const;
+        iterator1 find_first1 (int rank, size_type i, size_type j);
+        const_iterator1 find_last1 (int rank, size_type i, size_type j) const;
+        iterator1 find_last1 (int rank, size_type i, size_type j);
+        const_iterator2 find_first2 (int rank, size_type i, size_type j) const;
+        iterator2 find_first2 (int rank, size_type i, size_type j);
+        const_iterator2 find_last2 (int rank, size_type i, size_type j) const;
+        iterator2 find_last2 (int rank, size_type i, size_type j);
+
+        // Iterators simply are pointers.
+
+        class const_iterator1:
+            public container_const_reference&lt;coordinate_matrix&gt;,
+            public bidirectional_iterator_base&lt;const_iterator1, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename coordinate_matrix::difference_type difference_type;
+            typedef typename coordinate_matrix::value_type value_type;
+            typedef typename coordinate_matrix::const_reference reference;
+            typedef typename coordinate_matrix::const_pointer pointer;
+            typedef const_iterator2 dual_iterator_type;
+            typedef const_reverse_iterator2 dual_reverse_iterator_type;
+            typedef typename functor_type::functor1_type functor1_type;
+
+            // Construction and destruction
+            const_iterator1 ();
+            const_iterator1 (const coordinate_matrix &amp;m, int rank, size_type i, size_type j, const const_iterator_type &amp;it);
+            const_iterator1 (const iterator1 &amp;it);
+
+            // Arithmetic
+            const_iterator1 &amp;operator ++ ();
+            const_iterator1 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            const_iterator2 begin () const;
+            const_iterator2 end () const;
+            const_reverse_iterator2 rbegin () const;
+            const_reverse_iterator2 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            const_iterator1 &amp;operator = (const const_iterator1 &amp;it);
+
+            // Comparison
+            bool operator == (const const_iterator1 &amp;it) const;
+        };
+
+        const_iterator1 begin1 () const;
+        const_iterator1 end1 () const;
+
+        class iterator1:
+            public container_reference&lt;coordinate_matrix&gt;,
+            public bidirectional_iterator_base&lt;iterator1, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename coordinate_matrix::difference_type difference_type;
+            typedef typename coordinate_matrix::value_type value_type;
+            typedef typename coordinate_matrix::reference reference;
+            typedef typename coordinate_matrix::pointer pointer;
+            typedef iterator2 dual_iterator_type;
+            typedef reverse_iterator2 dual_reverse_iterator_type;
+            typedef typename functor_type::functor1_type functor1_type;
+
+            // Construction and destruction
+            iterator1 ();
+            iterator1 (coordinate_matrix &amp;m, int rank, size_type i, size_type j, const iterator_type &amp;it);
+
+            // Arithmetic
+            iterator1 &amp;operator ++ ();
+            iterator1 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            iterator2 begin () const;
+            iterator2 end () const;
+            reverse_iterator2 rbegin () const;
+            reverse_iterator2 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            iterator1 &amp;operator = (const iterator1 &amp;it);
+
+            // Comparison
+            bool operator == (const iterator1 &amp;it) const;
+        };
+
+        iterator1 begin1 ();
+        iterator1 end1 ();
+
+        class const_iterator2:
+            public container_const_reference&lt;coordinate_matrix&gt;,
+            public bidirectional_iterator_base&lt;const_iterator2, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename coordinate_matrix::difference_type difference_type;
+            typedef typename coordinate_matrix::value_type value_type;
+            typedef typename coordinate_matrix::const_reference reference;
+            typedef typename coordinate_matrix::const_pointer pointer;
+            typedef const_iterator1 dual_iterator_type;
+            typedef const_reverse_iterator1 dual_reverse_iterator_type;
+            typedef typename functor_type::functor2_type functor2_type;
+
+            // Construction and destruction
+            const_iterator2 ();
+            const_iterator2 (const coordinate_matrix &amp;m, int rank, size_type i, size_type j, const const_iterator_type &amp;it);
+            const_iterator2 (const iterator2 &amp;it);
+
+            // Arithmetic
+            const_iterator2 &amp;operator ++ ();
+            const_iterator2 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            const_iterator1 begin () const;
+            const_iterator1 end () const;
+            const_reverse_iterator1 rbegin () const;
+            const_reverse_iterator1 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            const_iterator2 &amp;operator = (const const_iterator2 &amp;it);
+
+            // Comparison
+            bool operator == (const const_iterator2 &amp;it) const;
+        };
+
+        const_iterator2 begin2 () const;
+        const_iterator2 end2 () const;
+
+        class iterator2:
+            public container_reference&lt;coordinate_matrix&gt;,
+            public bidirectional_iterator_base&lt;iterator2, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename coordinate_matrix::difference_type difference_type;
+            typedef typename coordinate_matrix::value_type value_type;
+            typedef typename coordinate_matrix::reference reference;
+            typedef typename coordinate_matrix::pointer pointer;
+            typedef iterator1 dual_iterator_type;
+            typedef reverse_iterator1 dual_reverse_iterator_type;
+            typedef typename functor_type::functor2_type functor2_type;
+
+            // Construction and destruction
+            iterator2 ();
+            iterator2 (coordinate_matrix &amp;m, int rank, size_type i, size_type j, const iterator_type &amp;it);
+
+            // Arithmetic
+            iterator2 &amp;operator ++ ();
+            iterator2 &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            iterator1 begin () const;
+            iterator1 end () const;
+            reverse_iterator1 rbegin () const;
+            reverse_iterator1 rend () const;
+
+            // Indices
+            size_type index1 () const;
+            size_type index2 () const;
+
+            // Assignment 
+            iterator2 &amp;operator = (const iterator2 &amp;it);
+
+            // Comparison
+            bool operator == (const iterator2 &amp;it) const;
+        };
+
+        iterator2 begin2 ();
+        iterator2 end2 ();
+
+        // Reverse iterators
+
+        const_reverse_iterator1 rbegin1 () const;
+        const_reverse_iterator1 rend1 () const;
+
+        reverse_iterator1 rbegin1 ();
+        reverse_iterator1 rend1 ();
+
+        const_reverse_iterator2 rbegin2 () const;
+        const_reverse_iterator2 rend2 () const;
+
+        reverse_iterator2 rbegin2 ();
+        reverse_iterator2 rend2 ();
+    };</code></pre>
+
 <hr>
 
 <p>Copyright (©) 2000-2002 Joerg Walter, Mathias Koch <br>
diff --git a/doc/overview.htm b/doc/overview.htm
index 062db06c..de8799f4 100644
--- a/doc/overview.htm
+++ b/doc/overview.htm
@@ -419,15 +419,25 @@ how BLAS calls may be mapped onto our classes. </p>
     </tr>
     <tr>
         <td><code>_t_sv</code></td>
-        <td><code>x = solve (A, x, tag) </code>or<br>
-        <code>x = solve (trans (A), x, tag) </code>or<code><br>
-        x = solve (herm (A), x, tag)</code></td>
-        <td><em>x &lt;- A</em><sup><em>-1</em></sup><em> x </em>or<em><br>
+        <td><code>y = solve (A, x, tag) </code>or<br>
+        <code>inplace_solve (A, x, tag) </code>or<br>
+        <code>y = solve (trans (A), x, tag) </code>or<code><br>
+        inplace_solve (trans (A), x, tag) </code>or<font
+        face="Courier New"><code><br>
+        </code></font><code>y = solve (herm (A), x, tag)</code>or<code><br>
+        inplace_solve (herm (A), x, tag)</code></td>
+        <td><em>y &lt;- A</em><sup><em>-1</em></sup><em> x </em>or<em><br>
+        x &lt;- A</em><sup><em>-1</em></sup><em> x </em>or<em><br>
+        y &lt;- A</em><sup><em>T</em></sup><sup><sup><em>-1</em></sup></sup><em>
+        x </em>or<em><br>
         x &lt;- A</em><sup><em>T</em></sup><sup><sup><em>-1</em></sup></sup><em>
         x </em>or<em><br>
+        y &lt;- A</em><sup><em>H</em></sup><sup><sup><em>-1</em></sup></sup><em>
+        x</em> or<em><br>
         x &lt;- A</em><sup><em>H</em></sup><sup><sup><em>-1</em></sup></sup><em>
         x</em></td>
-        <td>Solves a linear equation for a triangular matrix.</td>
+        <td>Solves a system of linear equations with triangular
+        form, i.e. <em>A </em>is triangular.</td>
     </tr>
     <tr>
         <td><code>_g_mv<br>
@@ -506,15 +516,23 @@ how BLAS calls may be mapped onto our classes. </p>
     </tr>
     <tr>
         <td><code>_t_sm</code></td>
-        <td><code>B = solve (A, B, tag) </code>or<br>
-        <code>B = solve (trans (A), B, tag) </code>or<code><br>
-        B = solve (herm (A), B, tag)</code></td>
-        <td><em>B &lt;- A</em><sup><em>-1</em></sup><em> B </em>or<em><br>
-        B &lt;- A</em><sup><em>T</em></sup><sup><sup><em>-1</em></sup></sup><em>
+        <td><code>C = solve (A, B, tag) </code>or<br>
+        <code>inplace_solve (A, B, tag) </code>or<br>
+        <code>C = solve (trans (A), B, tag) </code>or<code><br>
+        inplace_solve (trans (A), B, tag) </code>or<code><br>
+        C = solve (herm (A), B, tag)</code> or<code><br>
+        inplace_solve (herm (A), B, tag)</code></td>
+        <td><em>C &lt;- A</em><sup><em>-1</em></sup><em> B </em>or<em><br>
+        B &lt;- A</em><sup><em>-1</em></sup><em> B </em>or<em><br>
+        C &lt;- A</em><sup><em>T</em></sup><sup><sup><em>-1</em></sup></sup><em>
         B </em>or<em><br>
+        B &lt;- A</em><sup><em>-1</em></sup><em> B </em>or<em><br>
+        C &lt;- A</em><sup><em>H</em></sup><sup><sup><em>-1</em></sup></sup><em>
+        B</em> or<em><br>
         B &lt;- A</em><sup><em>H</em></sup><sup><sup><em>-1</em></sup></sup><em>
         B</em></td>
-        <td>Solves a linear equation for a triangular matrix.</td>
+        <td>Solves a system of linear equations with triangular
+        form, i.e. <em>A </em>is triangular.</td>
     </tr>
     <tr>
         <td><code>_g_mm<br>
@@ -634,6 +652,18 @@ most common constructions of vectors and matrices comes next.</p>
         <td>Constructs a sparse vector, storage is provided by a
         map array.</td>
     </tr>
+    <tr>
+        <td><code>compressed_vector&lt;T&gt; <br>
+        &nbsp;&nbsp;v (size, non_zeros)</code></td>
+        <td>Constructs a compressed vector.<br>
+        The storage layout usually is BLAS compliant.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_vector&lt;T&gt; <br>
+        &nbsp;&nbsp;v (size, non_zeros)</code></td>
+        <td>Constructs a coordinate vector.<br>
+        The storage layout usually is BLAS compliant.</td>
+    </tr>
     <tr>
         <td><code>vector_range&lt;V&gt; <br>
         &nbsp;&nbsp;vr (v, range)</code></td>
@@ -801,6 +831,38 @@ most common constructions of vectors and matrices comes next.</p>
         <td>Constructs a sparse matrix, orientation is column
         major, storage is provided by a map array.</td>
     </tr>
+    <tr>
+        <td><code>compressed_matrix&lt;T, <br>
+        &nbsp;row_major&gt; <br>
+        &nbsp;&nbsp;m (size1, size2, non_zeros)</code></td>
+        <td>Constructs a compressed matrix, orientation is row
+        major.<br>
+        The storage layout usually is BLAS compliant.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_matrix&lt;T, <br>
+        &nbsp;column_major&gt; <br>
+        &nbsp;&nbsp;m (size1, size2, non_zeros)</code></td>
+        <td>Constructs a compressed matrix, orientation is column
+        major.<br>
+        The storage layout usually is BLAS compliant.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix&lt;T, <br>
+        &nbsp;row_major&gt; <br>
+        &nbsp;&nbsp;m (size1, size2, non_zeros)</code></td>
+        <td>Constructs a coordinate matrix, orientation is row
+        major.<br>
+        The storage layout usually is BLAS compliant.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_matrix&lt;T, <br>
+        &nbsp;column_major&gt; <br>
+        &nbsp;&nbsp;m (size1, size2, non_zeros)</code></td>
+        <td>Constructs a coordinate matrix, orientation is column
+        major.<br>
+        The storage layout usually is BLAS compliant.</td>
+    </tr>
     <tr>
         <td><code>matrix_row&lt;M&gt; <br>
         &nbsp;&nbsp;mr (m, i)</code></td>
@@ -940,6 +1002,10 @@ config.hpp. But we suspect this optimization to be fragile.</p>
         type="circle">
             <li><a href="vector_sparse.htm#sparse_vector">Sparse
                 Vector</a></li>
+            <li><a href="vector_sparse.htm#compressed_vector">Compressed
+                Vector</a></li>
+            <li><a href="vector_sparse.htm#coordinate_vector">Coordinate
+                Vector</a></li>
         </ul>
     </li>
     <li><a href="vector_proxy.htm">Vector Proxies</a><ul
@@ -1003,6 +1069,10 @@ config.hpp. But we suspect this optimization to be fragile.</p>
         type="circle">
             <li><a href="matrix_sparse.htm#sparse_matrix">Sparse
                 Matrix</a></li>
+            <li><a href="matrix_sparse.htm#compressed_matrix">Compressed
+                Matrix</a></li>
+            <li><a href="matrix_sparse.htm#coordinate_matrix">Coordinate
+                Matrix</a></li>
         </ul>
     </li>
     <li><a href="matrix_proxy.htm">Matrix Proxies</a><ul
diff --git a/doc/vector.htm b/doc/vector.htm
index 86bae77a..7b5b1ce6 100644
--- a/doc/vector.htm
+++ b/doc/vector.htm
@@ -50,7 +50,7 @@ is mapped to the <em>i-</em>th element of the container.</p>
     </tr>
     <tr>
         <td><code>A</code></td>
-        <td>The type of the adapted array.</td>
+        <td>The type of the adapted array. <a href="#vector_1">[1]</a></td>
         <td><code>unbounded_array&lt;T&gt;</code></td>
     </tr>
 </table>
@@ -246,6 +246,12 @@ href="container.htm#vector">Vector</a>.</p>
     </tr>
 </table>
 
+<h4>Notes</h4>
+
+<p><a name="#vector_1">[1]</a> Supported parameters for the
+adapted array are <code>unbounded_array&lt;T&gt;</code>, <code>bounded_array&lt;T&gt;</code>
+and <code>std::vector&lt;T&gt;</code>. </p>
+
 <h4>Interface</h4>
 
 <pre><code>    // Array based vector class 
diff --git a/doc/vector_sparse.htm b/doc/vector_sparse.htm
index ea782e01..ceed0dbb 100644
--- a/doc/vector_sparse.htm
+++ b/doc/vector_sparse.htm
@@ -18,9 +18,13 @@ width="277" height="86">Sparse Vector</h1>
 
 <p>The templated class <code>sparse_vector&lt;T, A&gt; </code>is
 the base container adaptor for sparse vectors. For a <em>n</em>-dimensional
-sparse vector and <em>0 &lt;= i &lt; n </em>every non-zero
-element <em>v</em><sub><em>i</em></sub> is mapped to the <em>i-</em>th
-element of the container.</p>
+sparse vector and <em>0 &lt;= i &lt; n </em>the non-zero elements
+<em>v</em><sub><em>i</em></sub> are mapped to consecutive
+elements of the associative container, i.e. for elements <em>k</em>
+= <em>v</em><sub><em>i</em></sub><sub><sub><em>1</em></sub></sub><sub><em>
+</em></sub>and <em>k + 1 = v</em><sub><em>i</em></sub><sub><sub><em>2</em></sub></sub><sub><em>
+</em></sub>of the container holds <em>i</em><sub><em>1</em></sub><em>
+&lt; i</em><sub><em>2</em></sub>.</p>
 
 <h4>Example</h4>
 
@@ -49,14 +53,10 @@ element of the container.</p>
         <td>The type of object stored in the sparse vector. </td>
         <td>&nbsp;</td>
     </tr>
-    <tr>
-        <td><code>F</code></td>
-        <td>Functor describing the storage organization.</td>
-        <td><code>forward</code></td>
-    </tr>
     <tr>
         <td><code>A</code></td>
-        <td>The type of the adapted array.</td>
+        <td>The type of the adapted array. <a
+        href="#sparse_vector_1">[1]</a></td>
         <td><code>map_array&lt;std::size_t, T&gt;</code></td>
     </tr>
 </table>
@@ -260,6 +260,12 @@ href="container.htm#vector">Vector</a>.</p>
     </tr>
 </table>
 
+<h4>Notes</h4>
+
+<p><a name="#sparse_vector_1">[1]</a> Supported parameters for
+the adapted array are <code>map_array&lt;std::size_t, T&gt;</code>
+and <code>std::map&lt;std::size_t, T&gt;</code>. </p>
+
 <h4>Interface</h4>
 
 <pre><code>    // Array based sparse vector class 
@@ -274,7 +280,6 @@ href="container.htm#vector">Vector</a>.</p>
         typedef T &amp;reference;
         typedef const T *const_pointer;
         typedef T *pointer;
-        typedef F functor_type;
         typedef A array_type;
         typedef const A const_array_type;
         typedef const sparse_vector&lt;T, A&gt; const_self_type;
@@ -434,6 +439,904 @@ href="container.htm#vector">Vector</a>.</p>
         reverse_iterator rend ();
     };</code></pre>
 
+<h2><a name="compressed_vector"></a>Compressed Vector</h2>
+
+<h4>Description</h4>
+
+<p>The templated class <code>compressed_vector&lt;T, IB, IA,
+TA&gt; </code>is the base container adaptor for compressed
+vectors. For a <em>n</em>-dimensional compressed vector and <em>0
+&lt;= i &lt; n </em>the non-zero elements <em>v</em><sub><em>i</em></sub>
+are mapped to consecutive elements of the index and value
+container, i.e. for elements <em>k</em> = <em>v</em><sub><em>i</em></sub><sub><sub><em>1</em></sub></sub><sub><em>
+</em></sub>and <em>k + 1 = v</em><sub><em>i</em></sub><sub><sub><em>2</em></sub></sub><sub><em>
+</em></sub>of these containers holds <em>i</em><sub><em>1</em></sub><em>
+&lt; i</em><sub><em>2</em></sub>.</p>
+
+<h4>Example</h4>
+
+<pre>int main () {
+    using namespace boost::numeric::ublas;
+    compressed_vector&lt;double&gt; v (3, 3);
+    for (int i = 0; i &lt; v.size (); ++ i) 
+        v (i) = i;
+    std::cout &lt;&lt; v &lt;&lt; std::endl;
+}</pre>
+
+<h4>Definition</h4>
+
+<p>Defined in the header vector_sparse.hpp.</p>
+
+<h4>Template parameters</h4>
+
+<table border="1">
+    <tr>
+        <th>Parameter </th>
+        <th>Description </th>
+        <th>Default </th>
+    </tr>
+    <tr>
+        <td><code>T</code> </td>
+        <td>The type of object stored in the compressed vector. </td>
+        <td>&nbsp;</td>
+    </tr>
+    <tr>
+        <td><code>IB</code></td>
+        <td>The index base of the compressed vector. <a
+        href="#compressed_vector_1">[1]</a></td>
+        <td><code>0</code></td>
+    </tr>
+    <tr>
+        <td><code>IA</code></td>
+        <td>The type of the adapted array for indices. <a
+        href="#compressed_vector_2">[2]</a></td>
+        <td><code>unbounded_array&lt;std::size_t&gt;</code></td>
+    </tr>
+    <tr>
+        <td><code>TA</code></td>
+        <td>The type of the adapted array for values. <a
+        href="#compressed_vector_2">[2]</a></td>
+        <td><code>unbounded_array&lt;T&gt;</code></td>
+    </tr>
+</table>
+
+<h4>Model of</h4>
+
+<p><a href="container.htm#vector">Vector</a>. </p>
+
+<h4>Type requirements</h4>
+
+<p>None, except for those imposed by the requirements of <a
+href="container.htm#vector">Vector</a>.</p>
+
+<h4>Public base classes</h4>
+
+<p><code>vector_expression&lt;compressed_vector&lt;T, IB, IA,
+TA&gt; &gt;</code> </p>
+
+<h4>Members</h4>
+
+<table border="1">
+    <tr>
+        <th>Member </th>
+        <th>Description </th>
+    </tr>
+    <tr>
+        <td><code>compressed_vector ()</code> </td>
+        <td>Allocates a<code> compressed_vector </code>that holds
+        zero elements.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_vector (size_type size, size_type
+        non_zeros)</code></td>
+        <td>Allocates a<code> compressed_vector</code> that holds
+        at most <code>size</code> elements.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_vector (const compressed_vector
+        &amp;v)</code></td>
+        <td>The copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector (size_type non_zeros, const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>void resize (size_type size, size_type
+        non_zeros)</code></td>
+        <td>Reallocates a <code>compressed_vector </code>to hold
+        at most <code>size</code> elements. The content of the <code>compressed_vector
+        </code>is preserved.</td>
+    </tr>
+    <tr>
+        <td><code>size_type size () const</code></td>
+        <td>Returns the size of the <code>compressed_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reference operator () (size_type i) const</code></td>
+        <td>Returns the value of the <code>i</code>-th element. </td>
+    </tr>
+    <tr>
+        <td><code>reference operator () (size_type i)</code></td>
+        <td>Returns a reference of the <code>i</code>-th element.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reference operator [] (size_type i) const</code></td>
+        <td>Returns the value of the <code>i</code>-th element. </td>
+    </tr>
+    <tr>
+        <td><code>reference operator [] (size_type i)</code></td>
+        <td>Returns a reference of the <code>i</code>-th element.
+        </td>
+    </tr>
+    <tr>
+        <td><code>compressed_vector &amp;operator = (const
+        compressed_vector &amp;v)</code></td>
+        <td>The assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>compressed_vector &amp;assign_temporary
+        (compressed_vector &amp;v)</code></td>
+        <td>Assigns a temporary. May change the compressed vector
+        <code>v</code>.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector &amp;operator = (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector &amp;assign (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Assigns a vector expression to the compressed vector.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector &amp;operator += (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Adds the vector
+        expression to the compressed vector.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector &amp;plus_assign (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Adds a vector expression to the compressed vector.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector &amp;operator -= (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Subtracts the vector
+        expression from the compressed vector.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        compressed_vector &amp;minus_assign (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Subtracts a vector expression from the compressed
+        vector. Left and right hand side of the assignment should
+        be independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        compressed_vector &amp;operator *= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Multiplies the
+        compressed vector with a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        compressed_vector &amp;operator /= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Divides the
+        compressed vector through a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>void swap (compressed_vector &amp;v)</code></td>
+        <td>Swaps the contents of the compressed vectors. </td>
+    </tr>
+    <tr>
+        <td><code>void insert (size_type i, const_reference t)</code></td>
+        <td>Inserts the value <code>t</code> at the <code>i</code>-th
+        element.</td>
+    </tr>
+    <tr>
+        <td><code>void erase (size_type i)</code></td>
+        <td>Erases the value at the <code>i</code>-th element.</td>
+    </tr>
+    <tr>
+        <td><code>void clear ()</code></td>
+        <td>Clears the compressed vector.</td>
+    </tr>
+    <tr>
+        <td><code>const_iterator begin () const</code></td>
+        <td>Returns a <code>const_iterator</code> pointing to the
+        beginning of the <code>compressed_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator end () const</code></td>
+        <td>Returns a <code>const_iterator</code> pointing to the
+        end of the <code>compressed_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator begin () </code></td>
+        <td>Returns a <code>iterator</code> pointing to the
+        beginning of the <code>compressed_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator end () </code></td>
+        <td>Returns a <code>iterator</code> pointing to the end
+        of the <code>compressed_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator rbegin () const</code></td>
+        <td>Returns a <code>const_reverse_iterator</code>
+        pointing to the beginning of the reversed <code>compressed_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator rend () const</code></td>
+        <td>Returns a <code>const_reverse_iterator</code>
+        pointing to the end of the reversed <code>compressed_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator rbegin () </code></td>
+        <td>Returns a <code>reverse_iterator</code> pointing to
+        the beginning of the reversed <code>compressed_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator rend () </code></td>
+        <td>Returns a <code>reverse_iterator</code> pointing to
+        the end of the reversed <code>compressed_vector</code>. </td>
+    </tr>
+</table>
+
+<h4>Notes</h4>
+
+<p><a name="#compressed_vector_1">[1]</a> Supported parameters
+for the index base are <code>0 </code>and <code>1</code> at
+least. </p>
+
+<p><a name="#compressed_vector_2">[2]</a> Supported parameters
+for the adapted array are <code>unbounded_array&lt;&gt;</code>, <code>bounded_array&lt;&gt;</code>
+and <code>std::vector&lt;&gt;</code>. </p>
+
+<h4>Interface</h4>
+
+<pre><code>    // Array based sparse vector class 
+    template&lt;class T, std::size_t IB, class IA, class TA&gt;
+    class compressed_vector: 
+        public vector_expression&lt;compressed_vector&lt;T, IB, IA, TA&gt; &gt; {
+    public:      
+        typedef std::size_t size_type;
+        typedef std::ptrdiff_t difference_type;
+        typedef T value_type;
+        typedef const T &amp;const_reference;
+        typedef T &amp;reference;
+        typedef const T *const_pointer;
+        typedef T *pointer;
+        typedef IA index_array_type;
+        typedef TA value_array_type;
+        typedef const compressed_vector&lt;T, IB, IA, TA&gt; const_self_type;
+        typedef compressed_vector&lt;T, IB, IA, TA&gt; self_type;
+        typedef const vector_const_reference&lt;const_self_type&gt; const_closure_type;
+        typedef vector_reference&lt;self_type&gt; closure_type;
+        typedef typename IA::const_iterator const_iterator_type;
+        typedef typename IA::iterator iterator_type;
+        typedef sparse_tag storage_category;
+
+        // Construction and destruction
+        compressed_vector (); 
+        compressed_vector (size_type size, size_type non_zeros = 0);
+        compressed_vector (const compressed_vector &amp;v);
+        template&lt;class AE&gt;
+        compressed_vector (const vector_expression&lt;AE&gt; &amp;ae, size_type non_zeros = 0);
+
+        // Accessors
+        size_type size () const;
+        size_type non_zeros () const;
+        static size_type index_base ();
+        const index_array_type &amp;index_data () const;
+        index_array_type &amp;index_data ();
+        const value_array_type &amp;value_data () const;
+        value_array_type &amp;value_data ();
+
+        // Resizing
+        void resize (size_type size, size_type non_zeros = 0);
+
+        // Element access
+        const_reference operator () (size_type i) const;
+        reference operator () (size_type i);
+
+        const_reference operator [] (size_type i) const;
+        reference operator [] (size_type i);
+
+        // Assignment
+        compressed_vector &amp;operator = (const compressed_vector &amp;v);
+        compressed_vector &amp;assign_temporary (compressed_vector &amp;v);
+        template&lt;class AE&gt;
+        compressed_vector &amp;operator = (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_vector &amp;reset (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_vector &amp;assign (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_vector &amp;operator += (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_vector &amp;plus_assign (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_vector &amp;operator -= (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        compressed_vector &amp;minus_assign (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AT&gt;
+        compressed_vector &amp;operator *= (const AT &amp;at);
+        template&lt;class AT&gt;
+        compressed_vector &amp;operator /= (const AT &amp;at);
+
+        // Swapping
+        void swap (compressed_vector &amp;v);
+        friend void swap (compressed_vector &amp;v1, compressed_vector &amp;v2);
+
+        // Element insertion and erasure
+        void insert (size_type i, const_reference t);
+        void erase (size_type i);
+        void clear ();
+
+        class const_iterator;
+        class iterator;
+
+        // Element lookup
+        const_iterator find (size_type i) const;
+        iterator find (size_type i);
+        const_iterator find_first (size_type i) const;
+        iterator find_first (size_type i);
+        const_iterator find_last (size_type i) const;
+        iterator find_last (size_type i);
+
+        // Iterators simply are pointers.
+
+        class const_iterator:
+            public container_const_reference&lt;compressed_vector&gt;,
+            public bidirectional_iterator_base&lt;const_iterator, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename compressed_vector::difference_type difference_type;
+            typedef typename compressed_vector::value_type value_type;
+            typedef typename compressed_vector::const_reference reference;
+            typedef typename compressed_vector::const_pointer pointer;
+
+            // Construction and destruction
+            const_iterator ();
+            const_iterator (const compressed_vector &amp;v, const const_iterator_type &amp;it);
+            const_iterator (const iterator &amp;it);
+
+            // Arithmetic
+            const_iterator &amp;operator ++ ();
+            const_iterator &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            // Index
+            size_type index () const;
+
+            // Assignment
+            const_iterator &amp;operator = (const const_iterator &amp;it);
+
+            // Comparison
+            bool operator == (const const_iterator &amp;it) const;
+        };
+
+        const_iterator begin () const;
+        const_iterator end () const;
+
+        class iterator:
+            public container_reference&lt;compressed_vector&gt;,
+            public bidirectional_iterator_base&lt;iterator, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename compressed_vector::difference_type difference_type;
+            typedef typename compressed_vector::value_type value_type;
+            typedef typename compressed_vector::reference reference;
+            typedef typename compressed_vector::pointer pointer;
+
+            // Construction and destruction
+            iterator ();
+            iterator (compressed_vector &amp;v, const iterator_type &amp;it);
+
+            // Arithmetic
+            iterator &amp;operator ++ ();
+            iterator &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            // Index
+            size_type index () const;
+
+            // Assignment
+            iterator &amp;operator = (const iterator &amp;it);
+
+            // Comparison
+            bool operator == (const iterator &amp;it) const;
+        };
+
+        iterator begin ();
+        iterator end ();
+
+        // Reverse iterator
+
+        typedef reverse_iterator_base&lt;const_iterator&gt; const_reverse_iterator;
+
+        const_reverse_iterator rbegin () const;
+        const_reverse_iterator rend () const;
+
+        typedef reverse_iterator_base&lt;iterator&gt; reverse_iterator;
+
+        reverse_iterator rbegin ();
+        reverse_iterator rend ();
+    };</code></pre>
+
+<h2><a name="coordinate_vector"></a>Coordinate Vector</h2>
+
+<h4>Description</h4>
+
+<p>The templated class <code>coordinate_vector&lt;T, IB, IA,
+TA&gt; </code>is the base container adaptor for compressed
+vectors. For a <em>n</em>-dimensional sorted coordinate vector
+and <em>0 &lt;= i &lt; n </em>the non-zero elements <em>v</em><sub><em>i</em></sub>
+are mapped to consecutive elements of the index and value
+container, i.e. for elements <em>k</em> = <em>v</em><sub><em>i</em></sub><sub><sub><em>1</em></sub></sub><sub><em>
+</em></sub>and <em>k + 1 = v</em><sub><em>i</em></sub><sub><sub><em>2</em></sub></sub><sub><em>
+</em></sub>of these containers holds <em>i</em><sub><em>1</em></sub><em>
+&lt; i</em><sub><em>2</em></sub>.</p>
+
+<h4>Example</h4>
+
+<pre>int main () {
+    using namespace boost::numeric::ublas;
+    coordinate_vector&lt;double&gt; v (3, 3);
+    for (int i = 0; i &lt; v.size (); ++ i) 
+        v (i) = i;
+    std::cout &lt;&lt; v &lt;&lt; std::endl;
+}</pre>
+
+<h4>Definition</h4>
+
+<p>Defined in the header vector_sparse.hpp.</p>
+
+<h4>Template parameters</h4>
+
+<table border="1">
+    <tr>
+        <th>Parameter </th>
+        <th>Description </th>
+        <th>Default </th>
+    </tr>
+    <tr>
+        <td><code>T</code> </td>
+        <td>The type of object stored in the coordinate vector. </td>
+        <td>&nbsp;</td>
+    </tr>
+    <tr>
+        <td><code>IB</code></td>
+        <td>The index base of the coordinate vector. <a
+        href="#coordinate_vector_1">[1]</a></td>
+        <td><code>0</code></td>
+    </tr>
+    <tr>
+        <td><code>IA</code></td>
+        <td>The type of the adapted array for indices. <a
+        href="#coordinate_vector_2">[2]</a></td>
+        <td><code>unbounded_array&lt;std::size_t&gt;</code></td>
+    </tr>
+    <tr>
+        <td><code>TA</code></td>
+        <td>The type of the adapted array for values. <a
+        href="#coordinate_vector_2">[2]</a></td>
+        <td><code>unbounded_array&lt;T&gt;</code></td>
+    </tr>
+</table>
+
+<h4>Model of</h4>
+
+<p><a href="container.htm#vector">Vector</a>. </p>
+
+<h4>Type requirements</h4>
+
+<p>None, except for those imposed by the requirements of <a
+href="container.htm#vector">Vector</a>.</p>
+
+<h4>Public base classes</h4>
+
+<p><code>vector_expression&lt;coordinate_vector&lt;T, IB, IA,
+TA&gt; &gt;</code> </p>
+
+<h4>Members</h4>
+
+<table border="1">
+    <tr>
+        <th>Member </th>
+        <th>Description </th>
+    </tr>
+    <tr>
+        <td><code>coordinate_vector ()</code> </td>
+        <td>Allocates a<code> coordinate_vector </code>that holds
+        zero elements.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_vector (size_type size, size_type
+        non_zeros)</code></td>
+        <td>Allocates a<code> coordinate_vector</code> that holds
+        at most <code>size</code> elements.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_vector (const coordinate_vector
+        &amp;v)</code></td>
+        <td>The copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector (size_type non_zeros, const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended copy constructor.</td>
+    </tr>
+    <tr>
+        <td><code>void resize (size_type size, size_type
+        non_zeros)</code></td>
+        <td>Reallocates a <code>coordinate_vector </code>to hold
+        at most <code>size</code> elements. The content of the <code>coordinate_vector
+        </code>is preserved.</td>
+    </tr>
+    <tr>
+        <td><code>size_type size () const</code></td>
+        <td>Returns the size of the <code>coordinate_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reference operator () (size_type i) const</code></td>
+        <td>Returns the value of the <code>i</code>-th element. </td>
+    </tr>
+    <tr>
+        <td><code>reference operator () (size_type i)</code></td>
+        <td>Returns a reference of the <code>i</code>-th element.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reference operator [] (size_type i) const</code></td>
+        <td>Returns the value of the <code>i</code>-th element. </td>
+    </tr>
+    <tr>
+        <td><code>reference operator [] (size_type i)</code></td>
+        <td>Returns a reference of the <code>i</code>-th element.
+        </td>
+    </tr>
+    <tr>
+        <td><code>coordinate_vector &amp;operator = (const
+        coordinate_vector &amp;v)</code></td>
+        <td>The assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>coordinate_vector &amp;assign_temporary
+        (coordinate_vector &amp;v)</code></td>
+        <td>Assigns a temporary. May change the coordinate vector
+        <code>v</code>.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector &amp;operator = (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>The extended assignment operator.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector &amp;assign (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Assigns a vector expression to the coordinate vector.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector &amp;operator += (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Adds the vector
+        expression to the coordinate vector.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector &amp;plus_assign (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Adds a vector expression to the coordinate vector.
+        Left and right hand side of the assignment should be
+        independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector &amp;operator -= (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>A computed assignment operator. Subtracts the vector
+        expression from the coordinate vector.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AE&gt;<br>
+        coordinate_vector &amp;minus_assign (const
+        vector_expression&lt;AE&gt; &amp;ae)</code></td>
+        <td>Subtracts a vector expression from the coordinate
+        vector. Left and right hand side of the assignment should
+        be independent.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        coordinate_vector &amp;operator *= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Multiplies the
+        coordinate vector with a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>template&lt;class AT&gt;<br>
+        coordinate_vector &amp;operator /= (const AT &amp;at)</code></td>
+        <td>A computed assignment operator. Divides the
+        coordinate vector through a scalar.</td>
+    </tr>
+    <tr>
+        <td><code>void swap (coordinate_vector &amp;v)</code></td>
+        <td>Swaps the contents of the coordinate vectors. </td>
+    </tr>
+    <tr>
+        <td><code>void insert (size_type i, const_reference t)</code></td>
+        <td>Inserts the value <code>t</code> at the <code>i</code>-th
+        element.</td>
+    </tr>
+    <tr>
+        <td><code>void erase (size_type i)</code></td>
+        <td>Erases the value at the <code>i</code>-th element.</td>
+    </tr>
+    <tr>
+        <td><code>void clear ()</code></td>
+        <td>Clears the coordinate vector.</td>
+    </tr>
+    <tr>
+        <td><code>const_iterator begin () const</code></td>
+        <td>Returns a <code>const_iterator</code> pointing to the
+        beginning of the <code>coordinate_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_iterator end () const</code></td>
+        <td>Returns a <code>const_iterator</code> pointing to the
+        end of the <code>coordinate_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator begin () </code></td>
+        <td>Returns a <code>iterator</code> pointing to the
+        beginning of the <code>coordinate_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>iterator end () </code></td>
+        <td>Returns a <code>iterator</code> pointing to the end
+        of the <code>coordinate_vector</code>. </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator rbegin () const</code></td>
+        <td>Returns a <code>const_reverse_iterator</code>
+        pointing to the beginning of the reversed <code>coordinate_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>const_reverse_iterator rend () const</code></td>
+        <td>Returns a <code>const_reverse_iterator</code>
+        pointing to the end of the reversed <code>coordinate_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator rbegin () </code></td>
+        <td>Returns a <code>reverse_iterator</code> pointing to
+        the beginning of the reversed <code>coordinate_vector</code>.
+        </td>
+    </tr>
+    <tr>
+        <td><code>reverse_iterator rend () </code></td>
+        <td>Returns a <code>reverse_iterator</code> pointing to
+        the end of the reversed <code>coordinate_vector</code>. </td>
+    </tr>
+</table>
+
+<h4>Notes</h4>
+
+<p><a name="#coordinate_vector_1">[1]</a> Supported parameters
+for the index base are <code>0 </code>and <code>1</code> at
+least. </p>
+
+<p><a name="#coordinate_vector_2">[2]</a> Supported parameters
+for the adapted array are <code>unbounded_array&lt;&gt;</code>, <code>bounded_array&lt;&gt;</code>
+and <code>std::vector&lt;&gt;</code>. </p>
+
+<h4>Interface</h4>
+
+<pre><code>    // Array based sparse vector class 
+    template&lt;class T, std::size_t IB, class IA, class TA&gt;
+    class coordinate_vector: 
+        public vector_expression&lt;coordinate_vector&lt;T, IB, IA, TA&gt; &gt; {
+    public:      
+        typedef std::size_t size_type;
+        typedef std::ptrdiff_t difference_type;
+        typedef T value_type;
+        typedef const T &amp;const_reference;
+        typedef T &amp;reference;
+        typedef const T *const_pointer;
+        typedef T *pointer;
+        typedef IA index_array_type;
+        typedef TA value_array_type;
+        typedef const coordinate_vector&lt;T, IB, IA, TA&gt; const_self_type;
+        typedef coordinate_vector&lt;T, IB, IA, TA&gt; self_type;
+        typedef const vector_const_reference&lt;const_self_type&gt; const_closure_type;
+        typedef vector_reference&lt;self_type&gt; closure_type;
+        typedef typename IA::const_iterator const_iterator_type;
+        typedef typename IA::iterator iterator_type;
+        typedef sparse_tag storage_category;
+
+        // Construction and destruction
+        coordinate_vector (); 
+        coordinate_vector (size_type size, size_type non_zeros = 0);
+        coordinate_vector (const coordinate_vector &amp;v);
+        template&lt;class AE&gt;
+        coordinate_vector (const vector_expression&lt;AE&gt; &amp;ae, size_type non_zeros = 0);
+
+        // Accessors
+        size_type size () const;
+        size_type non_zeros () const;
+        static size_type index_base ();
+        const index_array_type &amp;index_data () const;
+        index_array_type &amp;index_data ();
+        const value_array_type &amp;value_data () const;
+        value_array_type &amp;value_data ();
+
+        // Resizing
+        void resize (size_type size, size_type non_zeros = 0);
+
+        // Element access
+        const_reference operator () (size_type i) const;
+        reference operator () (size_type i);
+
+        const_reference operator [] (size_type i) const;
+        reference operator [] (size_type i);
+
+        // Assignment
+        coordinate_vector &amp;operator = (const coordinate_vector &amp;v);
+        coordinate_vector &amp;assign_temporary (coordinate_vector &amp;v);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;operator = (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;reset (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;assign (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;operator += (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;plus_assign (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;operator -= (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AE&gt;
+        coordinate_vector &amp;minus_assign (const vector_expression&lt;AE&gt; &amp;ae);
+        template&lt;class AT&gt;
+        coordinate_vector &amp;operator *= (const AT &amp;at);
+        template&lt;class AT&gt;
+        coordinate_vector &amp;operator /= (const AT &amp;at);
+
+        // Swapping
+        void swap (coordinate_vector &amp;v);
+        friend void swap (coordinate_vector &amp;v1, coordinate_vector &amp;v2);
+
+        // Element insertion and erasure
+        void insert (size_type i, const_reference t);
+        void erase (size_type i);
+        void clear ();
+
+        class const_iterator;
+        class iterator;
+
+        // Element lookup
+        const_iterator find (size_type i) const;
+        iterator find (size_type i);
+        const_iterator find_first (size_type i) const;
+        iterator find_first (size_type i);
+        const_iterator find_last (size_type i) const;
+        iterator find_last (size_type i);
+
+        // Iterators simply are pointers.
+
+        class const_iterator:
+            public container_const_reference&lt;coordinate_vector&gt;,
+            public bidirectional_iterator_base&lt;const_iterator, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename coordinate_vector::difference_type difference_type;
+            typedef typename coordinate_vector::value_type value_type;
+            typedef typename coordinate_vector::const_reference reference;
+            typedef typename coordinate_vector::const_pointer pointer;
+
+            // Construction and destruction
+            const_iterator ();
+            const_iterator (const coordinate_vector &amp;v, const const_iterator_type &amp;it);
+            const_iterator (const iterator &amp;it);
+
+            // Arithmetic
+            const_iterator &amp;operator ++ ();
+            const_iterator &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            // Index
+            size_type index () const;
+
+            // Assignment
+            const_iterator &amp;operator = (const const_iterator &amp;it);
+
+            // Comparison
+            bool operator == (const const_iterator &amp;it) const;
+        };
+
+        const_iterator begin () const;
+        const_iterator end () const;
+
+        class iterator:
+            public container_reference&lt;coordinate_vector&gt;,
+            public bidirectional_iterator_base&lt;iterator, value_type&gt; {
+        public:
+            typedef sparse_bidirectional_iterator_tag iterator_category;
+            typedef typename coordinate_vector::difference_type difference_type;
+            typedef typename coordinate_vector::value_type value_type;
+            typedef typename coordinate_vector::reference reference;
+            typedef typename coordinate_vector::pointer pointer;
+
+            // Construction and destruction
+            iterator ();
+            iterator (coordinate_vector &amp;v, const iterator_type &amp;it);
+
+            // Arithmetic
+            iterator &amp;operator ++ ();
+            iterator &amp;operator -- ();
+
+            // Dereference
+            reference operator * () const;
+
+            // Index
+            size_type index () const;
+
+            // Assignment
+            iterator &amp;operator = (const iterator &amp;it);
+
+            // Comparison
+            bool operator == (const iterator &amp;it) const;
+        };
+
+        iterator begin ();
+        iterator end ();
+
+        // Reverse iterator
+
+        typedef reverse_iterator_base&lt;const_iterator&gt; const_reverse_iterator;
+
+        const_reverse_iterator rbegin () const;
+        const_reverse_iterator rend () const;
+
+        typedef reverse_iterator_base&lt;iterator&gt; reverse_iterator;
+
+        reverse_iterator rbegin ();
+        reverse_iterator rend ();
+    };</code></pre>
+
 <hr>
 
 <p>Copyright (©) 2000-2002 Joerg Walter, Mathias Koch <br>
diff --git a/include/boost/numeric/ublas/banded.hpp b/include/boost/numeric/ublas/banded.hpp
index 03244972..6cc53065 100644
--- a/include/boost/numeric/ublas/banded.hpp
+++ b/include/boost/numeric/ublas/banded.hpp
@@ -989,7 +989,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -1003,11 +1003,18 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename detail::ct_if<boost::is_const<M>::value,
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
+#endif
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
         typedef const banded_adaptor<M> const_self_type;
         typedef banded_adaptor<M> self_type;
-        typedef const matrix_const_reference<const_self_type> const_closure_type;
-        typedef matrix_reference<self_type> closure_type;
+        typedef const_self_type const_closure_type;
+        typedef self_type closure_type;
         typedef typename storage_restrict_traits<typename M::storage_category,
                                                  packed_proxy_tag>::storage_category storage_category;
         typedef typename M::orientation_category orientation_category;
@@ -1041,11 +1048,11 @@ namespace boost { namespace numeric { namespace ublas {
             return upper_;
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -1811,7 +1818,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         size_type lower_;
         size_type upper_;
         static matrix_type nil_;
diff --git a/include/boost/numeric/ublas/concepts.hpp b/include/boost/numeric/ublas/concepts.hpp
index 544a7ada..1a90d8de 100644
--- a/include/boost/numeric/ublas/concepts.hpp
+++ b/include/boost/numeric/ublas/concepts.hpp
@@ -872,6 +872,7 @@ namespace boost { namespace numeric { namespace ublas {
     ZeroElement (double) {
         return 0.;
     }
+#ifndef BOOST_MSVC
     template<>
     std::complex<float>
     ZeroElement (std::complex<float>) {
@@ -880,8 +881,9 @@ namespace boost { namespace numeric { namespace ublas {
     template<>
     std::complex<double>
     ZeroElement (std::complex<double>) {
-        return std::complex<float> (0.);
+        return std::complex<double> (0.);
     }
+#endif
     template<>
     vector<float> 
     ZeroElement (vector<float>) {
@@ -892,6 +894,7 @@ namespace boost { namespace numeric { namespace ublas {
     ZeroElement (vector<double>) {
         return zero_vector<double> ();
     }
+#ifndef BOOST_MSVC
     template<>
     vector<std::complex<float> > 
     ZeroElement (vector<std::complex<float> >) {
@@ -902,6 +905,7 @@ namespace boost { namespace numeric { namespace ublas {
     ZeroElement (vector<std::complex<double> >) {
         return zero_vector<std::complex<double> > ();
     }
+#endif
     template<>
     matrix<float> 
     ZeroElement (matrix<float>) {
@@ -912,6 +916,7 @@ namespace boost { namespace numeric { namespace ublas {
     ZeroElement (matrix<double>) {
         return zero_matrix<double> ();
     }
+#ifndef BOOST_MSVC
     template<>
     matrix<std::complex<float> > 
     ZeroElement (matrix<std::complex<float> >) {
@@ -922,6 +927,7 @@ namespace boost { namespace numeric { namespace ublas {
     ZeroElement (matrix<std::complex<double> >) {
         return zero_matrix<std::complex<double> > ();
     }
+#endif
 
     template<class T>
     T 
@@ -936,6 +942,7 @@ namespace boost { namespace numeric { namespace ublas {
     OneElement (double) {
         return 1.;
     }
+#ifndef BOOST_MSVC
     template<>
     std::complex<float>
     OneElement (std::complex<float>) {
@@ -946,6 +953,7 @@ namespace boost { namespace numeric { namespace ublas {
     OneElement (std::complex<double>) {
         return std::complex<double> (1.);
     }
+#endif
     template<>
     matrix<float> 
     OneElement (matrix<float>) {
@@ -956,6 +964,7 @@ namespace boost { namespace numeric { namespace ublas {
     OneElement (matrix<double>) {
         return identity_matrix<double> ();
     }
+#ifndef BOOST_MSVC
     template<>
     matrix<std::complex<float> > 
     OneElement (matrix<std::complex<float> >) {
@@ -966,6 +975,7 @@ namespace boost { namespace numeric { namespace ublas {
     OneElement (matrix<std::complex<double> >) {
         return identity_matrix<std::complex<double> > ();
     }
+#endif
 
     template<class E1, class E2>
     bool 
@@ -1641,6 +1651,7 @@ namespace boost { namespace numeric { namespace ublas {
         VectorSpaceConcept<double, matrix<double> >::constraints ();
         LinearOperatorConcept<double, vector<double>, matrix<double> >::constraints ();
 
+#ifndef BOOST_MSVC
         AdditiveAbelianGroupConcept<std::complex<float> >::constraints ();
         CommutativeRingWithIdentityConcept<std::complex<float> >::constraints ();
         FieldConcept<std::complex<float> >::constraints ();
@@ -1656,6 +1667,7 @@ namespace boost { namespace numeric { namespace ublas {
         RingWithIdentityConcept<matrix<std::complex<double> > >::constraints (0);
         VectorSpaceConcept<std::complex<double>, matrix<std::complex<double> > >::constraints ();
         LinearOperatorConcept<std::complex<double>, vector<std::complex<double> >, matrix<std::complex<double> > >::constraints ();
+#endif
 #endif
     }
 
diff --git a/include/boost/numeric/ublas/config.hpp b/include/boost/numeric/ublas/config.hpp
index fd0ddb4d..3912208c 100644
--- a/include/boost/numeric/ublas/config.hpp
+++ b/include/boost/numeric/ublas/config.hpp
@@ -24,122 +24,21 @@
 #include <boost/type_traits.hpp>
 #include <boost/limits.hpp>
 
-#ifdef BOOST_MSVC
-// Disable some MSVC specific warnings.
-#pragma warning (disable: 4355)
-#pragma warning (disable: 4503)
-#pragma warning (disable: 4786)
-#endif
-
-
-
-#ifdef BOOST_MSVC
-// MSVC doesn't always accept the 'typename' keyword.
-#define BOOST_UBLAS_TYPENAME
-// MSVC doesn't accept the 'using' keyword (at least for importing base members).
-#define BOOST_UBLAS_USING
-#else
-#define BOOST_UBLAS_TYPENAME typename
-#define BOOST_UBLAS_USING using
-#endif
-// This could be eliminated.
-#define BOOST_UBLAS_EXPLICIT explicit
-
-
-
-// #ifdef BOOST_MSVC
-// With MSVC we can could perform IO via basic_stream
-// #define BOOST_UBLAS_USE_BASIC_STREAM
-// #else
-// IO via streams
-#define BOOST_UBLAS_USE_STREAM
-// #endif
-
-
-
-// Enable different sparse element proxies
-// These fix a [1] = a [0] = 1, but probably won't work on broken compilers.
-// Thanks to Marc Duflot for spotting this.
-// #define BOOST_UBLAS_STRICT_STORAGE_SPARSE
-#define BOOST_UBLAS_STRICT_VECTOR_SPARSE
-#define BOOST_UBLAS_STRICT_MATRIX_SPARSE
-
-
-
-// Enable performance options in release mode
-#ifdef NDEBUG
-
-#ifdef BOOST_MSVC
-// MSVC has special inlining options
-#pragma inline_recursion (on)
-#pragma inline_depth (255)
-#pragma auto_inline (on)
-// #define BOOST_UBLAS_INLINE __forceinline
-#define BOOST_UBLAS_INLINE __inline
-#else
-#define BOOST_UBLAS_INLINE inline
-#endif
-
-// Do not check sizes!
-#define BOOST_UBLAS_USE_FAST_SAME
-
-// Use expression templates.
-#define BOOST_UBLAS_USE_ET
-
-// Disable performance options in debug mode
-#else
-
-#ifdef BOOST_MSVC
-// MSVC has special inlining options
-// #pragma inline_recursion (off)
-// #pragma inline_depth ()
-// #pragma auto_inline (off)
-#endif
-#define BOOST_UBLAS_INLINE
-
-#ifdef BOOST_MSVC
-// Use expression templates (otherwise we get many ICE's)
-#define BOOST_UBLAS_USE_ET
-#endif
-
-// Bounds check
-#define BOOST_UBLAS_BOUNDS_CHECK
-// Type check for non dense matrices
-#define BOOST_UBLAS_TYPE_CHECK
-
-#endif
-
-
-
-// Use invariant hoisting.
-#define BOOST_UBLAS_USE_INVARIANT_HOISTING
-
-// Use Duff's device
-// #define BOOST_UBLAS_USE_DUFF_DEVICE
-
-// Choose evaluation method for dense vectors and matrices
-#define BOOST_UBLAS_USE_INDEXING
-// #define BOOST_UBLAS_USE_ITERATING
-// #define BOOST_UBLAS_ITERATOR_THRESHOLD 0
-// #define BOOST_UBLAS_ITERATOR_THRESHOLD (std::numeric_limits<std::ptrdiff_t>::max ())
-
-// ET options
-#define BOOST_UBLAS_ET_VALUE
-// #define BOOST_UBLAS_ET_REFERENCE
-// #define BOOST_UBLAS_ET_CLOSURE_VALUE
-#define BOOST_UBLAS_ET_CLOSURE_REFERENCE
-
-// Use indexed iterators.
-// #define BOOST_UBLAS_USE_INDEXED_ITERATOR
-
 
 
 // Compiler specific problems
 #ifdef BOOST_MSVC
 
+// Disable some MSVC specific warnings.
+#pragma warning (disable: 4355)
+#pragma warning (disable: 4503)
+#pragma warning (disable: 4786)
+
 // Open problems:
 // MSVC lacks some specializations in <cmath>
 #define BOOST_UBLAS_C_MATH
+// MSVC lacks std::abort() in <cstdlib>
+#define BOOST_UBLAS_C_STDLIB
 // MSVC allows to implement free function as friends.
 #define BOOST_UBLAS_FRIEND_FUNCTION
 
@@ -247,6 +146,112 @@
 
 
 
+#ifdef BOOST_MSVC
+// MSVC doesn't always accept the 'typename' keyword.
+#define BOOST_UBLAS_TYPENAME
+// MSVC doesn't accept the 'using' keyword (at least for importing base members).
+#define BOOST_UBLAS_USING
+#else
+#define BOOST_UBLAS_TYPENAME typename
+#define BOOST_UBLAS_USING using
+#endif
+// This could be eliminated.
+#define BOOST_UBLAS_EXPLICIT explicit
+
+
+
+// #ifdef BOOST_MSVC
+// With MSVC we could perform IO via basic_stream
+// #define BOOST_UBLAS_USE_BASIC_STREAM
+// #else
+// IO via streams
+#define BOOST_UBLAS_USE_STREAM
+// #endif
+
+// Enable assignment of non conformant proxies
+// Thanks to Michael Stevens for spotting this.
+#define BOOST_UBLAS_NON_CONFORMANT_PROXIES
+
+// Enable different sparse element proxies
+// These fix a [1] = a [0] = 1, but probably won't work on broken compilers.
+// Thanks to Marc Duflot for spotting this.
+#ifndef BOOST_MSVC_STD_ITERATOR
+// #define BOOST_UBLAS_STRICT_STORAGE_SPARSE
+#define BOOST_UBLAS_STRICT_VECTOR_SPARSE
+#define BOOST_UBLAS_STRICT_MATRIX_SPARSE
+#endif
+
+// Enable compile time typedefs for proxies
+#ifndef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+// Doesn't work under Borland
+#ifndef __BORLANDC__
+#define BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
+#define BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+#endif
+#endif
+
+// Enable performance options in release mode
+#ifdef NDEBUG
+
+#ifdef BOOST_MSVC
+// MSVC has special inlining options
+#pragma inline_recursion (on)
+#pragma inline_depth (255)
+#pragma auto_inline (on)
+// #define BOOST_UBLAS_INLINE __forceinline
+#define BOOST_UBLAS_INLINE __inline
+#else
+#define BOOST_UBLAS_INLINE inline
+#endif
+
+// Do not check sizes!
+#define BOOST_UBLAS_USE_FAST_SAME
+
+// Use expression templates.
+#define BOOST_UBLAS_USE_ET
+
+// Disable performance options in debug mode
+#else
+
+#ifdef BOOST_MSVC
+// MSVC has special inlining options
+// #pragma inline_recursion (off)
+// #pragma inline_depth ()
+// #pragma auto_inline (off)
+#endif
+#define BOOST_UBLAS_INLINE
+
+#ifdef BOOST_MSVC
+// Use expression templates (otherwise we get many ICE's)
+#define BOOST_UBLAS_USE_ET
+#endif
+
+// Bounds check
+#define BOOST_UBLAS_BOUNDS_CHECK
+// Type check for non dense matrices
+#define BOOST_UBLAS_TYPE_CHECK
+
+#endif
+
+
+
+// Use invariant hoisting.
+#define BOOST_UBLAS_USE_INVARIANT_HOISTING
+
+// Use Duff's device
+// #define BOOST_UBLAS_USE_DUFF_DEVICE
+
+// Choose evaluation method for dense vectors and matrices
+#define BOOST_UBLAS_USE_INDEXING
+// #define BOOST_UBLAS_USE_ITERATING
+// #define BOOST_UBLAS_ITERATOR_THRESHOLD 0
+// #define BOOST_UBLAS_ITERATOR_THRESHOLD (std::numeric_limits<std::ptrdiff_t>::max ())
+
+// Use indexed iterators.
+// #define BOOST_UBLAS_USE_INDEXED_ITERATOR
+
+
+
 // Forward declarations
 namespace boost { namespace numeric { namespace ublas {
 
@@ -325,10 +330,10 @@ namespace boost { namespace numeric { namespace ublas {
     template<class T, class A = map_array<std::size_t, T> >
     class sparse_vector;
 
-    template<class T, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T>, std::size_t IB = 0>
+    template<class T, std::size_t IB = 0, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T> >
     class compressed_vector;
 
-    template<class T, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T>, std::size_t IB = 0>
+    template<class T, std::size_t IB = 0, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T> >
     class coordinate_vector;
 
     struct unknown_orientation_tag {};
@@ -390,10 +395,10 @@ namespace boost { namespace numeric { namespace ublas {
     template<class T, class F = row_major, class A = map_array<std::size_t, map_array<std::size_t, T> > >
     class sparse_vector_of_sparse_vector;
 
-    template<class T, class F = row_major, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T>, std::size_t IB = 0>
+    template<class T, class F = row_major, std::size_t IB = 0, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T> >
     class compressed_matrix;
 
-    template<class T, class F = row_major, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T>, std::size_t IB = 0>
+    template<class T, class F = row_major, std::size_t IB = 0, class IA = unbounded_array<std::size_t>, class TA = unbounded_array<T> >
     class coordinate_matrix;
 
     template<class V>
diff --git a/include/boost/numeric/ublas/exception.hpp b/include/boost/numeric/ublas/exception.hpp
index 0a572177..ff9e21bd 100644
--- a/include/boost/numeric/ublas/exception.hpp
+++ b/include/boost/numeric/ublas/exception.hpp
@@ -17,6 +17,7 @@
 #ifndef BOOST_UBLAS_EXCEPTION_H
 #define BOOST_UBLAS_EXCEPTION_H
 
+#include <cstdlib>
 #include <exception>
 #include <stdexcept>
 #include <iostream>
@@ -39,7 +40,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -57,7 +62,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -75,7 +84,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -93,7 +106,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -111,7 +128,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -129,7 +150,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -147,7 +172,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
@@ -165,7 +194,11 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_NO_EXCEPTIONS) && defined (BOOST_UBLAS_USE_EXCEPTIONS)
             throw *this;
 #else
-            abort ();
+#ifdef BOOST_UBLAS_C_STDLIB
+            ::abort ();
+#else
+            std::abort ();
+#endif
 #endif
         }
     };
diff --git a/include/boost/numeric/ublas/hermitian.hpp b/include/boost/numeric/ublas/hermitian.hpp
index 45fe6e2b..c2d22b4a 100644
--- a/include/boost/numeric/ublas/hermitian.hpp
+++ b/include/boost/numeric/ublas/hermitian.hpp
@@ -933,7 +933,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         // FIXME: no better way to not return the address of a temporary?
         // typedef typename M::const_reference const_reference;
         typedef typename M::value_type const_reference;
@@ -954,11 +954,18 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename detail::ct_if<boost::is_const<M>::value,
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
+#endif
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
         typedef const hermitian_adaptor<M, F> const_self_type;
         typedef hermitian_adaptor<M, F> self_type;
-        typedef const matrix_const_reference<const_self_type> const_closure_type;
-        typedef matrix_reference<self_type> closure_type;
+        typedef const_self_type const_closure_type;
+        typedef self_type closure_type;
         typedef typename storage_restrict_traits<typename M::storage_category,
                                                  packed_proxy_tag>::storage_category storage_category;
         typedef typename F::packed_category packed_category;
@@ -966,7 +973,7 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Construction and destruction
         BOOST_UBLAS_INLINE
-        hermitian_adaptor (): 
+        hermitian_adaptor ():
             data_ (nil_) {
             BOOST_UBLAS_CHECK (data_.size1 () == data_.size2 (), bad_size ());
         }
@@ -991,11 +998,11 @@ namespace boost { namespace numeric { namespace ublas {
             return data_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -1717,7 +1724,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         static matrix_type nil_;
         static value_type conj_;
     };
diff --git a/include/boost/numeric/ublas/matrix_assign.hpp b/include/boost/numeric/ublas/matrix_assign.hpp
index b537a330..ec28e152 100644
--- a/include/boost/numeric/ublas/matrix_assign.hpp
+++ b/include/boost/numeric/ublas/matrix_assign.hpp
@@ -695,6 +695,7 @@ namespace boost { namespace numeric { namespace ublas {
                         functor_type () (*it2, value_type ());
                         ++ it2;
                     } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                         // Sparse proxies don't need to be conformant.
                         // Thanks to Michael Stevens for suggesting this.
                         size_type index1 (it2e.index1 ()), index2 (it2e.index2 ());
@@ -702,8 +703,10 @@ namespace boost { namespace numeric { namespace ublas {
                         restart (m, index1, index2, it1, it1_end, it2, it2_end, row_major_tag ());
                         // The proxies could reference the same container.
                         restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
+#endif
                     }
                 }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 while (it2e != it2e_end) {
                     // Sparse proxies don't need to be conformant.
                     // Thanks to Michael Stevens for suggesting this.
@@ -712,6 +715,7 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
                 }
+#endif
                 while (it2 != it2_end) {
                     functor_type () (*it2, value_type ());
                     ++ it2;
@@ -726,6 +730,7 @@ namespace boost { namespace numeric { namespace ublas {
                 }
                 ++ it1;
             } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 typename E::const_iterator2 it2e (it1e.begin ());
                 typename E::const_iterator2 it2e_end (it1e.end ());
                 while (it2e != it2e_end) {
@@ -736,9 +741,11 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
                 }
+#endif
                 ++ it1e;
             }
         }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
         while (it1e != it1e_end) {
             typename E::const_iterator2 it2e (it1e.begin ());
             typename E::const_iterator2 it2e_end (it1e.end ());
@@ -752,6 +759,7 @@ namespace boost { namespace numeric { namespace ublas {
             }
             ++ it1e;
         }
+#endif
         while (it1 != it1_end) {
             typename M::iterator2 it2 (it1.begin ());
             typename M::iterator2 it2_end (it1.end ());
@@ -800,6 +808,7 @@ namespace boost { namespace numeric { namespace ublas {
                         functor_type () (*it1, value_type ());
                         ++ it1;
                     } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                         // Sparse proxies don't need to be conformant.
                         // Thanks to Michael Stevens for suggesting this.
                         size_type index1 (it1e.index1 ()), index2 (it1e.index2 ());
@@ -807,8 +816,10 @@ namespace boost { namespace numeric { namespace ublas {
                         restart (m, index1, index2, it2, it2_end, it1, it1_end, column_major_tag ());
                         // The proxies could reference the same container.
                         restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
+#endif
                     }
                 }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 while (it1e != it1e_end) {
                     // Sparse proxies don't need to be conformant.
                     // Thanks to Michael Stevens for suggesting this.
@@ -817,6 +828,7 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
                 }
+#endif
                 while (it1 != it1_end) {
                     functor_type () (*it1, value_type ());
                     ++ it1;
@@ -831,6 +843,7 @@ namespace boost { namespace numeric { namespace ublas {
                 }
                 ++ it2;
             } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 typename E::const_iterator1 it1e (it2e.begin ());
                 typename E::const_iterator1 it1e_end (it2e.end ());
                 while (it1e != it1e_end) {
@@ -841,9 +854,11 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
                 }
+#endif
                 ++ it2e;
             }
         }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
         while (it2e != it2e_end) {
             typename E::const_iterator1 it1e (it2e.begin ());
             typename E::const_iterator1 it1e_end (it2e.end ());
@@ -857,6 +872,7 @@ namespace boost { namespace numeric { namespace ublas {
             }
             ++ it2e;
         }
+#endif
         while (it2 != it2_end) {
             typename M::iterator1 it1 (it2.begin ());
             typename M::iterator1 it1_end (it2.end ());
@@ -1008,6 +1024,7 @@ namespace boost { namespace numeric { namespace ublas {
                         functor_type () (*it2, *it2e);
                         ++ it2, ++ it2e;
                     } else if (compare < 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                         // Sparse proxies don't need to be conformant.
                         // Thanks to Michael Stevens for suggesting this.
                         size_type index1 (it2.index1 ()), index2 (it2.index2 ());
@@ -1015,7 +1032,9 @@ namespace boost { namespace numeric { namespace ublas {
                         restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
                         // The proxies could reference the same container.
                         restart (m, index1, index2, it1, it1_end, it2, it2_end, row_major_tag ());
+#endif
                     } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                         // Sparse proxies don't need to be conformant.
                         // Thanks to Michael Stevens for suggesting this.
                         size_type index1 (it2e.index1 ()), index2 (it2e.index2 ());
@@ -1023,8 +1042,10 @@ namespace boost { namespace numeric { namespace ublas {
                         restart (m, index1, index2, it1, it1_end, it2, it2_end, row_major_tag ());
                         // The proxies could reference the same container.
                         restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
+#endif
                     }
                 }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 while (it2e != it2e_end) {
                     // Sparse proxies don't need to be conformant.
                     // Thanks to Michael Stevens for suggesting this.
@@ -1043,8 +1064,10 @@ namespace boost { namespace numeric { namespace ublas {
                     restart (m, index1, index2, it1, it1_end, it2, it2_end, row_major_tag ());
                     restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
                 }
+#endif
                 ++ it1, ++ it1e;
             } else if (compare < 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 typename M::iterator2 it2 (it1.begin ());
                 typename M::iterator2 it2_end (it1.end ());
                 while (it2 != it2_end) {
@@ -1055,8 +1078,10 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (m, index1, index2, it1, it1_end, it2, it2_end, row_major_tag ());
                 }
+#endif
                 ++ it1;
             } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 typename E::iterator2 it2e (it1e.begin ());
                 typename E::iterator2 it2e_end (it1e.end ());
                 while (it2e != it2e_end) {
@@ -1067,9 +1092,11 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (e, index1, index2, it1e, it1e_end, it2e, it2e_end, row_major_tag ());
                 }
+#endif
                 ++ it1e;
             }
         }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
         while (it1e != it1e_end) {
             typename E::iterator2 it2e (it1e.begin ());
             typename E::iterator2 it2e_end (it1e.end ());
@@ -1096,6 +1123,7 @@ namespace boost { namespace numeric { namespace ublas {
             }
             ++ it1;
         }
+#endif
     }
     // Sparse (proxy) column major case
     template<class F, class M, class E>
@@ -1124,6 +1152,7 @@ namespace boost { namespace numeric { namespace ublas {
                         functor_type () (*it1, *it1e);
                         ++ it1, ++ it1e;
                     } else if (compare < 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                         // Sparse proxies don't need to be conformant.
                         // Thanks to Michael Stevens for suggesting this.
                         size_type index1 (it1.index1 ()), index2 (it1.index2 ());
@@ -1131,7 +1160,9 @@ namespace boost { namespace numeric { namespace ublas {
                         restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
                         // The proxies could reference the same container.
                         restart (m, index1, index2, it2, it2_end, it1, it1_end, column_major_tag ());
+#endif
                     } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                         // Sparse proxies don't need to be conformant.
                         // Thanks to Michael Stevens for suggesting this.
                         size_type index1 (it1e.index1 ()), index2 (it1e.index2 ());
@@ -1139,8 +1170,10 @@ namespace boost { namespace numeric { namespace ublas {
                         restart (m, index1, index2, it2, it2_end, it1, it1_end, column_major_tag ());
                         // The proxies could reference the same container.
                         restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
+#endif
                     }
                 }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 while (it1e != it1e_end) {
                     // Sparse proxies don't need to be conformant.
                     // Thanks to Michael Stevens for suggesting this.
@@ -1159,8 +1192,10 @@ namespace boost { namespace numeric { namespace ublas {
                     restart (m, index1, index2, it2, it2_end, it1, it1_end, column_major_tag ());
                     restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
                 }
+#endif
                 ++ it2, ++ it2e;
             } else if (compare < 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 typename M::iterator1 it1 (it2.begin ());
                 typename M::iterator1 it1_end (it2.end ());
                 while (it1 != it1_end) {
@@ -1171,8 +1206,10 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (m, index1, index2, it2, it2_end, it1, it1_end, column_major_tag ());
                 }
+#endif
                 ++ it2;
             } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 typename E::iterator1 it1e (it2e.begin ());
                 typename E::iterator1 it1e_end (it2e.end ());
                 while (it1e != it1e_end) {
@@ -1183,9 +1220,11 @@ namespace boost { namespace numeric { namespace ublas {
                     // The proxies could reference the same container.
                     restart (e, index1, index2, it2e, it2e_end, it1e, it1e_end, column_major_tag ());
                 }
+#endif
                 ++ it2e;
             }
         }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
         while (it2e != it2e_end) {
             typename E::iterator1 it1e (it2e.begin ());
             typename E::iterator1 it1e_end (it2e.end ());
@@ -1212,6 +1251,7 @@ namespace boost { namespace numeric { namespace ublas {
             }
             ++ it2;
         }
+#endif
     }
 
     // Dispatcher
diff --git a/include/boost/numeric/ublas/matrix_expression.hpp b/include/boost/numeric/ublas/matrix_expression.hpp
index 6da89536..e96ea81b 100644
--- a/include/boost/numeric/ublas/matrix_expression.hpp
+++ b/include/boost/numeric/ublas/matrix_expression.hpp
@@ -243,7 +243,7 @@ namespace boost { namespace numeric { namespace ublas {
         // Element access
         BOOST_UBLAS_INLINE
         const_reference operator () (size_type i, size_type j) const {
-            return e_ (i, j);
+            return expression () (i, j);
         }
 
         typedef const_iterator1_type const_iterator1;
@@ -254,39 +254,39 @@ namespace boost { namespace numeric { namespace ublas {
         // Element lookup
         BOOST_UBLAS_INLINE
         const_iterator1 find_first1 (int rank, size_type i, size_type j) const {
-            return const_iterator1 (e_.find_first1 (rank, i, j));
+            return expression ().find_first1 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         const_iterator1 find_last1 (int rank, size_type i, size_type j) const {
-            return const_iterator1 (e_.find_last1 (rank, i, j));
+            return expression ().find_last1 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         const_iterator2 find_first2 (int rank, size_type i, size_type j) const {
-            return const_iterator2 (e_.find_first2 (rank, i, j));
+            return expression ().find_first2 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         const_iterator2 find_last2 (int rank, size_type i, size_type j) const {
-            return const_iterator2 (e_.find_last2 (rank, i, j));
+            return expression ().find_last2 (rank, i, j);
         }
 
         // Iterators are the iterators of the referenced expression.
 
         BOOST_UBLAS_INLINE
         const_iterator1 begin1 () const {
-            return e_.begin1 ();
+            return expression ().begin1 ();
         }
         BOOST_UBLAS_INLINE
         const_iterator1 end1 () const {
-            return e_.end1 ();
+            return expression ().end1 ();
         }
 
         BOOST_UBLAS_INLINE
         const_iterator2 begin2 () const {
-            return e_.begin2 ();
+            return expression ().begin2 ();
         }
         BOOST_UBLAS_INLINE
         const_iterator2 end2 () const {
-            return e_.end2 ();
+            return expression ().end2 ();
         }
 
         // Reverse iterators
@@ -380,17 +380,17 @@ namespace boost { namespace numeric { namespace ublas {
         // Resizing
         BOOST_UBLAS_INLINE
         void resize (size_type size1, size_type size2) {
-            e_.resize (size1, size2);
+            expression ().resize (size1, size2);
         }
 
         // Element access
         BOOST_UBLAS_INLINE
         const_reference operator () (size_type i, size_type j) const {
-            return e_ (i, j);
+            return expression () (i, j);
         }
         BOOST_UBLAS_INLINE
         reference operator () (size_type i, size_type j) {
-            return e_ (i, j);
+            return expression () (i, j);
         }
 
         typedef const_iterator1_type const_iterator1;
@@ -401,73 +401,73 @@ namespace boost { namespace numeric { namespace ublas {
         // Element lookup
         BOOST_UBLAS_INLINE
         const_iterator1 find_first1 (int rank, size_type i, size_type j) const {
-            return const_iterator1 (e_.find_first1 (rank, i, j));
+            return expression ().find_first1 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         iterator1 find_first1 (int rank, size_type i, size_type j) {
-            return iterator1 (e_.find_first1 (rank, i, j));
+            return expression ().find_first1 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         const_iterator1 find_last1 (int rank, size_type i, size_type j) const {
-            return const_iterator1 (e_.find_last1 (rank, i, j));
+            return expression ().find_last1 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         iterator1 find_last1 (int rank, size_type i, size_type j) {
-            return iterator1 (e_.find_last1 (rank, i, j));
+            return expression ().find_last1 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         const_iterator2 find_first2 (int rank, size_type i, size_type j) const {
-            return const_iterator2 (e_.find_first2 (rank, i, j));
+            return expression ().find_first2 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         iterator2 find_first2 (int rank, size_type i, size_type j) {
-            return iterator2 (e_.find_first2 (rank, i, j));
+            return expression ().find_first2 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         const_iterator2 find_last2 (int rank, size_type i, size_type j) const {
-            return const_iterator2 (e_.find_last2 (rank, i, j));
+            return expression ().find_last2 (rank, i, j);
         }
         BOOST_UBLAS_INLINE
         iterator2 find_last2 (int rank, size_type i, size_type j) {
-            return iterator2 (e_.find_last2 (rank, i, j));
+            return expression ().find_last2 (rank, i, j);
         }
 
         // Iterators are the iterators of the referenced expression.
 
         BOOST_UBLAS_INLINE
         const_iterator1 begin1 () const {
-            return e_.begin1 ();
+            return expression ().begin1 ();
         }
         BOOST_UBLAS_INLINE
         const_iterator1 end1 () const {
-            return e_.end1 ();
+            return expression ().end1 ();
         }
 
         BOOST_UBLAS_INLINE
         iterator1 begin1 () {
-            return e_.begin1 ();
+            return expression ().begin1 ();
         }
         BOOST_UBLAS_INLINE
         iterator1 end1 () {
-            return e_.end1 ();
+            return expression ().end1 ();
         }
 
         BOOST_UBLAS_INLINE
         const_iterator2 begin2 () const {
-            return e_.begin2 ();
+            return expression ().begin2 ();
         }
         BOOST_UBLAS_INLINE
         const_iterator2 end2 () const {
-            return e_.end2 ();
+            return expression ().end2 ();
         }
 
         BOOST_UBLAS_INLINE
         iterator2 begin2 () {
-            return e_.begin2 ();
+            return expression ().begin2 ();
         }
         BOOST_UBLAS_INLINE
         iterator2 end2 () {
-            return e_.end2 ();
+            return expression ().end2 ();
         }
 
         // Reverse iterators
@@ -557,7 +557,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const vector_matrix_binary<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef unknown_orientation_tag orientation_category;
         typedef typename E1::const_iterator const_iterator1_type;
         typedef typename E2::const_iterator const_iterator2_type;
@@ -581,11 +584,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e2_.size ();
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -1026,20 +1029,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class E1, class E2, class F>
     struct vector_matrix_binary_traits {
-        typedef vector_matrix_binary<typename E1::const_closure_type, 
-                                     typename E2::const_closure_type, F> expression_type;
+        typedef vector_matrix_binary<E1, E2, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -1073,7 +1069,9 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E::const_closure_type expression_closure_type;
         typedef const matrix_unary1<E, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E::orientation_category orientation_category;
         typedef typename E::const_iterator1 const_iterator1_type;
         typedef typename E::const_iterator2 const_iterator2_type;
@@ -1097,7 +1095,7 @@ namespace boost { namespace numeric { namespace ublas {
             return e_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const expression_type &expression () const {
+        const expression_closure_type &expression () const {
             return e_;
         }
 
@@ -1436,17 +1434,12 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression_type e_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression_type &e_;
-#endif
+        expression_closure_type e_;
     };
 
     template<class E, class F>
     struct matrix_unary1_traits {
-        typedef matrix_unary1<typename E::const_closure_type, F> expression_type;
+        typedef matrix_unary1<E, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -1506,7 +1499,9 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E::const_closure_type expression_closure_type;
         typedef const matrix_unary2<E, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E::orientation_category orientation_category;
         typedef typename E::const_iterator1 const_iterator2_type;
         typedef typename E::const_iterator2 const_iterator1_type;
@@ -1530,7 +1525,7 @@ namespace boost { namespace numeric { namespace ublas {
             return e_.size1 ();
         }
         BOOST_UBLAS_INLINE
-        const expression_type &expression () const {
+        const expression_closure_type &expression () const {
             return e_;
         }
 
@@ -1869,17 +1864,12 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression_type e_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression_type &e_;
-#endif
+        expression_closure_type e_;
     };
 
     template<class E, class F>
     struct matrix_unary2_traits {
-        typedef matrix_unary2<typename E::const_closure_type, F> expression_type;
+        typedef matrix_unary2<E, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -1922,7 +1912,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const matrix_binary<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef unknown_orientation_tag orientation_category;
         typedef typename E1::const_iterator1 const_iterator11_type;
         typedef typename E1::const_iterator2 const_iterator12_type;
@@ -1948,11 +1941,11 @@ namespace boost { namespace numeric { namespace ublas {
             return BOOST_UBLAS_SAME (e1_.size2 (), e2_.size2 ());
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -2571,20 +2564,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class E1, class E2, class F>
     struct matrix_binary_traits {
-        typedef matrix_binary<typename E1::const_closure_type, 
-                              typename E2::const_closure_type, F> expression_type;
+        typedef matrix_binary<E1, E2, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -2633,7 +2619,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef E1 expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const matrix_binary_scalar1<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E2::orientation_category orientation_category;
         typedef typename E1::value_type const_iterator1_type;
         typedef typename E2::const_iterator1 const_iterator21_type;
@@ -2658,11 +2647,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e2_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -3013,20 +3002,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class T1, class E2, class F>
     struct matrix_binary_scalar1_traits {
-        typedef matrix_binary_scalar1<scalar_const_reference<T1>,
-                                      typename E2::const_closure_type, F> expression_type;
+        typedef matrix_binary_scalar1<scalar_const_reference<T1>, E2, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type;
 #else
@@ -3061,7 +3043,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef E2 expression2_closure_type;
         typedef const matrix_binary_scalar2<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E1::orientation_category orientation_category;
         typedef typename E1::const_iterator1 const_iterator11_type;
         typedef typename E1::const_iterator2 const_iterator12_type;
@@ -3086,11 +3071,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e1_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -3441,20 +3426,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class E1, class T2, class F>
     struct matrix_binary_scalar2_traits {
-        typedef matrix_binary_scalar2<typename E1::const_closure_type,
-                                      scalar_const_reference<T2>, F> expression_type;
+        typedef matrix_binary_scalar2<E1, scalar_const_reference<T2>, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -3499,7 +3477,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const matrix_vector_binary1<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E1::const_iterator1 const_iterator1_type;
         typedef typename E2::const_iterator const_iterator2_type;
         typedef unknown_storage_tag storage_category;
@@ -3518,11 +3499,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e1_.size1 (); 
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -3755,14 +3736,8 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class T1, class E1, class T2, class E2>
@@ -3770,9 +3745,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<typename E1::const_closure_type,
-                                      typename E2::const_closure_type,
-                                      matrix_vector_prod1<T1, T2, promote_type> > expression_type;
+        typedef matrix_vector_binary1<E1, E2, matrix_vector_prod1<T1, T2, promote_type> > expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type;
 #else
@@ -3867,7 +3840,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const matrix_vector_binary2<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E1::const_iterator const_iterator1_type;
         typedef typename E2::const_iterator2 const_iterator2_type;
         typedef unknown_storage_tag storage_category;
@@ -3886,11 +3862,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e2_.size2 (); 
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -4124,14 +4100,8 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class T1, class E1, class T2, class E2>
@@ -4139,9 +4109,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<typename E1::const_closure_type,
-                                      typename E2::const_closure_type,
-                                      matrix_vector_prod2<T1, T2, promote_type> > expression_type;
+        typedef matrix_vector_binary2<E1, E2, matrix_vector_prod2<T1, T2, promote_type> > expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type;
 #else
@@ -4236,7 +4204,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const matrix_matrix_binary<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef unknown_orientation_tag orientation_category;
         typedef typename E1::const_iterator1 const_iterator11_type;
         typedef typename E1::const_iterator2 const_iterator12_type;
@@ -4262,11 +4233,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e2_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -4767,14 +4738,8 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class T1, class E1, class T2, class E2>
@@ -4782,9 +4747,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<typename E1::const_closure_type,
-                                     typename E2::const_closure_type,
-                                     matrix_matrix_prod<T1, T2, promote_type> > expression_type;
+        typedef matrix_matrix_binary<E1, E2, matrix_matrix_prod<T1, T2, promote_type> > expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type;
 #else
@@ -4869,6 +4832,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef E expression_type;
         typedef F functor_type;
         typedef typename F::result_type value_type;
+        typedef typename E::const_closure_type expression_closure_type;
 
         // Construction and destruction
         BOOST_UBLAS_INLINE
@@ -4880,7 +4844,7 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Accessors
         BOOST_UBLAS_INLINE
-        const expression_type &expression () const {
+        const expression_closure_type &expression () const {
             return e_;
         }
 
@@ -4890,17 +4854,12 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression_type e_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression_type &e_;
-#endif
+        expression_closure_type e_;
     };
     
     template<class E, class F>
     struct matrix_scalar_unary_traits {
-        typedef matrix_scalar_unary<typename E::const_closure_type, F> expression_type;
+        typedef matrix_scalar_unary<E, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
          typedef expression_type result_type; 
 #else
diff --git a/include/boost/numeric/ublas/matrix_proxy.hpp b/include/boost/numeric/ublas/matrix_proxy.hpp
index 5b9609a9..65a69dea 100644
--- a/include/boost/numeric/ublas/matrix_proxy.hpp
+++ b/include/boost/numeric/ublas/matrix_proxy.hpp
@@ -38,7 +38,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -53,15 +53,16 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const matrix_row<matrix_type> > const_closure_type;
-        typedef vector_reference<matrix_row<matrix_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_row<matrix_type> const_closure_type;
         typedef matrix_row<matrix_type> closure_type;
-#endif
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_iterator2 const_iterator_type;
         typedef typename M::iterator2 iterator_type;
 #else
@@ -94,11 +95,11 @@ namespace boost { namespace numeric { namespace ublas {
             return i_;
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -509,7 +510,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         size_type i_;
         static matrix_type nil_;
     };
@@ -526,7 +527,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M>
     BOOST_UBLAS_INLINE
-    const matrix_row<const M> row (const M &data, std::size_t i) {
+    matrix_row<const M> row (const M &data, std::size_t i) {
         return matrix_row<const M> (data, i);
     }
 #endif
@@ -544,7 +545,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -559,15 +560,16 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const matrix_column<matrix_type> > const_closure_type;
-        typedef vector_reference<matrix_column<matrix_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_column<matrix_type> const_closure_type;
         typedef matrix_column<matrix_type> closure_type;
-#endif
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_iterator1 const_iterator_type;
         typedef typename M::iterator1 iterator_type;
 #else
@@ -600,11 +602,11 @@ namespace boost { namespace numeric { namespace ublas {
             return j_;
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -1015,7 +1017,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         size_type j_;
         static matrix_type nil_;
     };
@@ -1032,7 +1034,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M>
     BOOST_UBLAS_INLINE
-    const matrix_column<const M> column (const M &data, std::size_t j) {
+    matrix_column<const M> column (const M &data, std::size_t j) {
         return matrix_column<const M> (data, j);
     }
 #endif
@@ -1050,7 +1052,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -1065,14 +1067,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const matrix_vector_range<matrix_type> > const_closure_type;
-        typedef vector_reference<matrix_vector_range<matrix_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_vector_range<matrix_type> const_closure_type;
         typedef matrix_vector_range<matrix_type> closure_type;
-#endif
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         typedef range<>::const_iterator const_iterator1_type;
         typedef range<>::const_iterator iterator1_type;
@@ -1089,7 +1092,7 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Construction and destruction
         BOOST_UBLAS_INLINE
-        matrix_vector_range (): 
+        matrix_vector_range ():
             data_ (nil_), r1_ (), r2_ () {}
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
@@ -1121,11 +1124,11 @@ namespace boost { namespace numeric { namespace ublas {
             return BOOST_UBLAS_SAME (r1_.size (), r2_.size ());
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -1528,7 +1531,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         range<> r1_;
         range<> r2_;
@@ -1555,7 +1558,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -1570,14 +1573,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const matrix_vector_slice<matrix_type> > const_closure_type;
-        typedef vector_reference<matrix_vector_slice<matrix_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_vector_slice<matrix_type> const_closure_type;
         typedef matrix_vector_slice<matrix_type> closure_type;
-#endif
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         typedef slice<>::const_iterator const_iterator1_type;
         typedef slice<>::const_iterator iterator1_type;
@@ -1629,14 +1633,14 @@ namespace boost { namespace numeric { namespace ublas {
         // Accessors
         BOOST_UBLAS_INLINE
         size_type size () const {
-            return BOOST_UBLAS_SAME (s1_.size (), s2_.size ()); 
+            return BOOST_UBLAS_SAME (s1_.size (), s2_.size ());
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -2051,7 +2055,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         slice<> s1_;
         slice<> s2_;
@@ -2080,7 +2084,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -2095,14 +2099,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const matrix_vector_indirect<matrix_type, indirect_array_type> > const_closure_type;
-        typedef vector_reference<matrix_vector_indirect<matrix_type, indirect_array_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_vector_indirect<matrix_type, indirect_array_type> const_closure_type;
         typedef matrix_vector_indirect<matrix_type, indirect_array_type> closure_type;
-#endif
         typedef typename IA::const_iterator const_iterator1_type;
         typedef typename IA::const_iterator iterator1_type;
         typedef typename IA::const_iterator const_iterator2_type;
@@ -2130,11 +2135,11 @@ namespace boost { namespace numeric { namespace ublas {
             return BOOST_UBLAS_SAME (ia1_.size (), ia2_.size ());
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
         BOOST_UBLAS_INLINE
@@ -2553,7 +2558,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         indirect_array_type ia1_;
         indirect_array_type ia2_;
         static matrix_type nil_;
@@ -2575,7 +2580,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -2590,15 +2595,16 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const matrix_const_reference<const matrix_range<matrix_type> > const_closure_type;
-        typedef matrix_reference<matrix_range<matrix_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_range<matrix_type> const_closure_type;
         typedef matrix_range<matrix_type> closure_type;
-#endif
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_iterator1 const_iterator1_type;
         typedef typename M::iterator1 iterator1_type;
         typedef typename M::const_iterator2 const_iterator2_type;
@@ -2631,6 +2637,15 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (r2_.start () <= data_.size2 () &&
                                r2_.start () + r2_.size () <= data_.size2 (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        matrix_range (const matrix_closure_type &data, const range<> &r1, const range<> &r2, int):
+            data_ (data), r1_ (r1.preprocess (data.size1 ())), r2_ (r2.preprocess (data.size2 ())) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (r1_.start () <= data_.size1 () &&
+                               r1_.start () + r1_.size () <= data_.size1 (), bad_index ());
+            BOOST_UBLAS_CHECK (r2_.start () <= data_.size2 () &&
+                               r2_.start () + r2_.size () <= data_.size2 (), bad_index ());
+        }
 #else
         BOOST_UBLAS_INLINE
         matrix_range (matrix_type &data, const range &r1, const range &r2):
@@ -2641,31 +2656,40 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (r2_.start () <= data_.size2 () &&
                                r2_.start () + r2_.size () <= data_.size2 (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        matrix_range (const matrix_closure_type &data, const range &r1, const range &r2, int):
+            data_ (data), r1_ (r1), r2_ (r2) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (r1_.start () <= data_.size1 () &&
+                               r1_.start () + r1_.size () <= data_.size1 (), bad_index ());
+            BOOST_UBLAS_CHECK (r2_.start () <= data_.size2 () &&
+                               r2_.start () + r2_.size () <= data_.size2 (), bad_index ());
+        }
 #endif
 
         // Accessors
         BOOST_UBLAS_INLINE
         size_type start1 () const {
-            return r1_.start (); 
+            return r1_.start ();
         }
         BOOST_UBLAS_INLINE
-        size_type size1 () const { 
+        size_type size1 () const {
             return r1_.size ();
         }
         BOOST_UBLAS_INLINE
-        size_type start2() const { 
-            return r2_.start (); 
+        size_type start2() const {
+            return r2_.start ();
         }
         BOOST_UBLAS_INLINE
         size_type size2 () const { 
             return r2_.size ();
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -2698,23 +2722,23 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
         matrix_range<matrix_type> project (const range<> &r1, const range<> &r2) const {
-            return matrix_range<matrix_type>  (data_, r1_.compose (r1), r2_.compose (r2));
+            return matrix_range<matrix_type>  (data_, r1_.compose (r1), r2_.compose (r2), 0);
         }
 #else
         BOOST_UBLAS_INLINE
         matrix_range<matrix_type> project (const range &r1, const range &r2) const {
-            return matrix_range<matrix_type>  (data_, r1_.compose (r1), r2_.compose (r2));
+            return matrix_range<matrix_type>  (data_, r1_.compose (r1), r2_.compose (r2), 0);
         }
 #endif
 
         // Assignment
         BOOST_UBLAS_INLINE
-        matrix_range &operator = (const matrix_range &mr) { 
-            matrix_assign (scalar_assign<value_type, value_type> (), *this, mr); 
+        matrix_range &operator = (const matrix_range &mr) {
+            matrix_assign (scalar_assign<value_type, value_type> (), *this, mr);
             return *this;
         }
         BOOST_UBLAS_INLINE
-        matrix_range &assign_temporary (matrix_range &mr) { 
+        matrix_range &assign_temporary (matrix_range &mr) {
             return *this = mr;
         }
         template<class AE>
@@ -2744,7 +2768,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class AE>
         BOOST_UBLAS_INLINE
         matrix_range& operator -= (const matrix_expression<AE> &ae) {
-            matrix_assign (scalar_assign<value_type, value_type> (), *this, matrix<value_type> (*this - ae)); 
+            matrix_assign (scalar_assign<value_type, value_type> (), *this, matrix<value_type> (*this - ae));
             return *this;
         }
         template<class AE>
@@ -2894,8 +2918,13 @@ namespace boost { namespace numeric { namespace ublas {
             public iterator_base_traits<typename const_iterator1_type::iterator_category>::template
                         iterator_base<const_iterator1, value_type>::type {
 #else
+#ifndef BOOST_MSVC
             public random_access_iterator_base<typename const_iterator1_type::iterator_category,
                                                const_iterator1, value_type> {
+#else
+            public random_access_iterator_base<typename M::const_iterator1::iterator_category,
+                                               const_iterator1, value_type> {
+#endif
 #endif                                               
         public:
             typedef typename const_iterator1_type::iterator_category iterator_category;
@@ -3019,8 +3048,13 @@ namespace boost { namespace numeric { namespace ublas {
             public iterator_base_traits<typename iterator1_type::iterator_category>::template
                         iterator_base<iterator1, value_type>::type {
 #else
+#ifndef BOOST_MSVC
             public random_access_iterator_base<typename iterator1_type::iterator_category,
                                                iterator1, value_type> {
+#else
+            public random_access_iterator_base<typename M::iterator1::iterator_category,
+                                               iterator1, value_type> {
+#endif
 #endif                                               
         public:
             typedef typename iterator1_type::iterator_category iterator_category;
@@ -3143,8 +3177,13 @@ namespace boost { namespace numeric { namespace ublas {
             public iterator_base_traits<typename const_iterator2_type::iterator_category>::template
                         iterator_base<const_iterator2, value_type>::type {
 #else
+#ifndef BOOST_MSVC
             public random_access_iterator_base<typename const_iterator2_type::iterator_category,
                                                const_iterator2, value_type> {
+#else
+            public random_access_iterator_base<typename M::const_iterator2::iterator_category,
+                                               const_iterator2, value_type> {
+#endif
 #endif                                               
         public:
             typedef typename const_iterator2_type::iterator_category iterator_category;
@@ -3268,8 +3307,13 @@ namespace boost { namespace numeric { namespace ublas {
             public iterator_base_traits<typename iterator2_type::iterator_category>::template
                         iterator_base<iterator2, value_type>::type {
 #else
+#ifndef BOOST_MSVC
             public random_access_iterator_base<typename iterator2_type::iterator_category,
                                                iterator2, value_type> {
+#else
+            public random_access_iterator_base<typename M::iterator2::iterator_category,
+                                               iterator2, value_type> {
+#endif
 #endif
         public:
             typedef typename iterator2_type::iterator_category iterator_category;
@@ -3424,7 +3468,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         range<> r1_;
         range<> r2_;
@@ -3448,7 +3492,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M>
     BOOST_UBLAS_INLINE
-    const matrix_range<const M> project (const M &data, const range<> &r1, const range<> &r2) {
+    matrix_range<const M> project (const M &data, const range<> &r1, const range<> &r2) {
         return matrix_range<const M> (data, r1, r2);
     }
     template<class M>
@@ -3466,7 +3510,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M>
     BOOST_UBLAS_INLINE
-    const matrix_range<const M> project (const M &data, const range &r1, const range &r2) {
+    matrix_range<const M> project (const M &data, const range &r1, const range &r2) {
         return matrix_range<const M> (data, r1, r2);
     }
     template<class M>
@@ -3490,7 +3534,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -3505,14 +3549,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const matrix_const_reference<const matrix_slice<matrix_type> > const_closure_type;
-        typedef matrix_reference<matrix_slice<matrix_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_slice<matrix_type> const_closure_type;
         typedef matrix_slice<matrix_type> closure_type;
-#endif
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         typedef slice<>::const_iterator const_iterator1_type;
         typedef slice<>::const_iterator iterator1_type;
@@ -3542,6 +3587,15 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (s2_.start () <= data_.size2 () &&
                                s2_.start () + s2_.stride () * (s2_.size () - (s2_.size () > 0)) <= data_.size2 (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        matrix_slice (const matrix_closure_type &data, const slice<> &s1, const slice<> &s2, int):
+            data_ (data), s1_ (s1.preprocess (data.size1 ())), s2_ (s2.preprocess (data.size2 ())) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (s1_.start () <= data_.size1 () &&
+                               s1_.start () + s1_.stride () * (s1_.size () - (s1_.size () > 0)) <= data_.size1 (), bad_index ());
+            BOOST_UBLAS_CHECK (s2_.start () <= data_.size2 () &&
+                               s2_.start () + s2_.stride () * (s2_.size () - (s2_.size () > 0)) <= data_.size2 (), bad_index ());
+        }
 #else
         BOOST_UBLAS_INLINE
         matrix_slice (matrix_type &data, const slice &s1, const slice &s2):
@@ -3552,6 +3606,15 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (s2_.start () <= data_.size2 () &&
                                s2_.start () + s2_.stride () * (s2_.size () - (s2_.size () > 0)) <= data_.size2 (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        matrix_slice (const matrix_closure_type &data, const slice &s1, const slice &s2, int):
+            data_ (data), s1_ (s1), s2_ (s2) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (s1_.start () <= data_.size1 () &&
+                               s1_.start () + s1_.stride () * (s1_.size () - (s1_.size () > 0)) <= data_.size1 (), bad_index ());
+            BOOST_UBLAS_CHECK (s2_.start () <= data_.size2 () &&
+                               s2_.start () + s2_.stride () * (s2_.size () - (s2_.size () > 0)) <= data_.size2 (), bad_index ());
+        }
 #endif
 
         // Accessors
@@ -3580,11 +3643,11 @@ namespace boost { namespace numeric { namespace ublas {
             return s2_.size ();
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -3617,20 +3680,20 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
         matrix_slice<matrix_type> project (const range<> &r1, const range<> &r2) const {
-            return matrix_slice<matrix_type>  (data_, s1_.compose (r1), s2_.compose (r2));
+            return matrix_slice<matrix_type>  (data_, s1_.compose (r1), s2_.compose (r2), 0);
         }
         BOOST_UBLAS_INLINE
         matrix_slice<matrix_type> project (const slice<> &s1, const slice<> &s2) const {
-            return matrix_slice<matrix_type>  (data_, s1_.compose (s1), s2_.compose (s2));
+            return matrix_slice<matrix_type>  (data_, s1_.compose (s1), s2_.compose (s2), 0);
         }
 #else
         BOOST_UBLAS_INLINE
         matrix_slice<matrix_type> project (const range &r1, const range &r2) const {
-            return matrix_slice<matrix_type>  (data_, s1_.compose (r1), s2_.compose (r2));
+            return matrix_slice<matrix_type>  (data_, s1_.compose (r1), s2_.compose (r2), 0);
         }
         BOOST_UBLAS_INLINE
         matrix_slice<matrix_type> project (const slice &s1, const slice &s2) const {
-            return matrix_slice<matrix_type>  (data_, s1_.compose (s1), s2_.compose (s2));
+            return matrix_slice<matrix_type>  (data_, s1_.compose (s1), s2_.compose (s2), 0);
         }
 #endif
 
@@ -3671,7 +3734,7 @@ namespace boost { namespace numeric { namespace ublas {
         template<class AE>
         BOOST_UBLAS_INLINE
         matrix_slice& operator -= (const matrix_expression<AE> &ae) {
-            matrix_assign (scalar_assign<value_type, value_type> (), *this, matrix<value_type> (*this - ae)); 
+            matrix_assign (scalar_assign<value_type, value_type> (), *this, matrix<value_type> (*this - ae));
             return *this;
         }
         template<class AE>
@@ -4343,7 +4406,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         slice<> s1_;
         slice<> s2_;
@@ -4374,7 +4437,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M>
     BOOST_UBLAS_INLINE
-    const matrix_slice<const M> project (const M &data, const slice<> &s1, const slice<> &s2) {
+    matrix_slice<const M> project (const M &data, const slice<> &s1, const slice<> &s2) {
         return matrix_slice<const M> (data, s1, s2);
     }
     template<class M>
@@ -4399,7 +4462,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M>
     BOOST_UBLAS_INLINE
-    const matrix_slice<const M> project (const M &data, const slice &s1, const slice &s2) {
+    matrix_slice<const M> project (const M &data, const slice &s1, const slice &s2) {
         return matrix_slice<const M> (data, s1, s2);
     }
     template<class M>
@@ -4427,7 +4490,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -4442,14 +4505,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const matrix_const_reference<const matrix_indirect<matrix_type, indirect_array_type> > const_closure_type;
-        typedef matrix_reference<matrix_indirect<matrix_type, indirect_array_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const matrix_indirect<matrix_type, indirect_array_type> const_closure_type;
         typedef matrix_indirect<matrix_type, indirect_array_type> closure_type;
-#endif
         typedef typename IA::const_iterator const_iterator1_type;
         typedef typename IA::const_iterator iterator1_type;
         typedef typename IA::const_iterator const_iterator2_type;
@@ -4469,10 +4533,16 @@ namespace boost { namespace numeric { namespace ublas {
         BOOST_UBLAS_INLINE
         matrix_indirect (matrix_type &data, const indirect_array_type &ia1, const indirect_array_type &ia2):
             data_ (data), ia1_ (ia1.preprocess (data.size1 ())), ia2_ (ia2.preprocess (data.size2 ())) {}
+        BOOST_UBLAS_INLINE
+        matrix_indirect (const matrix_closure_type &data, const indirect_array_type &ia1, const indirect_array_type &ia2, int):
+            data_ (data), ia1_ (ia1.preprocess (data.size1 ())), ia2_ (ia2.preprocess (data.size2 ())) {}
 #else
         BOOST_UBLAS_INLINE
         matrix_indirect (matrix_type &data, const indirect_array_type &ia1, const indirect_array_type &ia2):
             data_ (data), ia1_ (ia1), ia2_ (ia2) {}
+        BOOST_UBLAS_INLINE
+        matrix_indirect (const matrix_closure_type &data, const indirect_array_type &ia1, const indirect_array_type &ia2, int):
+            data_ (data), ia1_ (ia1), ia2_ (ia2) {}
 #endif
 
         // Accessors
@@ -4485,11 +4555,11 @@ namespace boost { namespace numeric { namespace ublas {
             return ia2_.size ();
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
         BOOST_UBLAS_INLINE
@@ -4538,25 +4608,25 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
         matrix_indirect<matrix_type, indirect_array_type> project (const range<> &r1, const range<> &r2) const {
-            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (r1), ia2_.compose (r2));
+            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (r1), ia2_.compose (r2), 0);
         }
         BOOST_UBLAS_INLINE
         matrix_indirect<matrix_type, indirect_array_type> project (const slice<> &s1, const slice<> &s2) const {
-            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (s1), ia2_.compose (s2));
+            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (s1), ia2_.compose (s2), 0);
         }
 #else
         BOOST_UBLAS_INLINE
         matrix_indirect<matrix_type, indirect_array_type> project (const range &r1, const range &r2) const {
-            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (r1), ia2_.compose (r2));
+            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (r1), ia2_.compose (r2), 0);
         }
         BOOST_UBLAS_INLINE
         matrix_indirect<matrix_type, indirect_array_type> project (const slice &s1, const slice &s2) const {
-            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (s1), ia2_.compose (s2));
+            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (s1), ia2_.compose (s2), 0);
         }
 #endif
         BOOST_UBLAS_INLINE
         matrix_indirect<matrix_type, indirect_array_type> project (const indirect_array_type &ia1, const indirect_array_type &ia2) const {
-            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (ia1), ia2_.compose (ia2));
+            return matrix_indirect<matrix_type, indirect_array_type> (data_, ia1_.compose (ia1), ia2_.compose (ia2), 0);
         }
 
         // Assignment
@@ -5268,7 +5338,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         indirect_array_type ia1_;
         indirect_array_type ia2_;
         static matrix_type nil_;
@@ -5314,7 +5384,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     // template<class M, class IA>
     // BOOST_UBLAS_INLINE
-    // const matrix_indirect<const M, IA> project (const M &data, const IA &ia1, const IA &ia2) {
+    // matrix_indirect<const M, IA> project (const M &data, const IA &ia1, const IA &ia2) {
     //     return matrix_indirect<const M, IA> (data, ia1, ia2);
     // }
     // template<class M, class IA>
@@ -5331,7 +5401,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class M, class A>
     BOOST_UBLAS_INLINE
-    const matrix_indirect<const M, indirect_array<A> > project (const M &data, const indirect_array<A> &ia1, const indirect_array<A> &ia2) {
+    matrix_indirect<const M, indirect_array<A> > project (const M &data, const indirect_array<A> &ia1, const indirect_array<A> &ia2) {
         return matrix_indirect<const M, indirect_array<A> > (data, ia1, ia2);
     }
     template<class M, class A>
diff --git a/include/boost/numeric/ublas/matrix_sparse.hpp b/include/boost/numeric/ublas/matrix_sparse.hpp
index 57ac8654..0576b168 100644
--- a/include/boost/numeric/ublas/matrix_sparse.hpp
+++ b/include/boost/numeric/ublas/matrix_sparse.hpp
@@ -185,7 +185,7 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
         typedef T &reference;
 #elif defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
-        typedef sparse_matrix_element<sparse_matrix> reference;
+        typedef sparse_matrix_element<sparse_matrix<T, F, A> > reference;
 #elif defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE)
         typedef typename map_traits<A>::reference reference;
 #endif
@@ -779,7 +779,7 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
+#if defined (BOOST_MSVC_STD_ITERATOR) || (! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE))
                     return (*it_).second;
 #elif defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return reference ((*this) (), &(*it_).second, index1 (), index2 ());
@@ -787,7 +787,7 @@ namespace boost { namespace numeric { namespace ublas {
                     return map_traits<array_type>::make_reference ((*this) ().data (), it_);
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -1058,7 +1058,7 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
+#if defined (BOOST_MSVC_STD_ITERATOR) || (! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE))
                     return (*it_).second;
 #elif defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return reference ((*this) (), &(*it_).second, index1 (), index2 ());
@@ -1066,7 +1066,7 @@ namespace boost { namespace numeric { namespace ublas {
                     return map_traits<array_type>::make_reference ((*this) ().data (), it_);
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -1215,7 +1215,7 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
         typedef T &reference;
 #elif defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
-        typedef sparse_matrix_element<sparse_vector_of_sparse_vector> reference;
+        typedef sparse_matrix_element<sparse_vector_of_sparse_vector<T, F, A> > reference;
 #elif defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE)
         typedef typename map_traits<typename A::data_value_type>::reference reference;
 #endif
@@ -1850,7 +1850,7 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
+#if defined (BOOST_MSVC_STD_ITERATOR) || (! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE))
                     return (*it_).second;
 #elif defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return reference ((*this) (), &(*it_).second, index1 (), index2 ());
@@ -1858,7 +1858,7 @@ namespace boost { namespace numeric { namespace ublas {
                     return map_traits<vector_data_value_type>::make_reference ((*itv_).second, it_);
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -2161,7 +2161,7 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
+#if defined (BOOST_MSVC_STD_ITERATOR) || (! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE))
                     return (*it_).second;
 #elif defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return reference ((*this) (), &(*it_).second, index1 (), index2 ());
@@ -2169,7 +2169,7 @@ namespace boost { namespace numeric { namespace ublas {
                     return map_traits<vector_data_value_type>::make_reference ((*itv_).second, it_);
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -2304,12 +2304,12 @@ namespace boost { namespace numeric { namespace ublas {
 
     // Array based sparse matrix class
     // Thanks to Kresimir Fresl for extending this to cover different index bases.
-    template<class T, class F, class IA, class TA, std::size_t IB>
+    template<class T, class F, std::size_t IB, class IA, class TA>
     class compressed_matrix:
-        public matrix_expression<compressed_matrix<T, F, IA, TA, IB> > {
+        public matrix_expression<compressed_matrix<T, F, IB, IA, TA> > {
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
-        BOOST_UBLAS_USING matrix_expression<compressed_matrix<T, F, IA, TA, IB> >::operator ();
+        BOOST_UBLAS_USING matrix_expression<compressed_matrix<T, F, IB, IA, TA> >::operator ();
 #endif
         typedef std::size_t size_type;
         typedef std::ptrdiff_t difference_type;
@@ -2319,15 +2319,15 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
         typedef T &reference;
 #else
-        typedef sparse_matrix_element<compressed_matrix> reference;
+        typedef sparse_matrix_element<compressed_matrix<T, F, IB, IA, TA> > reference;
 #endif
         typedef const T *const_pointer;
         typedef T *pointer;
         typedef IA index_array_type;
         typedef TA value_array_type;
         typedef F functor_type;
-        typedef const compressed_matrix<T, F, IA, TA, IB> const_self_type;
-        typedef compressed_matrix<T, F, IA, TA, IB> self_type;
+        typedef const compressed_matrix<T, F, IB, IA, TA> const_self_type;
+        typedef compressed_matrix<T, F, IB, IA, TA> self_type;
         typedef const matrix_const_reference<const_self_type> const_closure_type;
         typedef matrix_reference<self_type> closure_type;
         typedef typename IA::const_iterator vector_const_iterator_type;
@@ -2393,6 +2393,10 @@ namespace boost { namespace numeric { namespace ublas {
             return filled2_;
         }
         BOOST_UBLAS_INLINE
+        static size_type index_base () {
+            return index_base_;
+        }
+        BOOST_UBLAS_INLINE
         const index_array_type &index1_data () const {
             return index1_data_;
         }
@@ -2416,10 +2420,6 @@ namespace boost { namespace numeric { namespace ublas {
         value_array_type &value_data () {
             return value_data_;
         }
-        BOOST_UBLAS_INLINE
-        size_type index_base () {
-            return index_base_;
-        }
 
         // Resizing
         BOOST_UBLAS_INLINE
@@ -2484,7 +2484,7 @@ namespace boost { namespace numeric { namespace ublas {
                 insert (i, j, value_type ());
             vector_const_iterator_type itv (index1_data ().begin () + element1);
             const_iterator_type it_begin (index2_data ().begin () + zero_based (*itv));
-            const_iterator_type it_end (index2_data ().begin () + zero_based (*(itv + 1));
+            const_iterator_type it_end (index2_data ().begin () + zero_based (*(itv + 1)));
             const_iterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
             if (it == it_end || *it != k_based (element2)) {
                 insert (i, j, value_type ());
@@ -3100,13 +3100,13 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
+#if defined (BOOST_MSVC_STD_ITERATOR) || ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return (*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()];
 #else
                     return reference ((*this) (), &(*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()], index1 (), index2 ());
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -3389,13 +3389,13 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
+#if defined (BOOST_MSVC_STD_ITERATOR) || ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return (*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()];
 #else
                     return reference ((*this) (), &(*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()], index1 (), index2 ());
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -3544,17 +3544,17 @@ namespace boost { namespace numeric { namespace ublas {
         friend class const_iterator2;
     };
 
-    template<class T, class F, class IA, class TA, std::size_t IB>
-    typename compressed_matrix<T, F, IA, TA, IB>::value_type compressed_matrix<T, F, IA, TA, IB>::zero_ = 0;
+    template<class T, class F, std::size_t IB, class IA, class TA>
+    typename compressed_matrix<T, F, IB, IA, TA>::value_type compressed_matrix<T, F, IB, IA, TA>::zero_ = 0;
 
     // Array based sparse matrix class
     // Thanks to Kresimir Fresl for extending this to cover different index bases.
-    template<class T, class F, class IA, class TA, std::size_t IB>
+    template<class T, class F, std::size_t IB, class IA, class TA>
     class coordinate_matrix:
-        public matrix_expression<coordinate_matrix<T, F, IA, TA, IB> > {
+        public matrix_expression<coordinate_matrix<T, F, IB, IA, TA> > {
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
-        BOOST_UBLAS_USING matrix_expression<coordinate_matrix<T, F, IA, TA, IB> >::operator ();
+        BOOST_UBLAS_USING matrix_expression<coordinate_matrix<T, F, IB, IA, TA> >::operator ();
 #endif
         typedef std::size_t size_type;
         typedef std::ptrdiff_t difference_type;
@@ -3564,15 +3564,15 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
         typedef T &reference;
 #else
-        typedef sparse_matrix_element<coordinate_matrix> reference;
+        typedef sparse_matrix_element<coordinate_matrix<T, F, IB, IA, TA> > reference;
 #endif
         typedef const T *const_pointer;
         typedef T *pointer;
         typedef IA index_array_type;
         typedef TA value_array_type;
         typedef F functor_type;
-        typedef const coordinate_matrix<T, F, IA, TA, IB> const_self_type;
-        typedef coordinate_matrix<T, F, IA, TA, IB> self_type;
+        typedef const coordinate_matrix<T, F, IB, IA, TA> const_self_type;
+        typedef coordinate_matrix<T, F, IB, IA, TA> self_type;
         typedef const matrix_const_reference<const_self_type> const_closure_type;
         typedef matrix_reference<self_type> closure_type;
         typedef typename IA::const_iterator vector_const_iterator_type;
@@ -3632,6 +3632,10 @@ namespace boost { namespace numeric { namespace ublas {
             return filled_;
         }
         BOOST_UBLAS_INLINE
+        static size_type index_base () {
+            return index_base_;
+        }
+        BOOST_UBLAS_INLINE
         const index_array_type &index1_data () const {
             return index1_data_;
         }
@@ -3655,10 +3659,6 @@ namespace boost { namespace numeric { namespace ublas {
         value_array_type &value_data () {
             return value_data_;
         }
-        BOOST_UBLAS_INLINE
-        size_type index_base () {
-            return index_base_;
-        }
 
         // Resizing
         BOOST_UBLAS_INLINE
@@ -3878,10 +3878,10 @@ namespace boost { namespace numeric { namespace ublas {
         BOOST_UBLAS_INLINE
         void sort () const {
             if (! sorted_) {
-                typedef std::pair<typename index_array_type::value_type,
-                                  typename index_array_type::value_type> pair_type;
+                typedef std::pair<BOOST_UBLAS_TYPENAME index_array_type::value_type,
+                                  BOOST_UBLAS_TYPENAME index_array_type::value_type> pair_type;
                 typedef std::pair<pair_type,
-                                  typename value_array_type::value_type> triple_type;
+                                  BOOST_UBLAS_TYPENAME value_array_type::value_type> triple_type;
                 unbounded_array<triple_type> a (filled_);
                 size_type i;
                 for (i = 0; i < filled_; ++ i)
@@ -4337,13 +4337,13 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
+#if defined (BOOST_MSVC_STD_ITERATOR) || ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return (*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()];
 #else
                     return reference ((*this) (), &(*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()], index1 (), index2 ());
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -4626,13 +4626,13 @@ namespace boost { namespace numeric { namespace ublas {
                 if (rank_ == 1) {
                     BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
                     BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
-#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
+#if defined (BOOST_MSVC_STD_ITERATOR) || ! defined (BOOST_UBLAS_STRICT_MATRIX_SPARSE)
                     return (*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()];
 #else
                     return reference ((*this) (), &(*this) ().value_data () [it_ - (*this) ().index2_data ().begin ()], index1 (), index2 ());
 #endif
                 } else {
-                    return (*this) () (i_, j_);
+                    return reference ((*this) () (i_, j_));
                 }
             }
 
@@ -4781,8 +4781,8 @@ namespace boost { namespace numeric { namespace ublas {
         friend class const_iterator2;
     };
 
-    template<class T, class F, class IA, class TA, std::size_t IB>
-    typename coordinate_matrix<T, F, IA, TA, IB>::value_type coordinate_matrix<T, F, IA, TA, IB>::zero_ = 0;
+    template<class T, class F, std::size_t IB, class IA, class TA>
+    typename coordinate_matrix<T, F, IB, IA, TA>::value_type coordinate_matrix<T, F, IB, IA, TA>::zero_ = 0;
 
 }}}
 
diff --git a/include/boost/numeric/ublas/storage_sparse.hpp b/include/boost/numeric/ublas/storage_sparse.hpp
index d797e88d..cc002e0d 100644
--- a/include/boost/numeric/ublas/storage_sparse.hpp
+++ b/include/boost/numeric/ublas/storage_sparse.hpp
@@ -152,11 +152,13 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Element access
         // FIXME: GCC 3.1 warn's, if enabled
-        // BOOST_UBLAS_INLINE
-        // const sparse_storage_element_traits<data_value_type>::data_const_reference
-        // operator [] (typename sparse_storage_element_traits<data_value_type>::index_type i) const {
-        //     return d_ [i];
-        // }
+#ifndef __GNUC__
+        BOOST_UBLAS_INLINE
+        typename sparse_storage_element_traits<data_value_type>::data_const_reference
+        operator [] (typename sparse_storage_element_traits<data_value_type>::index_type i) const {
+            return d_ [i];
+        }
+#endif
         BOOST_UBLAS_INLINE
         typename sparse_storage_element_traits<data_value_type>::data_reference
         operator [] (typename sparse_storage_element_traits<data_value_type>::index_type i) {
@@ -220,12 +222,14 @@ namespace boost { namespace numeric { namespace ublas {
             return *this;
         }
 
-        // FIXME: GCC 3.1 warn's, if enabled
         // Conversion
+        // FIXME: GCC 3.1 warn's, if enabled
+#ifndef __GNUC__
         BOOST_UBLAS_INLINE
-        // operator const data_const_reference () const {
-        //     return d_;
-        // }
+        operator data_const_reference () const {
+            return d_;
+        }
+#endif
         BOOST_UBLAS_INLINE
         operator data_reference () {
             return d_;
diff --git a/include/boost/numeric/ublas/symmetric.hpp b/include/boost/numeric/ublas/symmetric.hpp
index 23a1be44..c213383b 100644
--- a/include/boost/numeric/ublas/symmetric.hpp
+++ b/include/boost/numeric/ublas/symmetric.hpp
@@ -913,7 +913,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -927,11 +927,18 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename detail::ct_if<boost::is_const<M>::value,
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
+#endif
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
         typedef const symmetric_adaptor<M, F> const_self_type;
         typedef symmetric_adaptor<M, F> self_type;
-        typedef const matrix_const_reference<const_self_type> const_closure_type;
-        typedef matrix_reference<self_type> closure_type;
+        typedef const_self_type const_closure_type;
+        typedef self_type closure_type;
         typedef typename storage_restrict_traits<typename M::storage_category,
                                                  packed_proxy_tag>::storage_category storage_category;
         typedef typename F::packed_category packed_category;
@@ -964,11 +971,11 @@ namespace boost { namespace numeric { namespace ublas {
             return data_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -1675,7 +1682,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         static matrix_type nil_;
     };
 
diff --git a/include/boost/numeric/ublas/triangular.hpp b/include/boost/numeric/ublas/triangular.hpp
index eda6880b..2ab6cf3f 100644
--- a/include/boost/numeric/ublas/triangular.hpp
+++ b/include/boost/numeric/ublas/triangular.hpp
@@ -56,7 +56,7 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Construction and destruction
         BOOST_UBLAS_INLINE
-        triangular_matrix (): 
+        triangular_matrix ():
             size1_ (0), size2_ (0),
             data_ (0) {}
         BOOST_UBLAS_INLINE
@@ -912,7 +912,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename M::size_type size_type;
         typedef typename M::difference_type difference_type;
         typedef typename M::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename M::const_reference const_reference;
         typedef typename M::reference reference;
         typedef typename M::const_pointer const_pointer;
@@ -926,11 +926,18 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename detail::ct_if<boost::is_const<M>::value,
                                        typename M::const_pointer,
                                        typename M::pointer>::type pointer;
+#endif
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename M::closure_type matrix_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<M>::value,
+                                       typename M::const_closure_type,
+                                       typename M::closure_type>::type matrix_closure_type;
 #endif
         typedef const triangular_adaptor<M, F> const_self_type;
         typedef triangular_adaptor<M, F> self_type;
-        typedef const matrix_const_reference<const_self_type> const_closure_type;
-        typedef matrix_reference<self_type> closure_type;
+        typedef const_self_type const_closure_type;
+        typedef self_type closure_type;
         typedef typename storage_restrict_traits<typename M::storage_category,
                                                  packed_proxy_tag>::storage_category storage_category;
         typedef typename F::packed_category packed_category;
@@ -938,13 +945,13 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Construction and destruction
         BOOST_UBLAS_INLINE
-        triangular_adaptor (): 
+        triangular_adaptor ():
             data_ (nil_) {}
         BOOST_UBLAS_INLINE
-        triangular_adaptor (matrix_type &data): 
+        triangular_adaptor (matrix_type &data):
             data_ (data) {}
         BOOST_UBLAS_INLINE
-        triangular_adaptor (const triangular_adaptor &m): 
+        triangular_adaptor (const triangular_adaptor &m):
             data_ (m.data_) {}
 
         // Accessors
@@ -957,11 +964,11 @@ namespace boost { namespace numeric { namespace ublas {
             return data_.size2 ();
         }
         BOOST_UBLAS_INLINE
-        const_matrix_type &data () const {
+        const matrix_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        matrix_type &data () {
+        matrix_closure_type &data () {
             return data_;
         }
 
@@ -1686,7 +1693,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        matrix_type &data_;
+        matrix_closure_type data_;
         static matrix_type nil_;
         static value_type zero_;
         static value_type one_;
@@ -1851,6 +1858,13 @@ namespace boost { namespace numeric { namespace ublas {
         inplace_solve (e1, e2, upper_tag (), vector_tag (), dispatch_category ());
     }
 
+    template<class E1, class E2, class C>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1,
+                        vector_expression<E2> &e2,
+                        C) {
+        inplace_solve (e1, e2, C (), vector_tag ());
+    }
     template<class E1, class E2, class C>
     BOOST_UBLAS_INLINE
     typename matrix_vector_solve_traits<E1, E2>::result_type
@@ -2006,6 +2020,13 @@ namespace boost { namespace numeric { namespace ublas {
         inplace_solve (e1, e2, upper_tag (), vector_tag (), dispatch_category ());
     }
 
+    template<class E1, class E2, class C>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (vector_expression<E1> &e1,
+                        const matrix_expression<E2> &e2,
+                        C) {
+        inplace_solve (e1, e2, vector_tag (), C ());
+    }
     template<class E1, class E2, class C>
     BOOST_UBLAS_INLINE
     typename matrix_vector_solve_traits<E1, E2>::result_type
@@ -2187,6 +2208,13 @@ namespace boost { namespace numeric { namespace ublas {
         inplace_solve (e1, e2, upper_tag (), matrix_tag (), dispatch_category ());
     }
 
+    template<class E1, class E2, class C>
+    BOOST_UBLAS_INLINE
+    void inplace_solve (const matrix_expression<E1> &e1,
+                        matrix_expression<E2> &e2,
+                        C) {
+        inplace_solve (e1, e2, C (), matrix_tag ());
+    }
     template<class E1, class E2, class C>
     BOOST_UBLAS_INLINE
     typename matrix_matrix_solve_traits<E1, E2>::result_type
diff --git a/include/boost/numeric/ublas/vector_assign.hpp b/include/boost/numeric/ublas/vector_assign.hpp
index c380c2b4..fc66ca53 100644
--- a/include/boost/numeric/ublas/vector_assign.hpp
+++ b/include/boost/numeric/ublas/vector_assign.hpp
@@ -347,6 +347,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef F functor_type;
         typedef typename V::size_type size_type;
         typedef typename V::value_type value_type;
+        typedef typename V::reference reference;
 #ifdef BOOST_UBLAS_TYPE_CHECK
         vector<value_type> cv (v.size ());
         indexing_vector_assign (scalar_assign<value_type, value_type> (), cv, v);
@@ -365,6 +366,7 @@ namespace boost { namespace numeric { namespace ublas {
                 functor_type () (*it, value_type ());
                 ++ it;
             } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 // Sparse proxies don't need to be conformant.
                 // Thanks to Michael Stevens for suggesting this.
                 size_type index (ite.index ());
@@ -372,8 +374,10 @@ namespace boost { namespace numeric { namespace ublas {
                 restart (v, index, it, it_end);
                 // The proxies could reference the same container.
                 restart (e, index, ite, ite_end);
+#endif
             }
         }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
         while (ite != ite_end) {
             // Sparse proxies don't need to be conformant.
             // Thanks to Michael Stevens for suggesting this.
@@ -382,6 +386,7 @@ namespace boost { namespace numeric { namespace ublas {
             // The proxies could reference the same container.
             restart (e, index, ite, ite_end);
         }
+#endif
         while (it != it_end) {
             functor_type () (*it, value_type ());
             ++ it;
@@ -462,6 +467,7 @@ namespace boost { namespace numeric { namespace ublas {
                 functor_type () (*it, *ite);
                 ++ it, ++ ite;
             } else if (compare < 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 // Sparse proxies don't need to be conformant.
                 // Thanks to Michael Stevens for suggesting this.
                 size_type index (it.index ());
@@ -469,7 +475,9 @@ namespace boost { namespace numeric { namespace ublas {
                 restart (v, index, it, it_end);
                 // The proxies could reference the same container.
                 restart (e, index, ite, ite_end);
+#endif
             } else if (compare > 0) {
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
                 // Sparse proxies don't need to be conformant.
                 // Thanks to Michael Stevens for suggesting this.
                 size_type index (ite.index ());
@@ -477,8 +485,10 @@ namespace boost { namespace numeric { namespace ublas {
                 restart (e, index, ite, ite_end);
                 // The proxies could reference the same container.
                 restart (v, index, it, it_end);
+#endif
             }
         }
+#ifdef BOOST_UBLAS_NON_CONFORMANT_PROXIES
         while (ite != ite_end) {
             // Sparse proxies don't need to be conformant.
             // Thanks to Michael Stevens for suggesting this.
@@ -495,6 +505,7 @@ namespace boost { namespace numeric { namespace ublas {
             // The proxies could reference the same container.
             restart (v, index, it, it_end);
         }
+#endif
     }
 
     // Dispatcher
diff --git a/include/boost/numeric/ublas/vector_expression.hpp b/include/boost/numeric/ublas/vector_expression.hpp
index 75857b86..99f08afc 100644
--- a/include/boost/numeric/ublas/vector_expression.hpp
+++ b/include/boost/numeric/ublas/vector_expression.hpp
@@ -261,12 +261,12 @@ namespace boost { namespace numeric { namespace ublas {
         // Element access
         BOOST_UBLAS_INLINE
         const_reference operator () (size_type i) const {
-            return e_ (i);
+            return expression () (i);
         }
 
         BOOST_UBLAS_INLINE
         const_reference operator [] (size_type i) const {
-            return e_ [i];
+            return expression () [i];
         }
 
         typedef const_iterator_type const_iterator;
@@ -275,22 +275,22 @@ namespace boost { namespace numeric { namespace ublas {
         // Element lookup
         BOOST_UBLAS_INLINE
         const_iterator find_first (size_type i) const {
-            return e_.find_first (i);
+            return expression ().find_first (i);
         }
         BOOST_UBLAS_INLINE
         const_iterator find_last (size_type i) const {
-            return e_.find_last (i);
+            return expression ().find_last (i);
         }
 
         // Iterator is the iterator of the referenced expression.
 
         BOOST_UBLAS_INLINE
         const_iterator begin () const {
-            return e_.begin ();
+            return expression ().begin ();
         }
         BOOST_UBLAS_INLINE
         const_iterator end () const {
-            return e_.end ();
+            return expression ().end ();
         }
 
         // Reverse iterator
@@ -362,26 +362,26 @@ namespace boost { namespace numeric { namespace ublas {
         // Resizing
         BOOST_UBLAS_INLINE
         void resize (size_type size) {
-            e_.resize (size);
+            expression ().resize (size);
         }
 
         // Element access
         BOOST_UBLAS_INLINE
-        const_reference operator () (size_type i) const { 
-            return e_ (i); 
+        const_reference operator () (size_type i) const {
+            return expression () (i);
         }
         BOOST_UBLAS_INLINE
-        reference operator () (size_type i) { 
-            return e_ (i); 
+        reference operator () (size_type i) {
+            return expression () (i);
         }
 
         BOOST_UBLAS_INLINE
         const_reference operator [] (size_type i) const {
-            return e_ [i];
+            return expression () [i];
         }
         BOOST_UBLAS_INLINE
         reference operator [] (size_type i) {
-            return e_ [i];
+            return expression () [i];
         }
 
         typedef const_iterator_type const_iterator;
@@ -390,39 +390,39 @@ namespace boost { namespace numeric { namespace ublas {
         // Element lookup
         BOOST_UBLAS_INLINE
         const_iterator find_first (size_type i) const {
-            return e_.find_first (i);
+            return expression ().find_first (i);
         }
         BOOST_UBLAS_INLINE
         iterator find_first (size_type i) {
-            return e_.find_first (i);
+            return expression ().find_first (i);
         }
         BOOST_UBLAS_INLINE
         const_iterator find_last (size_type i) const {
-            return e_.find_last (i);
+            return expression ().find_last (i);
         }
         BOOST_UBLAS_INLINE
         iterator find_last (size_type i) {
-            return e_.find_last (i);
+            return expression ().find_last (i);
         }
 
         // Iterator is the iterator of the referenced expression.
 
         BOOST_UBLAS_INLINE
         const_iterator begin () const {
-            return e_.begin ();
+            return expression ().begin ();
         }
         BOOST_UBLAS_INLINE
         const_iterator end () const {
-            return e_.end (); 
+            return expression ().end ();
         }
 
         BOOST_UBLAS_INLINE
         iterator begin () {
-            return e_.begin ();
+            return expression ().begin ();
         }
         BOOST_UBLAS_INLINE
         iterator end () {
-            return e_.end (); 
+            return expression ().end ();
         }
 
         // Reverse iterator
@@ -481,7 +481,9 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E::const_closure_type expression_closure_type;
         typedef const vector_unary<E, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E::const_iterator const_iterator_type;
         typedef unknown_storage_tag storage_category;
 
@@ -499,7 +501,7 @@ namespace boost { namespace numeric { namespace ublas {
             return e_.size ();
         }
         BOOST_UBLAS_INLINE
-        const expression_type &expression () const {
+        const expression_closure_type &expression () const {
             return e_;
         }
 
@@ -664,17 +666,12 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression_type e_; 
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression_type &e_;
-#endif
+        expression_closure_type e_;
     };
 
     template<class E, class F>
     struct vector_unary_traits {
-        typedef vector_unary<typename E::const_closure_type, F> expression_type;
+        typedef vector_unary<E, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -701,7 +698,7 @@ namespace boost { namespace numeric { namespace ublas {
     }
 
     // (real v) [i] = real (v [i])
-    template<class E> 
+    template<class E>
     BOOST_UBLAS_INLINE
     typename vector_unary_traits<E, scalar_real<typename E::value_type> >::result_type
     real (const vector_expression<E> &e) {
@@ -710,7 +707,7 @@ namespace boost { namespace numeric { namespace ublas {
     }
 
     // (imag v) [i] = imag (v [i])
-    template<class E> 
+    template<class E>
     BOOST_UBLAS_INLINE
     typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::result_type
     imag (const vector_expression<E> &e) {
@@ -719,7 +716,7 @@ namespace boost { namespace numeric { namespace ublas {
     }
 
     // (trans v) [i] = v [i]
-    template<class E> 
+    template<class E>
     BOOST_UBLAS_INLINE
     typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::result_type
     trans (const vector_expression<E> &e) {
@@ -728,7 +725,7 @@ namespace boost { namespace numeric { namespace ublas {
     }
 
     // (herm v) [i] = conj (v [i])
-    template<class E> 
+    template<class E>
     BOOST_UBLAS_INLINE
     typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
     herm (const vector_expression<E> &e) {
@@ -753,7 +750,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const vector_binary<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E1::const_iterator const_iterator1_type;
         typedef typename E2::const_iterator const_iterator2_type;
         typedef unknown_storage_tag storage_category;
@@ -772,11 +772,11 @@ namespace boost { namespace numeric { namespace ublas {
             return BOOST_UBLAS_SAME (e1_.size (), e2_.size ()); 
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -1064,20 +1064,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_; 
-        expression2_type e2_; 
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_; 
-        const expression2_type &e2_; 
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class E1, class E2, class F>
     struct vector_binary_traits {
-        typedef vector_binary<typename E1::const_closure_type, 
-                              typename E2::const_closure_type, F> expression_type;
+        typedef vector_binary<E1, E2, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type; 
 #else
@@ -1126,7 +1119,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef E1 expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef const vector_binary_scalar1<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E1::value_type const_iterator1_type;
         typedef typename E2::const_iterator const_iterator2_type;
         typedef unknown_storage_tag storage_category;
@@ -1145,11 +1141,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e2_.size ();
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -1320,22 +1316,15 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_; 
-        expression2_type e2_; 
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_; 
-        const expression2_type &e2_; 
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class T1, class E2, class F>
     struct vector_binary_scalar1_traits {
-        typedef vector_binary_scalar1<scalar_const_reference<T1>, 
-                                      typename E2::const_closure_type, F> expression_type;
+        typedef vector_binary_scalar1<scalar_const_reference<T1>, E2, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
-        typedef expression_type result_type; 
+        typedef expression_type result_type;
 #else
         typedef vector<typename F::result_type> result_type;
 #endif
@@ -1345,7 +1334,7 @@ namespace boost { namespace numeric { namespace ublas {
     template<class T1, class E2>
     BOOST_UBLAS_INLINE
     typename vector_binary_scalar1_traits<T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
-    operator * (const T1 &e1, 
+    operator * (const T1 &e1,
                 const vector_expression<E2> &e2) {
         typedef BOOST_UBLAS_TYPENAME vector_binary_scalar1_traits<T1, E2, scalar_multiplies<T1, BOOST_UBLAS_TYPENAME E2::value_type> >::expression_type expression_type;
         return expression_type (e1, e2 ());
@@ -1368,7 +1357,10 @@ namespace boost { namespace numeric { namespace ublas {
         typedef const_reference reference;
         typedef const value_type *const_pointer;
         typedef const_pointer pointer;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef E2 expression2_closure_type;
         typedef const vector_binary_scalar2<E1, E2, F> const_closure_type;
+        typedef const_closure_type closure_type;
         typedef typename E1::const_iterator const_iterator1_type;
         typedef typename E2::value_type const_iterator2_type;
         typedef unknown_storage_tag storage_category;
@@ -1387,11 +1379,11 @@ namespace boost { namespace numeric { namespace ublas {
             return e1_.size (); 
         }
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -1562,20 +1554,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_; 
-        expression2_type e2_; 
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_; 
-        const expression2_type &e2_; 
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class E1, class T2, class F>
     struct vector_binary_scalar2_traits {
-        typedef vector_binary_scalar2<typename E1::const_closure_type,        
-                                      scalar_const_reference<T2>, F> expression_type;
+        typedef vector_binary_scalar2<E1, scalar_const_reference<T2>, F> expression_type;
 #ifdef BOOST_UBLAS_USE_ET
         typedef expression_type result_type;
 #else
@@ -1612,6 +1597,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename F::size_type size_type;
         typedef typename F::difference_type difference_type;
         typedef typename F::result_type value_type;
+        typedef typename E::const_closure_type expression_closure_type;
         typedef typename E::const_iterator::iterator_category iterator_category;
         typedef unknown_storage_tag storage_category;
 
@@ -1625,7 +1611,7 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Accessors
         BOOST_UBLAS_INLINE
-        const expression_type &expression () const {
+        const expression_closure_type &expression () const {
             return e_;
         }
 
@@ -1664,17 +1650,12 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression_type e_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression_type &e_;
-#endif
+        expression_closure_type e_;
     };
-    
+
     template<class E, class F>
     struct vector_scalar_unary_traits {
-        typedef vector_scalar_unary<typename E::const_closure_type, F> expression_type;
+        typedef vector_scalar_unary<E, F> expression_type;
 #if defined (BOOST_UBLAS_USE_ET) && defined (BOOST_UBLAS_USE_SCALAR_ET)
          typedef expression_type result_type; 
 #else
@@ -1740,6 +1721,8 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename F::size_type size_type;
         typedef typename F::difference_type difference_type;
         typedef typename F::result_type value_type;
+        typedef typename E1::const_closure_type expression1_closure_type;
+        typedef typename E2::const_closure_type expression2_closure_type;
         typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
                                                   typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
         typedef unknown_storage_tag storage_category;
@@ -1754,11 +1737,11 @@ namespace boost { namespace numeric { namespace ublas {
 
         // Accessors
         BOOST_UBLAS_INLINE
-        const expression1_type &expression1 () const {
+        const expression1_closure_type &expression1 () const {
             return e1_;
         }
         BOOST_UBLAS_INLINE
-        const expression2_type &expression2 () const {
+        const expression2_closure_type &expression2 () const {
             return e2_;
         }
 
@@ -1797,20 +1780,13 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-#ifdef BOOST_UBLAS_ET_VALUE
-        expression1_type e1_;
-        expression2_type e2_;
-#endif
-#ifdef BOOST_UBLAS_ET_REFERENCE
-        const expression1_type &e1_;
-        const expression2_type &e2_;
-#endif
+        expression1_closure_type e1_;
+        expression2_closure_type e2_;
     };
 
     template<class E1, class E2, class F>
     struct vector_scalar_binary_traits {
-        typedef vector_scalar_binary<typename E1::const_closure_type, 
-                                     typename E2::const_closure_type, F> expression_type;
+        typedef vector_scalar_binary<E1, E2, F> expression_type;
 #if defined (BOOST_UBLAS_USE_ET) && defined (BOOST_UBLAS_USE_SCALAR_ET)
         typedef expression_type result_type; 
 #else
diff --git a/include/boost/numeric/ublas/vector_proxy.hpp b/include/boost/numeric/ublas/vector_proxy.hpp
index 53f8dc60..be621940 100644
--- a/include/boost/numeric/ublas/vector_proxy.hpp
+++ b/include/boost/numeric/ublas/vector_proxy.hpp
@@ -37,7 +37,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename V::size_type size_type;
         typedef typename V::difference_type difference_type;
         typedef typename V::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename V::const_reference const_reference;
         typedef typename V::reference reference;
         typedef typename V::const_pointer const_pointer;
@@ -52,15 +52,16 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename V::const_pointer,
                                        typename V::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const vector_range<vector_type> > const_closure_type;
-        typedef vector_reference<vector_range<vector_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename V::closure_type vector_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<V>::value,
+                                       typename V::const_closure_type,
+                                       typename V::closure_type>::type vector_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const vector_range<vector_type> const_closure_type;
         typedef vector_range<vector_type> closure_type;
-#endif
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename V::const_iterator const_iterator_type;
         typedef typename V::iterator iterator_type;
 #else
@@ -84,6 +85,13 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (r_.start () <= data_.size () &&
                                r_.start () + r_.size () <= data_.size (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        vector_range (const vector_closure_type &data, const range<> &r, int):
+            data_ (data), r_ (r.preprocess (data.size ())) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (r_.start () <= data_.size () &&
+                               r_.start () + r_.size () <= data_.size (), bad_index ());
+        }
 #else
         BOOST_UBLAS_INLINE
         vector_range (vector_type &data, const range &r):
@@ -92,6 +100,13 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (r_.start () <= data_.size () &&
                                r_.start () + r_.size () <= data_.size (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        vector_range (const vector_closure_type &data, const range &r, int):
+            data_ (data), r_ (r) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (r_.start () <= data_.size () &&
+                               r_.start () + r_.size () <= data_.size (), bad_index ());
+        }
 #endif
 
         // Accessors
@@ -104,11 +119,11 @@ namespace boost { namespace numeric { namespace ublas {
             return r_.size ();
         }
         BOOST_UBLAS_INLINE
-        const_vector_type &data () const {
+        const vector_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        vector_type &data () {
+        vector_closure_type &data () {
             return data_;
         }
 
@@ -149,12 +164,12 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
         vector_range<vector_type> project (const range<> &r) const {
-            return vector_range<vector_type> (data_, r_.compose (r));
+            return vector_range<vector_type> (data_, r_.compose (r), 0);
         }
 #else
         BOOST_UBLAS_INLINE
         vector_range<vector_type> project (const range &r) const {
-            return vector_range<vector_type> (data_, r_.compose (r));
+            return vector_range<vector_type> (data_, r_.compose (r), 0);
         }
 #endif
 
@@ -298,8 +313,13 @@ namespace boost { namespace numeric { namespace ublas {
             public iterator_base_traits<typename const_iterator_type::iterator_category>::template
                         iterator_base<const_iterator, value_type>::type {
 #else
+#ifndef BOOST_MSVC
             public random_access_iterator_base<typename const_iterator_type::iterator_category,
                                                const_iterator, value_type> {
+#else
+            public random_access_iterator_base<typename V::const_iterator::iterator_category,
+                                               const_iterator, value_type> {
+#endif
 #endif
         public:
             typedef typename const_iterator_type::iterator_category iterator_category;
@@ -364,7 +384,7 @@ namespace boost { namespace numeric { namespace ublas {
                 return it_.index () - (*this) ().start ();
             }
 
-            // Assignment 
+            // Assignment
             BOOST_UBLAS_INLINE
             const_iterator &operator = (const const_iterator &it) {
                 container_const_reference<vector_range>::assign (&it ());
@@ -405,8 +425,13 @@ namespace boost { namespace numeric { namespace ublas {
             public iterator_base_traits<typename iterator_type::iterator_category>::template
                         iterator_base<iterator, value_type>::type {
 #else
+#ifndef BOOST_MSVC
             public random_access_iterator_base<typename iterator_type::iterator_category,
                                                iterator, value_type> {
+#else
+            public random_access_iterator_base<typename V::iterator::iterator_category,
+                                               iterator, value_type> {
+#endif
 #endif
         public:
             typedef typename iterator_type::iterator_category iterator_category;
@@ -531,7 +556,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        vector_type &data_;
+        vector_closure_type data_;
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         range<> r_;
 #else
@@ -553,7 +578,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class V>
     BOOST_UBLAS_INLINE
-    const vector_range<const V> project (const V &data, const range<> &r) {
+    vector_range<const V> project (const V &data, const range<> &r) {
         return vector_range<const V> (data, r);
     }
     template<class V>
@@ -571,7 +596,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class V>
     BOOST_UBLAS_INLINE
-    const vector_range<const V> project (const V &data, const range &r) {
+    vector_range<const V> project (const V &data, const range &r) {
         return vector_range<const V> (data, r);
     }
     template<class V>
@@ -595,7 +620,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename V::size_type size_type;
         typedef typename V::difference_type difference_type;
         typedef typename V::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename V::const_reference const_reference;
         typedef typename V::reference reference;
         typedef typename V::const_pointer const_pointer;
@@ -610,14 +635,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename V::const_pointer,
                                        typename V::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const vector_slice<vector_type> > const_closure_type;
-        typedef vector_reference<vector_slice<vector_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename V::closure_type vector_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<V>::value,
+                                       typename V::const_closure_type,
+                                       typename V::closure_type>::type vector_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const vector_slice<vector_type> const_closure_type;
         typedef vector_slice<vector_type> closure_type;
-#endif
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         typedef slice<>::const_iterator const_iterator_type;
         typedef slice<>::const_iterator iterator_type;
@@ -640,6 +666,13 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (s_.start () <= data_.size () &&
                                s_.start () + s_.stride () * (s_.size () - (s_.size () > 0)) <= data_.size (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        vector_slice (const vector_closure_type &data, const slice<> &s, int):
+            data_ (data), s_ (s.preprocess (data.size ())) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (s_.start () <= data_.size () &&
+                               s_.start () + s_.stride () * (s_.size () - (s_.size () > 0)) <= data_.size (), bad_index ());
+        }
 #else
         BOOST_UBLAS_INLINE
         vector_slice (vector_type &data, const slice &s):
@@ -648,6 +681,13 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_CHECK (s_.start () <= data_.size () &&
                                s_.start () + s_.stride () * (s_.size () - (s_.size () > 0)) <= data_.size (), bad_index ());
         }
+        BOOST_UBLAS_INLINE
+        vector_slice (const vector_closure_type &data, const slice &s, int):
+            data_ (data), s_ (s) {
+            // Early checking of preconditions.
+            BOOST_UBLAS_CHECK (s_.start () <= data_.size () &&
+                               s_.start () + s_.stride () * (s_.size () - (s_.size () > 0)) <= data_.size (), bad_index ());
+        }
 #endif
 
         // Accessors
@@ -664,11 +704,11 @@ namespace boost { namespace numeric { namespace ublas {
             return s_.size (); 
         }
         BOOST_UBLAS_INLINE
-        const_vector_type &data () const {
+        const vector_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        vector_type &data () {
+        vector_closure_type &data () {
             return data_;
         }
 
@@ -709,20 +749,20 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
         vector_slice<vector_type> project (const range<> &r) const {
-            return vector_slice<vector_type>  (data_, s_.compose (r));
+            return vector_slice<vector_type>  (data_, s_.compose (r), 0);
         }
         BOOST_UBLAS_INLINE
         vector_slice<vector_type> project (const slice<> &s) const {
-            return vector_slice<vector_type>  (data_, s_.compose (s));
+            return vector_slice<vector_type>  (data_, s_.compose (s), 0);
         }
 #else
         BOOST_UBLAS_INLINE
         vector_slice<vector_type> project (const range &r) const {
-            return vector_slice<vector_type>  (data_, s_.compose (r));
+            return vector_slice<vector_type>  (data_, s_.compose (r), 0);
         }
         BOOST_UBLAS_INLINE
         vector_slice<vector_type> project (const slice &s) const {
-            return vector_slice<vector_type>  (data_, s_.compose (s));
+            return vector_slice<vector_type>  (data_, s_.compose (s), 0);
         }
 #endif
 
@@ -1095,7 +1135,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        vector_type &data_;
+        vector_closure_type data_;
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         slice<> s_;
 #else
@@ -1124,7 +1164,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class V>
     BOOST_UBLAS_INLINE
-    const vector_slice<const V> project (const V &data, const slice<> &s) {
+    vector_slice<const V> project (const V &data, const slice<> &s) {
         return vector_slice<const V> (data, s);
     }
     template<class V>
@@ -1149,7 +1189,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class V>
     BOOST_UBLAS_INLINE
-    const vector_slice<const V> project (const V &data, const slice &s) {
+    vector_slice<const V> project (const V &data, const slice &s) {
         return vector_slice<const V> (data, s);
     }
     template<class V>
@@ -1177,7 +1217,7 @@ namespace boost { namespace numeric { namespace ublas {
         typedef typename V::size_type size_type;
         typedef typename V::difference_type difference_type;
         typedef typename V::value_type value_type;
-#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION
+#ifndef BOOST_UBLAS_CT_PROXY_BASE_TYPEDEFS
         typedef typename V::const_reference const_reference;
         typedef typename V::reference reference;
         typedef typename V::const_pointer const_pointer;
@@ -1192,14 +1232,15 @@ namespace boost { namespace numeric { namespace ublas {
                                        typename V::const_pointer,
                                        typename V::pointer>::type pointer;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_REFERENCE
-        typedef const vector_const_reference<const vector_indirect<vector_type, indirect_array_type> > const_closure_type;
-        typedef vector_reference<vector_indirect<vector_type, indirect_array_type> > closure_type;
+#ifndef BOOST_UBLAS_CT_PROXY_CLOSURE_TYPEDEFS
+        typedef typename V::closure_type vector_closure_type;
+#else
+        typedef typename detail::ct_if<boost::is_const<V>::value,
+                                       typename V::const_closure_type,
+                                       typename V::closure_type>::type vector_closure_type;
 #endif
-#ifdef BOOST_UBLAS_ET_CLOSURE_VALUE
         typedef const vector_indirect<vector_type, indirect_array_type> const_closure_type;
         typedef vector_indirect<vector_type, indirect_array_type> closure_type;
-#endif
         typedef typename IA::const_iterator const_iterator_type;
         typedef typename IA::const_iterator iterator_type;
         typedef typename storage_restrict_traits<typename V::storage_category,
@@ -1216,10 +1257,16 @@ namespace boost { namespace numeric { namespace ublas {
         BOOST_UBLAS_INLINE
         vector_indirect (vector_type &data, const indirect_array_type &ia):
             data_ (data), ia_ (ia.preprocess (data.size ())) {}
+        BOOST_UBLAS_INLINE
+        vector_indirect (const vector_closure_type &data, const indirect_array_type &ia, int):
+            data_ (data), ia_ (ia.preprocess (data.size ())) {}
 #else
         BOOST_UBLAS_INLINE
         vector_indirect (vector_type &data, const indirect_array_type &ia):
             data_ (data), ia_ (ia) {}
+        BOOST_UBLAS_INLINE
+        vector_indirect (const vector_closure_type &data, const indirect_array_type &ia, int):
+            data_ (data), ia_ (ia) {}
 #endif
 
         // Accessors
@@ -1228,11 +1275,11 @@ namespace boost { namespace numeric { namespace ublas {
             return ia_.size ();
         }
         BOOST_UBLAS_INLINE
-        const_vector_type &data () const {
+        const vector_closure_type &data () const {
             return data_;
         }
         BOOST_UBLAS_INLINE
-        vector_type &data () {
+        vector_closure_type &data () {
             return data_;
         }
         BOOST_UBLAS_INLINE
@@ -1281,25 +1328,25 @@ namespace boost { namespace numeric { namespace ublas {
 #ifdef BOOST_UBLAS_ENABLE_INDEX_SET_ALL
         BOOST_UBLAS_INLINE
         vector_indirect<vector_type, indirect_array_type> project (const range<> &r) const {
-            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (r));
+            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (r), 0);
         }
         BOOST_UBLAS_INLINE
         vector_indirect<vector_type, indirect_array_type> project (const slice<> &s) const {
-            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (s));
+            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (s), 0);
         }
 #else
         BOOST_UBLAS_INLINE
         vector_indirect<vector_type, indirect_array_type> project (const range &r) const {
-            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (r));
+            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (r), 0);
         }
         BOOST_UBLAS_INLINE
         vector_indirect<vector_type, indirect_array_type> project (const slice &s) const {
-            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (s));
+            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (s), 0);
         }
 #endif
         BOOST_UBLAS_INLINE
         vector_indirect<vector_type, indirect_array_type> project (const indirect_array_type &ia) const {
-            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (ia));
+            return vector_indirect<vector_type, indirect_array_type> (data_, ia_.compose (ia), 0);
         }
 
         // Assignment
@@ -1671,7 +1718,7 @@ namespace boost { namespace numeric { namespace ublas {
         }
 
     private:
-        vector_type &data_;
+        vector_closure_type data_;
         indirect_array_type ia_;
         static vector_type nil_;
     };
@@ -1716,7 +1763,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     // template<class V, class IA>
     // BOOST_UBLAS_INLINE
-    // const vector_indirect<const V, IA> project (const V &data, const IA &ia) {
+    // vector_indirect<const V, IA> project (const V &data, const IA &ia) {
     //     return vector_indirect<const V, IA> (data, ia);
     // }
     // template<class V, class IA>
@@ -1733,7 +1780,7 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
     template<class V, class A>
     BOOST_UBLAS_INLINE
-    const vector_indirect<const V, indirect_array<A> > project (const V &data, const indirect_array<A> &ia) {
+    vector_indirect<const V, indirect_array<A> > project (const V &data, const indirect_array<A> &ia) {
         return vector_indirect<const V, indirect_array<A> > (data, ia);
     }
     template<class V, class A>
diff --git a/include/boost/numeric/ublas/vector_sparse.hpp b/include/boost/numeric/ublas/vector_sparse.hpp
index cef33b67..970c4028 100644
--- a/include/boost/numeric/ublas/vector_sparse.hpp
+++ b/include/boost/numeric/ublas/vector_sparse.hpp
@@ -151,7 +151,7 @@ namespace boost { namespace numeric { namespace ublas {
         // FIXME: GCC 3.1 warns, if enabled
 #ifndef __GNUC__
         BOOST_UBLAS_INLINE
-        operator const const_reference () const {
+        operator const_reference () const {
             return d_;
         }
 #endif
@@ -184,7 +184,7 @@ namespace boost { namespace numeric { namespace ublas {
 #if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE)
         typedef T &reference;
 #elif defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE)
-        typedef sparse_vector_element<sparse_vector> reference;
+        typedef sparse_vector_element<sparse_vector<T, A> > reference;
 #elif defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE)
         typedef typename map_traits<A>::reference reference;
 #endif
@@ -577,7 +577,7 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_INLINE
             reference operator * () const {
                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
-#if ! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE)
+#if defined (BOOST_MSVC_STD_ITERATOR) || (! defined (BOOST_UBLAS_STRICT_STORAGE_SPARSE) && ! defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE))
                 return (*it_).second;
 #elif defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE)
                 return reference ((*this) (), &(*it_).second, index ());
@@ -666,12 +666,12 @@ namespace boost { namespace numeric { namespace ublas {
 
     // Array based sparse vector class
     // Thanks to Kresimir Fresl for extending this to cover different index bases.
-    template<class T, class IA, class TA, std::size_t IB>
+    template<class T, std::size_t IB, class IA, class TA>
     class compressed_vector:
-        public vector_expression<compressed_vector<T, IA, TA, IB> > {
+        public vector_expression<compressed_vector<T, IB, IA, TA> > {
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
-        BOOST_UBLAS_USING vector_expression<compressed_vector<T, IA, TA, IB> >::operator ();
+        BOOST_UBLAS_USING vector_expression<compressed_vector<T, IB, IA, TA> >::operator ();
 #endif
         typedef std::size_t size_type;
         typedef std::ptrdiff_t difference_type;
@@ -681,14 +681,14 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
         typedef T &reference;
 #else
-        typedef sparse_vector_element<compressed_vector> reference;
+        typedef sparse_vector_element<compressed_vector<T, IB, IA, TA> > reference;
 #endif
         typedef const T *const_pointer;
         typedef T *pointer;
         typedef IA index_array_type;
         typedef TA value_array_type;
-        typedef const compressed_vector<T, IA, TA, IB> const_self_type;
-        typedef compressed_vector<T, IA, TA, IB> self_type;
+        typedef const compressed_vector<T, IB, IA, TA> const_self_type;
+        typedef compressed_vector<T, IB, IA, TA> self_type;
         typedef const vector_const_reference<const_self_type> const_closure_type;
         typedef vector_reference<self_type> closure_type;
         typedef typename IA::const_iterator const_iterator_type;
@@ -733,6 +733,10 @@ namespace boost { namespace numeric { namespace ublas {
             return filled_;
         }
         BOOST_UBLAS_INLINE
+        static size_type index_base () {
+            return index_base_;
+        }
+        BOOST_UBLAS_INLINE
         const index_array_type &index_data () const {
             return index_data_;
         }
@@ -748,10 +752,6 @@ namespace boost { namespace numeric { namespace ublas {
         value_array_type &value_data () {
             return value_data_;
         }
-        BOOST_UBLAS_INLINE
-        size_type index_base () {
-            return index_base_;
-        }
 
         // Resizing
         BOOST_UBLAS_INLINE
@@ -1132,7 +1132,7 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_INLINE
             reference operator * () const {
                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
-#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
+#if defined (BOOST_MSVC_STD_ITERATOR) || ! defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE)
                 return (*this) ().value_data () [it_ - (*this) ().index_data ().begin ()];
 #else
                 return reference ((*this) (), &(*this) ().value_data () [it_ - (*this) ().index_data ().begin ()], index ());
@@ -1229,17 +1229,17 @@ namespace boost { namespace numeric { namespace ublas {
         friend class const_iterator;
     };
 
-    template<class T, class IA, class TA, std::size_t IB>
-    typename compressed_vector<T, IA, TA, IB>::value_type compressed_vector<T, IA, TA, IB>::zero_ = 0;
+    template<class T, std::size_t IB, class IA, class TA>
+    typename compressed_vector<T, IB, IA, TA>::value_type compressed_vector<T, IB, IA, TA>::zero_ = 0;
 
     // Array based sparse vector class
     // Thanks to Kresimir Fresl for extending this to cover different index bases.
-    template<class T, class IA, class TA, std::size_t IB>
+    template<class T, std::size_t IB, class IA, class TA>
     class coordinate_vector:
-        public vector_expression<coordinate_vector<T, IA, TA, IB> > {
+        public vector_expression<coordinate_vector<T, IB, IA, TA> > {
     public:
 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
-        BOOST_UBLAS_USING vector_expression<coordinate_vector<T, IA, TA, IB> >::operator ();
+        BOOST_UBLAS_USING vector_expression<coordinate_vector<T, IB, IA, TA> >::operator ();
 #endif
         typedef std::size_t size_type;
         typedef std::ptrdiff_t difference_type;
@@ -1249,14 +1249,14 @@ namespace boost { namespace numeric { namespace ublas {
 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
         typedef T &reference;
 #else
-        typedef sparse_vector_element<coordinate_vector> reference;
+        typedef sparse_vector_element<coordinate_vector<T, IB, IA, TA> > reference;
 #endif
         typedef const T *const_pointer;
         typedef T *pointer;
         typedef IA index_array_type;
         typedef TA value_array_type;
-        typedef const coordinate_vector<T, IA, TA, IB> const_self_type;
-        typedef coordinate_vector<T, IA, TA, IB> self_type;
+        typedef const coordinate_vector<T, IB, IA, TA> const_self_type;
+        typedef coordinate_vector<T, IB, IA, TA> self_type;
         typedef const vector_const_reference<const_self_type> const_closure_type;
         typedef vector_reference<self_type> closure_type;
         typedef typename IA::const_iterator const_iterator_type;
@@ -1301,6 +1301,10 @@ namespace boost { namespace numeric { namespace ublas {
             return filled_;
         }
         BOOST_UBLAS_INLINE
+        static size_type index_base () {
+            return index_base_;
+        }
+        BOOST_UBLAS_INLINE
         const index_array_type &index_data () const {
             return index_data_;
         }
@@ -1316,10 +1320,6 @@ namespace boost { namespace numeric { namespace ublas {
         value_array_type &value_data () {
             return value_data_;
         }
-        BOOST_UBLAS_INLINE
-        size_type index_base () {
-            return index_base_;
-        }
 
         // Resizing
         BOOST_UBLAS_INLINE
@@ -1505,8 +1505,8 @@ namespace boost { namespace numeric { namespace ublas {
         BOOST_UBLAS_INLINE
         void sort () const {
             if (! sorted_) {
-                typedef std::pair<typename index_array_type::value_type,
-                                  typename value_array_type::value_type> pair_type;
+                typedef std::pair<BOOST_UBLAS_TYPENAME index_array_type::value_type,
+                                  BOOST_UBLAS_TYPENAME value_array_type::value_type> pair_type;
                 unbounded_array<pair_type> a (filled_);
                 size_type i;
                 for (i = 0; i < filled_; ++ i)
@@ -1722,7 +1722,7 @@ namespace boost { namespace numeric { namespace ublas {
             BOOST_UBLAS_INLINE
             reference operator * () const {
                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
-#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
+#if defined (BOOST_MSVC_STD_ITERATOR) || ! defined (BOOST_UBLAS_STRICT_VECTOR_SPARSE)
                 return (*this) ().value_data () [it_ - (*this) ().index_data ().begin ()];
 #else
                 return reference ((*this) (), &(*this) ().value_data () [it_ - (*this) ().index_data ().begin ()], index ());
@@ -1820,8 +1820,8 @@ namespace boost { namespace numeric { namespace ublas {
         friend class const_iterator;
     };
 
-    template<class T, class IA, class TA, std::size_t IB>
-    typename coordinate_vector<T, IA, TA, IB>::value_type coordinate_vector<T, IA, TA, IB>::zero_ = 0;
+    template<class T, std::size_t IB, class IA, class TA>
+    typename coordinate_vector<T, IB, IA, TA>::value_type coordinate_vector<T, IB, IA, TA>::zero_ = 0;
 
 }}}
 
diff --git a/test1/test1.dsp b/test1/test1.dsp
index d56610b3..5745284a 100644
--- a/test1/test1.dsp
+++ b/test1/test1.dsp
@@ -150,10 +150,6 @@ SOURCE=..\..\..\..\boost\numeric\ublas\iterator.hpp
 # End Source File
 # Begin Source File
 
-SOURCE=..\..\..\..\boost\numeric\ublas\math.hpp
-# End Source File
-# Begin Source File
-
 SOURCE=..\..\..\..\boost\numeric\ublas\matrix.hpp
 # End Source File
 # Begin Source File
diff --git a/test2/test2.dsp b/test2/test2.dsp
index 2beeaac3..17f4b9bb 100644
--- a/test2/test2.dsp
+++ b/test2/test2.dsp
@@ -147,10 +147,6 @@ SOURCE=..\..\..\..\boost\numeric\ublas\iterator.hpp
 # End Source File
 # Begin Source File
 
-SOURCE=..\..\..\..\boost\numeric\ublas\math.hpp
-# End Source File
-# Begin Source File
-
 SOURCE=..\..\..\..\boost\numeric\ublas\matrix.hpp
 # End Source File
 # Begin Source File
diff --git a/test2/test2.hpp b/test2/test2.hpp
index b67cb9a2..6a98e39b 100644
--- a/test2/test2.hpp
+++ b/test2/test2.hpp
@@ -67,7 +67,9 @@ struct test_blas_3 {
 
 #define USE_FLOAT
 #define USE_DOUBLE
+#ifndef BOOST_MSVC
 #define USE_STD_COMPLEX
+#endif
 
 #endif
 
diff --git a/test3/Jamfile b/test3/Jamfile
index 6e43adef..d70b77a3 100644
--- a/test3/Jamfile
+++ b/test3/Jamfile
@@ -7,6 +7,6 @@ exe test3
     : <include>$(BOOST_ROOT)
       <define>BOOST_UBLAS_SAFE
       <define>USE_RANGE
-      <borland><*><cxxflags>"-w-8026 -w-8027 -w-8057 -w-8084 -w-8092"
+      <borland><*><cxxflags>"-w-8026 -w-8027 -w-8030 -w-8057 -w-8084 -w-8092"
     ;
 
diff --git a/test3/test3.dsp b/test3/test3.dsp
index be2cd916..3e3238fd 100644
--- a/test3/test3.dsp
+++ b/test3/test3.dsp
@@ -148,10 +148,6 @@ SOURCE=..\..\..\..\boost\numeric\ublas\iterator.hpp
 # End Source File
 # Begin Source File
 
-SOURCE=..\..\..\..\boost\numeric\ublas\math.hpp
-# End Source File
-# Begin Source File
-
 SOURCE=..\..\..\..\boost\numeric\ublas\matrix.hpp
 # End Source File
 # Begin Source File
diff --git a/test4/test4.dsp b/test4/test4.dsp
index 3c9da3f4..61e6c99a 100644
--- a/test4/test4.dsp
+++ b/test4/test4.dsp
@@ -148,10 +148,6 @@ SOURCE=..\..\..\..\boost\numeric\ublas\iterator.hpp
 # End Source File
 # Begin Source File
 
-SOURCE=..\..\..\..\boost\numeric\ublas\math.hpp
-# End Source File
-# Begin Source File
-
 SOURCE=..\..\..\..\boost\numeric\ublas\matrix.hpp
 # End Source File
 # Begin Source File
diff --git a/test5/test5.dsp b/test5/test5.dsp
index 730b744c..70f598f7 100644
--- a/test5/test5.dsp
+++ b/test5/test5.dsp
@@ -148,10 +148,6 @@ SOURCE=..\..\..\..\boost\numeric\ublas\iterator.hpp
 # End Source File
 # Begin Source File
 
-SOURCE=..\..\..\..\boost\numeric\ublas\math.hpp
-# End Source File
-# Begin Source File
-
 SOURCE=..\..\..\..\boost\numeric\ublas\matrix.hpp
 # End Source File
 # Begin Source File
diff --git a/test6/test6.dsp b/test6/test6.dsp
index 543a8992..ace806a0 100644
--- a/test6/test6.dsp
+++ b/test6/test6.dsp
@@ -148,10 +148,6 @@ SOURCE=..\..\..\..\boost\numeric\ublas\iterator.hpp
 # End Source File
 # Begin Source File
 
-SOURCE=..\..\..\..\boost\numeric\ublas\math.hpp
-# End Source File
-# Begin Source File
-
 SOURCE=..\..\..\..\boost\numeric\ublas\matrix.hpp
 # End Source File
 # Begin Source File