Index: operation.hpp
===================================================================
RCS file: /cvsroot/boost/boost/boost/numeric/ublas/operation.hpp,v
retrieving revision 1.19
diff -u -r1.19 operation.hpp
--- operation.hpp	23 Mar 2005 20:26:57 -0000	1.19
+++ operation.hpp	11 Jun 2005 22:18:26 -0000
@@ -512,7 +512,36 @@
     axpy_prod (const matrix_expression<E1> &e1,
                const matrix_expression<E2> &e2,
                M &m, TRI,
-               dense_proxy_tag, row_major_tag) {
+               dense_proxy_tag, row_major_tag, column_major_tag) {
+        typedef M matrix_type;
+        typedef const E1 expression1_type;
+        typedef const E2 expression2_type;
+        typedef typename M::size_type size_type;
+        typedef typename M::value_type value_type;
+
+#if BOOST_UBLAS_TYPE_CHECK
+        matrix<value_type, row_major> cm (m);
+        typedef typename type_traits<value_type>::real_type real_type;
+        real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
+        indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
+#endif
+        size_type size1 (m.size1 ());
+        size_type size2 (m.size2 ());
+        for (size_type i = 0; i < size1; ++ i)
+            for (size_type j = 0; j < size2; ++ j)
+ 			    m(i,j) += inner_prod( row( e1 (), i ), column (e2 (), j) );
+#if BOOST_UBLAS_TYPE_CHECK
+        BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
+#endif
+        return m;
+    }
+    template<class M, class E1, class E2, class TRI>
+    BOOST_UBLAS_INLINE
+    M &
+    axpy_prod (const matrix_expression<E1> &e1,
+               const matrix_expression<E2> &e2,
+               M &m, TRI,
+               dense_proxy_tag, row_major_tag, row_major_tag) {
         typedef M matrix_type;
         typedef const E1 expression1_type;
         typedef const E2 expression2_type;
@@ -535,6 +564,18 @@
 #endif
         return m;
     }
+    // Dispatcher
+    template<class M, class E1, class E2, class TRI>
+    BOOST_UBLAS_INLINE
+    M &
+    axpy_prod (const matrix_expression<E1> &e1,
+               const matrix_expression<E2> &e2,
+               M &m, TRI,
+               dense_proxy_tag, row_major_tag) {
+	  return axpy_prod( e1, e2, m, TRI(), dense_proxy_tag(), row_major_tag(), 
+						typename E2::orientation_category () );
+	}
+
     template<class M, class E1, class E2, class TRI>
     BOOST_UBLAS_INLINE
     M &
@@ -591,7 +632,7 @@
     axpy_prod (const matrix_expression<E1> &e1,
                const matrix_expression<E2> &e2,
                M &m, TRI,
-               dense_proxy_tag, column_major_tag) {
+               dense_proxy_tag, column_major_tag, column_major_tag) {
         typedef M matrix_type;
         typedef const E1 expression1_type;
         typedef const E2 expression2_type;
@@ -614,12 +655,52 @@
 #endif
         return m;
     }
+  template<class M, class E1, class E2, class TRI>
+    BOOST_UBLAS_INLINE
+    M &
+    axpy_prod (const matrix_expression<E1> &e1,
+               const matrix_expression<E2> &e2,
+               M &m, TRI,
+               dense_proxy_tag, column_major_tag, row_major_tag) {
+        typedef M matrix_type;
+        typedef const E1 expression1_type;
+        typedef const E2 expression2_type;
+        typedef typename M::size_type size_type;
+        typedef typename M::value_type value_type;
+
+#if BOOST_UBLAS_TYPE_CHECK
+        matrix<value_type, column_major> cm (m);
+        typedef typename type_traits<value_type>::real_type real_type;
+        real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
+        indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
+#endif
+        size_type size1 (m.size1 ());
+        size_type size2 (m.size2 ());
+        for (size_type j = 0; j < size2; ++ j)
+            for (size_type i = 0; i < size1; ++ i)
+			    m(i, j) += inner_prod( row( e1(), i ), column( e2(), j ) );
+#if BOOST_UBLAS_TYPE_CHECK
+        BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
+#endif
+        return m;
+    }
+    // Dispatcher
     template<class M, class E1, class E2, class TRI>
     BOOST_UBLAS_INLINE
     M &
     axpy_prod (const matrix_expression<E1> &e1,
                const matrix_expression<E2> &e2,
                M &m, TRI,
+               dense_proxy_tag, column_major_tag) {
+	  return axpy_prod( e1, e2, m, TRI(), dense_proxy_tag(), column_major_tag(), 
+						typename E1::orientation_category () );
+	}
+  template<class M, class E1, class E2, class TRI>
+    BOOST_UBLAS_INLINE
+    M &
+    axpy_prod (const matrix_expression<E1> &e1,
+               const matrix_expression<E2> &e2,
+               M &m, TRI,
                sparse_proxy_tag, column_major_tag) {
         typedef M matrix_type;
         typedef TRI triangular_restriction;