[selgen]

[Stefan Seefeld <stefan@xxxxxxxxxxxxxxxx>]

The attached patch implements all functions from section 9.4 ([selgen])
of the spec, i.e.

* first()
* indexbool()
* gather()
* scatter()
* clip()
* invclip()
* swap()

together with unit tests. It also contains some bits and pieces I
submitted earlier today to cleanup header dependencies etc., as I
wasn't able to easily separate the two.

Regards,
		Stefan

Index: src/vsip/dense.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/dense.hpp,v
retrieving revision 1.28
diff -u -r1.28 dense.hpp
--- src/vsip/dense.hpp	16 Sep 2005 21:51:08 -0000	1.28
+++ src/vsip/dense.hpp	27 Sep 2005 20:36:01 -0000
@@ -23,7 +23,6 @@
 #include <vsip/impl/layout.hpp>
 #include <vsip/impl/extdata.hpp>
 #include <vsip/impl/block-traits.hpp>
-#include <vsip/impl/view_traits.hpp>
 #include <vsip/impl/point.hpp>
 
 /// Complex storage format for dense blocks.
Index: src/vsip/matrix.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/matrix.hpp,v
retrieving revision 1.26
diff -u -r1.26 matrix.hpp
--- src/vsip/matrix.hpp	16 Sep 2005 22:03:20 -0000	1.26
+++ src/vsip/matrix.hpp	27 Sep 2005 20:36:01 -0000
@@ -40,11 +40,9 @@
 {
 
 /// View which appears as a two-dimensional, read-only matrix.
-template <typename T = VSIP_DEFAULT_VALUE_TYPE,
-	  typename Block = Dense<2, T> >
-class const_Matrix :
-  public impl::Non_assignable,
-  public  vsip::impl_const_View<vsip::const_Matrix,Block>
+template <typename T, typename Block>
+class const_Matrix : public impl::Non_assignable,
+		     public  vsip::impl_const_View<vsip::const_Matrix,Block>
 {
   typedef vsip::impl_const_View<vsip::const_Matrix,Block> impl_base_type;
   typedef typename impl::Lvalue_factory_type<Block>::type impl_factory_type;
@@ -156,10 +154,8 @@
 /// inherits from const_Matrix, so only the members that const_Matrix
 /// does not carry, or that are different, need be specified.
 
-template <typename T = VSIP_DEFAULT_VALUE_TYPE,
-	  typename Block = Dense<2, T> >
-class Matrix : 
-  public vsip::impl_View<vsip::Matrix,Block>
+template <typename T, typename Block>
+class Matrix : public vsip::impl_View<vsip::Matrix,Block>
 {
   typedef vsip::impl_View<vsip::Matrix,Block> impl_base_type;
   typedef typename impl::Lvalue_factory_type<Block>::type impl_factory_type;
@@ -374,6 +370,20 @@
   static bool const value = true;
 };
 
+template <typename T, typename Block>
+T
+get(const_Matrix<T, Block> view, Index<2> const &i)
+{
+  return view.get(i[0], i[1]);
+}
+
+template <typename T, typename Block>
+void
+put(Matrix<T, Block> view, Index<2> const &i, T value)
+{
+  view.put(i[0], i[1], value);
+}
+
 } // namespace vsip::impl
 
 } // namespace vsip
Index: src/vsip/selgen.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/selgen.hpp,v
retrieving revision 1.1
diff -u -r1.1 selgen.hpp
--- src/vsip/selgen.hpp	26 Sep 2005 20:23:29 -0000	1.1
+++ src/vsip/selgen.hpp	27 Sep 2005 20:36:01 -0000
@@ -1,7 +1,7 @@
 /* Copyright (c) 2005 by CodeSourcery, LLC.  All rights reserved. */
 
 /** @file    vsip/selgen.hpp
-    @author  Jules Bergmann
+    @author  Jules Bergmann, Stefan Seefeld
     @date    2005-08-15
     @brief   VSIPL++ Library: Selection functions.
 
@@ -15,9 +15,10 @@
 ***********************************************************************/
 
 #include <vsip/impl/block-traits.hpp>
+#include <vsip/impl/expr_unary_block.hpp>
 #include <vsip/impl/expr_generator_block.hpp>
 #include <vsip/vector.hpp>
-
+#include <vsip/matrix.hpp>
 
 
 /***********************************************************************
@@ -30,6 +31,29 @@
 namespace impl
 {
 
+template <typename T, typename B1, typename B2>
+length_type
+indexbool(const_Vector<T, B1> source, Vector<Index<1>, B2> indices)
+{
+  index_type cursor = 0;
+  for (index_type i = 0; i != source.size(); ++i)
+    if (source.get(i))
+      indices.put(cursor++, Index<1>(i));
+  return cursor;
+}
+
+template <typename T, typename B1, typename B2>
+length_type
+indexbool(const_Matrix<T, B1> source, Vector<Index<2>, B2> indices)
+{
+  index_type cursor = 0;
+  for (index_type r = 0; r != source.size(0); ++r)
+    for (index_type c = 0; c != source.size(1); ++c)
+      if (source.get(r, c))
+	indices.put(cursor++, Index<2>(r, c));
+  return cursor;
+}
+
 /// Generator functor for ramp.
 
 template <typename T>
@@ -59,7 +83,56 @@
 
 } // namespace vsip::impl
 
+template <typename Predicate,
+	  typename T1, typename T2,
+	  typename B1, typename B2>
+index_type
+first(index_type begin, Predicate p, 
+      const_Vector<T1, B1> v, const_Vector<T2, B2> w) VSIP_NOTHROW
+{
+  assert(impl::view_domain(v).element_conformant(impl::view_domain(w)));
+  if (begin > v.size()) 
+    return begin;
+  for (index_type i = begin, end = v.size(); i != end; ++i)
+    if (p(v.get(i), w.get(i)))
+      return i;
+  return v.size();
+}
+
+template <template <typename, typename> class const_View,
+	  typename T, typename B1, typename B2>
+length_type
+indexbool(const_View<T, B1> source, 
+	  Vector<Index<const_View<T, B2>::dim>, B2> indices)
+VSIP_NOTHROW
+{
+  return impl::indexbool(source, indices);
+}
+
+template <template <typename, typename> class const_View,
+	  typename T, typename B1, typename B2>
+Vector<T, Dense<1, T> >
+gather(const_View<T, B1> source,
+       const_Vector<Index<const_View<T, B2>::dim>, B2> indices)
+VSIP_NOTHROW
+{
+  Vector<T, Dense<1, T> > result(indices.size());
+  for (index_type i = 0; i != indices.size(); ++i)
+    result.put(i, get(source, indices.get(i)));
+  return result;
+}
 
+template <template <typename, typename> class View,
+	  typename T, typename B1, typename B2, typename B3>
+void
+scatter(const_Vector<T, B1> source,
+	const_Vector<Index<View<T, B3>::dim>, B2>  indices,
+	View<T, B3> destination)
+VSIP_NOTHROW
+{
+  for (index_type i = 0; i != indices.size(); ++i)
+    put(destination, indices.get(i), source.get(i));
+}
 
 /// Generate a linear ramp.
 
@@ -85,6 +158,126 @@
   return const_Vector<T, block_type>(block);
 }
 
+namespace impl
+{
+template <typename Tout, typename Tin1>
+struct clip_wrapper
+{
+  template <typename Tin0>
+  struct clip_functor
+  {
+    typedef Tout result_type;
+    result_type operator()(Tin0 t) const 
+    {
+      return t <= lower_threshold ? lower_clip_value 
+	: t < upper_threshold ? t
+	: upper_clip_value;
+    }
+
+    Tin1 lower_threshold;
+    Tin1 upper_threshold;
+    result_type lower_clip_value;
+    result_type upper_clip_value;
+  };
+  template <typename Tin0>
+  struct invclip_functor
+  {
+    typedef Tout result_type;
+    result_type operator()(Tin0 t) const 
+    {
+      return t < lower_threshold ? t
+	: t < middle_threshold ? lower_clip_value
+	: t <= upper_threshold ? upper_clip_value
+	: t;
+    }
+
+    Tin1 lower_threshold;
+    Tin1 middle_threshold;
+    Tin1 upper_threshold;
+    result_type lower_clip_value;
+    result_type upper_clip_value;
+  };
+};
+  
+}
+
+template <typename Tout, typename Tin0, typename Tin1,
+	  template <typename, typename> class const_View,
+	  typename Block>
+const_View<Tout,
+  impl::Unary_expr_block<const_View<Tin0, Block>::dim,
+			 impl::clip_wrapper<Tout, Tin1>::template clip_functor,
+			 Block, Tin0> const>
+clip(const_View<Tin0, Block> v, Tin1 lower_threshold, Tin1 upper_threshold,
+     Tout lower_clip_value, Tout upper_clip_value)
+{
+  typedef impl::Unary_expr_block<const_View<Tin0, Block>::dim,
+    impl::clip_wrapper<Tout, Tin1>::template clip_functor,
+    Block, Tin0> block_type;
+
+  typename impl::clip_wrapper<Tout, Tin1>::template clip_functor<Tin0> functor;
+  functor.lower_threshold = lower_threshold;
+  functor.upper_threshold = upper_threshold;
+  functor.lower_clip_value = lower_clip_value;
+  functor.upper_clip_value = upper_clip_value;
+
+  return const_View<Tout, block_type const>(block_type(v.block(), functor));
+}
+
+template <typename Tout, typename Tin0, typename Tin1,
+	  template <typename, typename> class const_View,
+	  typename Block>
+const_View<Tout,
+  impl::Unary_expr_block<const_View<Tin0, Block>::dim,
+			 impl::clip_wrapper<Tout, Tin1>::template invclip_functor,
+			 Block, Tin0> const>
+invclip(const_View<Tin0, Block> v,
+	Tin1 lower_threshold, Tin1 middle_threshold, Tin1 upper_threshold,
+	Tout lower_clip_value, Tout upper_clip_value)
+{
+  typedef impl::Unary_expr_block<const_View<Tin0, Block>::dim,
+    impl::clip_wrapper<Tout, Tin1>::template invclip_functor,
+    Block, Tin0> block_type;
+
+  typename impl::clip_wrapper<Tout, Tin1>::template invclip_functor<Tin0> functor;
+  functor.lower_threshold = lower_threshold;
+  functor.middle_threshold = middle_threshold;
+  functor.upper_threshold = upper_threshold;
+  functor.lower_clip_value = lower_clip_value;
+  functor.upper_clip_value = upper_clip_value;
+
+  return const_View<Tout, block_type const>(block_type(v.block(), functor));
+}
+
+namespace impl
+{
+/// Generic swapping of the content of two blocks.
+template <typename Block1, typename Block2>
+struct Swap
+{
+  static void apply(Block1 &block1, Block2 &block2)
+  {
+    assert(block1.size() == block2.size());
+    for (index_type i = 0; i != block1.size(); ++i)
+    {
+      typename Block1::value_type tmp = block1.get(i);
+      block1.put(i, block2.get(i));
+      block2.put(i, tmp);
+    }
+
+  }
+};
+}
+
+template <typename T1, typename T2,
+	  template <typename, typename> class View,
+	  typename Block1, typename Block2>
+inline void
+swap(View<T1, Block1> v, View<T2, Block2> w)
+{
+  impl::Swap<Block1, Block2>::apply(v.block(), w.block());
+}
+
 } // namespace vsip
 
 #endif // VSIP_SELGEN_HPP
Index: src/vsip/tensor.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/tensor.hpp,v
retrieving revision 1.19
diff -u -r1.19 tensor.hpp
--- src/vsip/tensor.hpp	16 Sep 2005 22:03:20 -0000	1.19
+++ src/vsip/tensor.hpp	27 Sep 2005 20:36:01 -0000
@@ -42,11 +42,9 @@
 enum whole_domain_type { whole_domain};
 
 /// View which appears as a three-dimensional, read-only tensor.
-template <typename T = VSIP_DEFAULT_VALUE_TYPE,
-	  typename Block = Dense<3, T> >
-class const_Tensor :
-  public impl::Non_assignable,
-  public  vsip::impl_const_View<vsip::const_Tensor,Block>
+template <typename T, typename Block>
+class const_Tensor : public impl::Non_assignable,
+		     public  vsip::impl_const_View<vsip::const_Tensor,Block>
 {
   typedef vsip::impl_const_View<vsip::const_Tensor,Block> impl_base_type;
   typedef typename impl::Lvalue_factory_type<Block>::type impl_factory_type;
@@ -263,8 +261,7 @@
 /// inherits from const_Tensor, so only the members that const_Tensor
 /// does not carry, or that are different, need be specified.
 
-template <typename T = VSIP_DEFAULT_VALUE_TYPE,
-	  typename Block = Dense<3, T> >
+template <typename T, typename Block>
 class Tensor : public vsip::impl_View<vsip::Tensor, Block>
 {
   typedef vsip::impl_View<vsip::Tensor, Block> impl_base_type;
@@ -604,6 +601,20 @@
   static bool const value = true;
 };
 
+template <typename T, typename Block>
+T
+get(const_Tensor<T, Block> view, Index<3> const &i)
+{
+  return view.get(i[0], i[1], i[2]);
+}
+
+template <typename T, typename Block>
+void
+put(Tensor<T, Block> view, Index<3> const &i, T value)
+{
+  view.put(i[0], i[1], i[2], value);
+}
+
 } // namespace vsip::impl
 
 } // namespace vsip
Index: src/vsip/vector.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/vector.hpp,v
retrieving revision 1.31
diff -u -r1.31 vector.hpp
--- src/vsip/vector.hpp	26 Sep 2005 20:11:05 -0000	1.31
+++ src/vsip/vector.hpp	27 Sep 2005 20:36:02 -0000
@@ -39,11 +39,9 @@
 {
 
 /// View which appears as a one-dimensional, read-only vector.
-template <typename T = VSIP_DEFAULT_VALUE_TYPE,
-	  typename Block = Dense<1, T> >
-class const_Vector :
-	public impl::Non_assignable,
-	public vsip::impl_const_View<vsip::const_Vector,Block>
+template <typename T, typename Block>
+class const_Vector : public impl::Non_assignable,
+		     public vsip::impl_const_View<vsip::const_Vector,Block>
 {
   typedef vsip::impl_const_View<vsip::const_Vector,Block> impl_base_type;
   typedef typename impl::Lvalue_factory_type<Block>::type impl_factory_type;
@@ -156,8 +154,7 @@
 /// inherits from const_Vector, so only the members that const_Vector
 /// does not carry, or that are different, need be specified.
 
-template <typename T = VSIP_DEFAULT_VALUE_TYPE,
-	  typename Block = Dense<1, T> >
+template <typename T, typename Block>
 class Vector : public vsip::impl_View<vsip::Vector,Block>
 {
   typedef vsip::impl_View<vsip::Vector,Block> impl_base_type;
@@ -321,7 +318,19 @@
   static bool const value = true;
 };
 
+template <typename T, typename Block>
+T
+get(const_Vector<T, Block> view, Index<1> const &i)
+{
+  return view.get(i[0]);
+}
 
+template <typename T, typename Block>
+void
+put(Vector<T, Block> view, Index<1> const &i, T value)
+{
+  view.put(i[0], value);
+}
 
 // Return the view extent as a domain.
 
Index: src/vsip/impl/expr_functor.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/impl/expr_functor.hpp,v
retrieving revision 1.6
diff -u -r1.6 expr_functor.hpp
--- src/vsip/impl/expr_functor.hpp	15 Sep 2005 14:49:25 -0000	1.6
+++ src/vsip/impl/expr_functor.hpp	27 Sep 2005 20:36:02 -0000
@@ -18,6 +18,7 @@
 #include <vsip/impl/expr_unary_block.hpp>
 #include <vsip/impl/expr_binary_block.hpp>
 #include <vsip/impl/expr_ternary_block.hpp>
+#include <vsip/impl/expr_binary_operators.hpp>
 
 namespace vsip
 {
Index: src/vsip/impl/matvec.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/impl/matvec.hpp,v
retrieving revision 1.1
diff -u -r1.1 matvec.hpp
--- src/vsip/impl/matvec.hpp	19 Sep 2005 21:06:46 -0000	1.1
+++ src/vsip/impl/matvec.hpp	27 Sep 2005 20:36:02 -0000
@@ -17,6 +17,8 @@
 
 #include <vsip/vector.hpp>
 #include <vsip/matrix.hpp>
+#include <vsip/impl/promote.hpp>
+#include <vsip/impl/fns_elementwise.hpp>
 
 
 namespace vsip
Index: src/vsip/impl/view_traits.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/impl/view_traits.hpp,v
retrieving revision 1.9
diff -u -r1.9 view_traits.hpp
--- src/vsip/impl/view_traits.hpp	28 Aug 2005 00:22:39 -0000	1.9
+++ src/vsip/impl/view_traits.hpp	27 Sep 2005 20:36:02 -0000
@@ -15,6 +15,7 @@
 ***********************************************************************/
 
 #include <vsip/support.hpp>
+#include <vsip/dense.hpp>
 #include <vsip/impl/subblock.hpp>
 #include <complex>
 
@@ -51,12 +52,18 @@
 
 } // impl
 
-template <typename T, typename B> struct Vector;
-template <typename T, typename B> struct Matrix;
-template <typename T, typename B> struct Tensor;
-template <typename T, typename B> struct const_Vector;
-template <typename T, typename B> struct const_Matrix;
-template <typename T, typename B> struct const_Tensor;
+template <typename T = VSIP_DEFAULT_VALUE_TYPE,
+	  typename B = Dense<1, T> > struct Vector;
+template <typename T = VSIP_DEFAULT_VALUE_TYPE,
+	  typename B = Dense<2, T> > struct Matrix;
+template <typename T = VSIP_DEFAULT_VALUE_TYPE,
+	  typename B = Dense<3, T> > struct Tensor;
+template <typename T = VSIP_DEFAULT_VALUE_TYPE,
+	  typename B = Dense<1, T> > struct const_Vector;
+template <typename T = VSIP_DEFAULT_VALUE_TYPE,
+	  typename B = Dense<2, T> > struct const_Matrix;
+template <typename T = VSIP_DEFAULT_VALUE_TYPE,
+	  typename B = Dense<3, T> > struct const_Tensor;
 
 namespace impl
 {
Index: tests/test.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/test.hpp,v
retrieving revision 1.6
diff -u -r1.6 test.hpp
--- tests/test.hpp	23 Aug 2005 19:32:53 -0000	1.6
+++ tests/test.hpp	27 Sep 2005 20:36:02 -0000
@@ -16,8 +16,11 @@
 
 #include <cstdlib>
 
+#include <vsip/support.hpp>
 #include <vsip/complex.hpp>
-
+#include <vsip/vector.hpp>
+#include <vsip/matrix.hpp>
+#include <vsip/tensor.hpp>
 
 
 /***********************************************************************
@@ -119,7 +122,43 @@
          equal(val1.imag(), val2.imag());
 }
 
+template <typename T, typename Block1, typename Block2>
+inline bool
+view_equal(vsip::const_Vector<T, Block1> v, vsip::const_Vector<T, Block2> w)
+{
+  if (v.size() != w.size()) return false;
+  for (vsip::length_type i = 0; i != v.size(); ++i)
+    if (!equal(v.get(i), w.get(i)))
+      return false;
+  return true;
+}
 
+template <typename T, typename Block1, typename Block2>
+inline bool
+view_equal(vsip::const_Matrix<T, Block1> v, vsip::const_Matrix<T, Block2> w)
+{
+  if (v.size(0) != w.size(0) || v.size(1) != w.size(1)) return false;
+  for (vsip::length_type i = 0; i != v.size(0); ++i)
+    for (vsip::length_type j = 0; j != v.size(1); ++j)
+      if (!equal(v.get(i, j), w.get(i, j)))
+	return false;
+  return true;
+}
+
+template <typename T, typename Block1, typename Block2>
+inline bool
+view_equal(vsip::const_Tensor<T, Block1> v, vsip::const_Tensor<T, Block2> w)
+{
+  if (v.size(0) != w.size(0) ||
+      v.size(1) != w.size(1) ||
+      v.size(2) != w.size(2)) return false;
+  for (vsip::length_type i = 0; i != v.size(0); ++i)
+    for (vsip::length_type j = 0; j != v.size(1); ++j)
+      for (vsip::length_type k = 0; k != v.size(2); ++k)
+	if (!equal(v.get(i, j, k), w.get(i, j, k)))
+	  return false;
+  return true;
+}
 
 /// Use a variable.  Useful for tests that must create a variable but
 /// do not otherwise use it.