Hi,
I'm currently trying to implement matrix and vector proxies in
PyViennaCL, and I can't get my matrices to look right. Suppose I have
the following arbitrary 5x5 matrix, as displayed in Python:
>>> m.value
array([[ 1., 2., 3., 4., 0.],
[ 5., 6., 7., 0., 8.],
[ 9., 10., 0., 11., 12.],
[ 13., 0., 14., 15., 16.],
[ 0., 17., 18., 19., 20.]])
And I want to extract the following arbitrary slices (indices from 0):
>>> m.value[1:4:2,1:4:2]
array([[ 6., 0.],
[ 0., 15.]])
This works fine, because I'm using NumPy's array slicing. But suppose I
use the ViennaCL slice/proxy interface:
>>> m[1:4:2,1:4:2].value
start: 1, 1
stride: 2, 2
size: 2, 2
array([[ 1., 2.],
[ 3., 4.]])
(The 'start', 'stride', and 'size' lines are debug output I put in
matrix_proxy.hpp). This clearly isn't right, and yet my slice objects
seem to be constructed correctly! If I shift my start/stop indices one
place to the left, then I get the same incorrect matrix, even though the
slice objects are still constructed as expected:
>>> m[0:3:2,0:3:2].value
start: 0, 0
stride: 2, 2
size: 2, 2
array([[ 1., 2.],
[ 3., 4.]])
I dug around in the matrix_proxy.hpp and vector_proxy.hpp files, and I
came across a couple of apparent inconsistencies. I patched those (see
below), but I don't know if my patch is right or wrong, since I'm not
creating proxies of proxies, and thus the behaviour hasn't changed. I
don't really know where to dig next to unearth the root of this
behaviour: whether it's a bug in my code, or in ViennaCL..
Oh, and the vector code seems to work fine, even though it's called in
almost exactly the same way as the matrix code!
Can you help? :)
T
diff --git a/viennacl/matrix_proxy.hpp b/viennacl/matrix_proxy.hpp
index d008888..ea31ffc 100644
--- a/viennacl/matrix_proxy.hpp
+++ b/viennacl/matrix_proxy.hpp
@@ -316,8 +316,8 @@ namespace viennacl
matrix_slice(MatrixType & A,
slice const & row_slice,
slice const & col_slice) : base_type(A.handle(),
- row_slice.size(), row_slice.start(), row_slice.stride(), A.internal_size1(),
- col_slice.size(), col_slice.start(), col_slice.stride(), A.internal_size2()) {}
+ row_slice.size(), A.start1() + A.stride1() * row_slice.start(), A.stride1() * row_slice.stride(), A.internal_size1(),
+ col_slice.size(), A.start2() + A.stride2() * col_slice.start(), A.stride2() * col_slice.stride(), A.internal_size2()) {}
using base_type::operator=;
diff --git a/viennacl/vector.hpp b/viennacl/vector.hpp
index c762a09..b18d5c0 100644
--- a/viennacl/vector.hpp
+++ b/viennacl/vector.hpp
@@ -2166,8 +2166,6 @@ namespace viennacl
}
};
-
-
// generic x = vec_expr1 + vec_expr2:
template <typename T, typename LHS, typename RHS>
struct op_executor<vector_base<T>, op_assign, vector_expression<const LHS, const RHS, op_add> >
diff --git a/viennacl/vector_proxy.hpp b/viennacl/vector_proxy.hpp
index 6d9acbf..57e2829 100644
--- a/viennacl/vector_proxy.hpp
+++ b/viennacl/vector_proxy.hpp
@@ -140,7 +140,7 @@ namespace viennacl
vector_range<VectorType> project(viennacl::vector_range<VectorType> & vec, viennacl::range const & r1)
{
assert(r1.size() <= vec.size() && bool("Size of range invalid!"));
- return vector_range<VectorType>(vec.get(), viennacl::range(vec.start() + r1.start(), vec.start() + r1.start() + r1.size()));
+ return vector_range<VectorType>(vec, viennacl::range(vec.start() + r1.start(), vec.start() + r1.start() + r1.size()));
}
//
@@ -244,7 +244,7 @@ namespace viennacl
vector_slice<VectorType> project(viennacl::vector_slice<VectorType> & vec, viennacl::slice const & s1)
{
assert(s1.size() <= vec.size() && bool("Size of slice larger than vector proxy!"));
- return vector_slice<VectorType>(vec.get(), viennacl::slice(vec.start() + s1.start(), vec.stride() * s1.stride(), s1.size()));
+ return vector_slice<VectorType>(vec, viennacl::slice(vec.start() + s1.start(), vec.stride() * s1.stride(), s1.size()));
}
// interaction with range and vector_range:
@@ -253,14 +253,14 @@ namespace viennacl
vector_slice<VectorType> project(viennacl::vector_slice<VectorType> & vec, viennacl::range const & r1)
{
assert(r1.size() <= vec.size() && bool("Size of slice larger than vector proxy!"));
- return vector_slice<VectorType>(vec.get(), viennacl::slice(vec.start() + r1.start(), vec.stride(), r1.size()));
+ return vector_slice<VectorType>(vec, viennacl::slice(vec.start() + r1.start(), vec.stride(), r1.size()));
}
template <typename VectorType>
vector_slice<VectorType> project(viennacl::vector_range<VectorType> & vec, viennacl::slice const & s1)
{
assert(s1.size() <= vec.size() && bool("Size of slice larger than vector proxy!"));
- return vector_slice<VectorType>(vec.get(), viennacl::range(vec.start() + s1.start(), s1.stride(), s1.size()));
+ return vector_slice<VectorType>(vec, viennacl::range(vec.start() + s1.start(), s1.stride(), s1.size()));
}
------------------------------------------------------------------------------
Get 100% visibility into Java/.NET code with AppDynamics Lite!
It's a free troubleshooting tool designed for production.
Get down to code-level detail for bottlenecks, with <2% overhead.
Download for free and get started troubleshooting in minutes.
http://pubads.g.doubleclick.net/gampad/clk?id=48897031&iu=/4140/ostg.clktrk
_______________________________________________
ViennaCL-devel mailing list
ViennaCL-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/viennacl-devel
