Ticket #4917: ublas-pow

Index: boost/numeric/ublas/vector_expression.hpp
===================================================================
--- boost/numeric/ublas/vector_expression.hpp   (revision 66890)
+++ boost/numeric/ublas/vector_expression.hpp   (working copy)
@@ -950,6 +950,17 @@
         return expression_type (e1 (), e2 ());
     }

+    // pow (v1, v2) [i] = pow( v1 [i], v2 [i] )
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    typename vector_binary_traits<E1, E2, scalar_pow<typename E1::value_type,
+                                                     typename E2::value_type> >::result_type
+    element_pow (const vector_expression<E1> &e1,
+                 const vector_expression<E2> &e2) {
+        typedef typename vector_binary_traits<E1, E2, scalar_pow<typename E1::value_type,
+                                                                 typename E2::value_type> >::expression_type expression_type;
+        return expression_type (e1 (), e2 ());
+    }

     template<class E1, class E2, class F>
     class vector_binary_scalar1:
@@ -1416,6 +1427,15 @@
         return expression_type (e1 (), e2);
     }

+    // pow (v, t) [i] = pow( v [i], t )
+    template<class E1, class T2>
+    BOOST_UBLAS_INLINE
+    typename vector_binary_scalar2_traits<E1, const T2, scalar_pow<typename E1::value_type, T2> >::result_type
+    pow (const vector_expression<E1> &e1,
+                const T2 &e2) {
+        typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_pow<typename E1::value_type, T2> >::expression_type expression_type;
+        return expression_type (e1 (), e2);
+    }

     template<class E, class F>
     class vector_scalar_unary:
Index: boost/numeric/ublas/matrix_expression.hpp
===================================================================
--- boost/numeric/ublas/matrix_expression.hpp   (revision 66890)
+++ boost/numeric/ublas/matrix_expression.hpp   (working copy)
@@ -2518,6 +2518,18 @@
         return expression_type (e1 (), e2 ());
     }

+    // pow (m1, m2) [i] [j] = pow (m1 [i] [j], m2 [i] [j])
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    typename matrix_binary_traits<E1, E2, scalar_pow<typename E1::value_type,
+                                                     typename E2::value_type> >::result_type
+    element_pow (const matrix_expression<E1> &e1,
+                 const matrix_expression<E2> &e2) {
+        typedef typename matrix_binary_traits<E1, E2, scalar_pow<typename E1::value_type,
+                                                                 typename E2::value_type> >::expression_type expression_type;
+        return expression_type (e1 (), e2 ());
+    }
+
     template<class E1, class E2, class F>
     class matrix_binary_scalar1:
         public matrix_expression<matrix_binary_scalar1<E1, E2, F> > {
@@ -3394,6 +3406,15 @@
         return expression_type (e1 (), e2);
     }

+    // pow (m, p) [i] [j] = pow (m [i] [j], p)
+    template<class E1, class T2>
+    BOOST_UBLAS_INLINE
+    typename matrix_binary_scalar2_traits<E1, const T2, scalar_pow<typename E1::value_type, T2> >::result_type
+    pow (const matrix_expression<E1> &e1,
+                const T2 &e2) {
+        typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_pow<typename E1::value_type, T2> >::expression_type expression_type;
+        return expression_type (e1 (), e2);
+    }

     template<class E1, class E2, class F>
     class matrix_vector_binary1:
Index: boost/numeric/ublas/functional.hpp
===================================================================
--- boost/numeric/ublas/functional.hpp  (revision 66890)
+++ boost/numeric/ublas/functional.hpp  (working copy)
@@ -172,7 +172,19 @@
             return t1 / t2;
         }
     };
+    template<class T, class U>
+    struct scalar_pow:
+        public scalar_binary_functor<T, U> {
+        typedef typename scalar_binary_functor<T, U>::argument1_type argument1_type;
+        typedef typename scalar_binary_functor<T, U>::argument2_type argument2_type;
+        typedef typename boost::remove_reference< argument1_type >::type result_type;

+        static BOOST_UBLAS_INLINE
+        result_type apply (argument1_type t, argument2_type p) {
+            return type_traits<T>::type_pow (t,p);
+        }
+    };
+
     template<class T1, class T2>
     struct scalar_binary_assign_functor {
         // ISSUE Remove reference to avoid reference to reference problems
Index: boost/numeric/ublas/traits.hpp
===================================================================
--- boost/numeric/ublas/traits.hpp      (revision 66890)
+++ boost/numeric/ublas/traits.hpp      (working copy)
@@ -41,6 +41,11 @@
     // we'll find either std::abs or else another version via ADL:
     return abs (t);
   }
+  template<class T, class U> T boost_numeric_ublas_pow (const T& t, const U& p) {
+    using namespace std;
+    // we'll find either std::pow or else another version via ADL:
+    return pow (t,p);
+  }
 }

 namespace boost { namespace numeric { namespace ublas {
@@ -186,6 +191,11 @@
         real_type type_abs (const_reference t) {
             return boost_numeric_ublas_abs (t);
         }
+        template<class U> static
+        BOOST_UBLAS_INLINE
+        value_type type_pow (const_reference t, const U& p) {
+            return boost_numeric_ublas_pow (t,p);
+        }
         static
         BOOST_UBLAS_INLINE
         value_type type_sqrt (const_reference t) {
@@ -297,6 +307,11 @@
         real_type type_abs (const_reference t) {
                 return abs (t);
         }
+        template<class U> static
+        BOOST_UBLAS_INLINE
+        value_type type_pow (const_reference t, const U& p) {
+            return boost_numeric_ublas_pow (t,p);
+        }
         static
         BOOST_UBLAS_INLINE
         value_type type_sqrt (const_reference t) {
Index: libs/numeric/ublas/test/test_pow.cpp
===================================================================
--- libs/numeric/ublas/test/test_pow.cpp        (revision 0)
+++ libs/numeric/ublas/test/test_pow.cpp        (revision 0)
@@ -0,0 +1,118 @@
+#include "test_pow.hpp"
+
+template<class V, int N>
+struct test_vector_pow {
+    typedef typename V::value_type value_type;
+    typedef typename V::size_type size_type;
+    typedef typename ublas::type_traits<value_type>::real_type real_type;
+
+    template<class VP>
+    void test_pow_with (VP &v1, VP &v2, VP &v3) const {
+        initialize_vector (v1);
+        v2 = v1;
+        std::cout << "v2 = v1 = " << v1 << std::endl;
+        v2 = pow( v1, 2 );
+        std::cout << "pow( v1, 2 ) = " << v2 << std::endl;
+        v1 = pow( v2, 0.5 );
+        std::cout << "pow( v2, 0.5 ) = " << v1 << std::endl;
+
+        initialize_vector (v1);
+        v2 = v1;
+        std::cout << "v2 = v1 = " << v1 << std::endl;
+        v3 = element_pow( v1, v2 );
+        std::cout << "element_pow( v1, v2 ) = " << v3 << std::endl;
+    }
+
+    void operator () () const {
+        V v1 (N);
+        V v2 (N);
+        V v3 (N);
+        test_pow_with (v1, v2, v3);
+    }
+};
+
+template<class M, int N>
+struct test_matrix_pow {
+    typedef typename M::value_type value_type;
+
+    template<class MP>
+    void test_pow_with (MP &m1, MP &m2, MP &m3) const {
+        initialize_matrix (m1);
+        m2 = m1;
+        std::cout << "m2 = m1 = " << m1 << std::endl;
+        m2 = pow( m1, 2 );
+        std::cout << "pow( m1, 2 ) = " << m2 << std::endl;
+        m1 = pow( m2, 0.5 );
+        std::cout << "pow( m2, 0.5 ) = " << m1 << std::endl;
+
+        initialize_matrix (m1);
+        m2 = m1;
+        std::cout << "m2 = m1 = " << m1 << std::endl;
+        m3 = element_pow( m1, m2 );
+        std::cout << "element_pow( m1, m2 ) = " << m3 << std::endl;
+    }
+    void operator () () const {
+        M m1 (N, N);
+        M m2 (N, N);
+        M m3 (N, N);
+        test_pow_with (m1, m2, m3);
+    }
+};
+
+void test_vector () {
+    std::cout << "test_vector" << std::endl;
+
+#ifdef USE_FLOAT
+    std::cout << "float" << std::endl;
+    test_vector_pow<ublas::vector<float>, 3> () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "double" << std::endl;
+    test_vector_pow<ublas::vector<double>, 3> () ();
+#endif
+
+#ifdef USE_STD_COMPLEX
+#ifdef USE_FLOAT
+    std::cout << "std::complex<float>" << std::endl;
+    test_vector_pow<ublas::vector<std::complex<float> >, 3> () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "std::complex<double>" << std::endl;
+    test_vector_pow<ublas::vector<std::complex<double> >, 3> () ();
+#endif
+#endif
+}
+
+void test_matrix () {
+    std::cout << "test_matrix" << std::endl;
+
+#ifdef USE_FLOAT
+    std::cout << "float" << std::endl;
+    test_matrix_pow<ublas::matrix<float>, 3> () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "double" << std::endl;
+    test_matrix_pow<ublas::matrix<double>, 3> () ();
+#endif
+
+#ifdef USE_STD_COMPLEX
+#ifdef USE_FLOAT
+    std::cout << "std::complex<float>" << std::endl;
+    test_matrix_pow<ublas::matrix<std::complex<float> >, 3> () ();
+#endif
+
+#ifdef USE_DOUBLE
+    std::cout << "std::complex<double>" << std::endl;
+    test_matrix_pow<ublas::matrix<std::complex<double> >, 3> () ();
+#endif
+#endif
+}
+
+int main () {
+    test_vector ();
+    test_matrix ();
+    return 0;
+}
Index: libs/numeric/ublas/test/test_pow.hpp
===================================================================
--- libs/numeric/ublas/test/test_pow.hpp        (revision 0)
+++ libs/numeric/ublas/test/test_pow.hpp        (revision 0)
@@ -0,0 +1,20 @@
+#ifndef TEST_POW_H
+#define TEST_POW_H
+
+#include <iostream>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/io.hpp>
+
+namespace ublas = boost::numeric::ublas;
+
+#include "common/init.hpp"
+
+void test_vector ();
+void test_matrix ();
+
+
+#endif
Index: libs/numeric/ublas/test/Jamfile.v2
===================================================================
--- libs/numeric/ublas/test/Jamfile.v2  (revision 66890)
+++ libs/numeric/ublas/test/Jamfile.v2  (working copy)
@@ -132,6 +132,10 @@
       ]
       [ run test_complex_norms.cpp
       ]
+      [ run test_pow.cpp
+        : : :
+            <define>$(UBLAS_TESTSET)
+      ]
        [ run test_assignment.cpp
        ]
     ;