<em>guwi17_at_[hidden] wrote:<br><br></em>&gt;BTW. I would not expect a big difference between a banded matrix with 
<br>&gt;compile time known bandwidth compared to the current implementation. I 
<br>&gt;still wait for someone who proves the opposite ;-)<br>&gt;<br>&gt;Do you already have a benchmark to measure the speed of banded_matrix? 
<br><br>I am doing a very simple test with the (undocumented?) diagonal_matrix in banded.hpp<br><br>A matrix product of two diagonal_matrix should return another diagonal_matrix and should have linear complexity O(size). The same goes for an element-wise vector product. Storage should not be an issue, as the data layout of the diagonal_matrix should match that of the vector.
<br><br>The test performs one of both, and time measurements for a size of 1E7 double elements shows that the vector product is 10 times faster. The test also uses a matrix_vector_slice to get a view of the diagonal_matrix as a vector and does an element-wise product. This is much faster that the matrix product, but vector product is still twice as fast.
<br><br>IMO compile time knowledge of the bandwith of a banded_matrix *IS* important, specially for large sizes and small bandwidth (being a large diagonal_matrix the extreme case).<br><br>I suppose that implementing a banded_matrix_c with compile time bandwith is not trivial, at least not if optimal traversal of elements is required...
<br><br><br><br>test code:<br><br><br>#include &lt;boost/numeric/ublas/matrix.hpp&gt;<br>#include &lt;boost/numeric/ublas/matrix_proxy.hpp&gt;<br>#include &lt;boost/numeric/ublas/banded.hpp&gt;<br><br>#ifdef NDEBUG<br>&nbsp;&nbsp;&nbsp; const int SIZE = 10000000;
<br>#else<br>&nbsp;&nbsp;&nbsp; const int SIZE = 10;<br>#endif<br><br>void diagonal()<br>{<br>&nbsp;&nbsp;&nbsp; using namespace boost::numeric::ublas;<br>&nbsp;&nbsp;&nbsp; diagonal_matrix&lt;double&gt; d1 ( SIZE );<br>&nbsp;&nbsp;&nbsp; diagonal_matrix&lt;double&gt; d2 ( SIZE );
<br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; // how to do this with diagonal_adaptor&lt;scalar_matrix&gt; ?<br>&nbsp;&nbsp;&nbsp; std::fill(d1.data().begin(), d1.data().end(), 2.0);<br>&nbsp;&nbsp;&nbsp; std::fill(d2.data().begin(), d2.data().end(), 10.0);<br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; diagonal_matrix&lt;double&gt; d3 = prod(d1,d2);
<br>}<br><br>void diagonal_as_vector()<br>{<br>&nbsp;&nbsp;&nbsp; using namespace boost::numeric::ublas;<br>&nbsp;&nbsp;&nbsp; diagonal_matrix&lt;double&gt; d1 ( SIZE );<br>&nbsp;&nbsp;&nbsp; diagonal_matrix&lt;double&gt; d2 ( SIZE );<br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; // how to do this with diagonal_adaptor&lt;scalar_matrix&gt; ?
<br>&nbsp;&nbsp;&nbsp; std::fill(d1.data().begin(), d1.data().end(), 2.0);<br>&nbsp;&nbsp;&nbsp; std::fill(d2.data().begin(), d2.data().end(), 10.0);<br><br>&nbsp;&nbsp;&nbsp; matrix_vector_slice&lt; diagonal_matrix&lt;double&gt; &gt; dx1(d1, slice(0,1,SIZE), slice(0,1,SIZE));
<br>&nbsp;&nbsp;&nbsp; matrix_vector_slice&lt; diagonal_matrix&lt;double&gt; &gt; dx2(d2, slice(0,1,SIZE), slice(0,1,SIZE));<br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; // any vector_expression to diagonal_matrix ?<br>&nbsp;&nbsp;&nbsp; vector&lt;double&gt; d3 = element_prod(dx1, dx2);
<br>}<br><br>void vectors()<br>{<br>&nbsp;&nbsp;&nbsp; using namespace boost::numeric::ublas;<br>&nbsp;&nbsp;&nbsp; vector&lt;double&gt; v1(scalar_vector&lt;double&gt;(SIZE, 2.0));<br>&nbsp;&nbsp;&nbsp; vector&lt;double&gt; v2(scalar_vector&lt;double&gt;(SIZE, 10.0));
<br>&nbsp;&nbsp;&nbsp; vector&lt;double&gt; v3 = element_prod(v1, v2);<br>}<br><br>int main(int argc, char* argv[])<br>{<br>&nbsp;&nbsp;&nbsp; if (argc != 2) return 1;<br>&nbsp;&nbsp;&nbsp; if (*argv[1] == &#39;d&#39;) diagonal();<br>&nbsp;&nbsp;&nbsp; else if (*argv[1] == &#39;v&#39;) vectors();
<br>&nbsp;&nbsp;&nbsp; else if (*argv[1] == &#39;x&#39;) diagonal_as_vector();<br>&nbsp;&nbsp;&nbsp; return 0;<br>}<br><br>