[patch] matvec: remaining prod functions

[Don McCoy <don@xxxxxxxxxxxxxxxx>]

The attached implements the last of the functions needed for [math.matvec].

--
Don McCoy
CodeSourcery, LLC

2005-09-28  Don McCoy  <don@xxxxxxxxxxxxxxxx>
	
	* src/vsip/impl/matvec-prod.hpp: added prod3, prod4,
	  prodh, prodj and prodt.
	* tests/matvec-prod.cpp: added tests for same.
	

Index: src/vsip/impl/matvec-prod.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/impl/matvec-prod.hpp,v
retrieving revision 1.1
diff -c -p -r1.1 matvec-prod.hpp
*** src/vsip/impl/matvec-prod.hpp	13 Sep 2005 16:39:45 -0000	1.1
--- src/vsip/impl/matvec-prod.hpp	28 Sep 2005 17:30:50 -0000
***************
*** 16,22 ****
  
  #include <vsip/vector.hpp>
  #include <vsip/matrix.hpp>
! 
  
  
  /***********************************************************************
--- 16,22 ----
  
  #include <vsip/vector.hpp>
  #include <vsip/matrix.hpp>
! #include <vsip/impl/matvec.hpp>
  
  
  /***********************************************************************
*************** prod(
*** 187,192 ****
--- 187,355 ----
    return r;
  }
  
+ 
+ 
+ /// [3x3] Matrix-matrix product dispatch.
+ 
+ template <typename T0,
+ 	  typename T1,
+ 	  typename Block0,
+ 	  typename Block1>
+ const_Matrix<typename Promotion<T0, T1>::type>
+ prod3(
+   const_Matrix<T0, Block0> a,
+   const_Matrix<T1, Block1> b)
+ {
+   assert( a.size(0) == 3 );
+   assert( a.size(1) == 3 );
+   typedef typename Promotion<T0, T1>::type return_type;
+ 
+   Matrix<return_type> r(a.size(0), b.size(1));
+ 
+   impl::generic_prod(a, b, r);
+ 
+   return r;
+ }
+ 
+ 
+ /// [3x3] Matrix-vector product dispatch.
+ 
+ template <typename T0,
+ 	  typename T1,
+ 	  typename Block0,
+ 	  typename Block1>
+ const_Vector<typename Promotion<T0, T1>::type>
+ prod3(
+   const_Matrix<T0, Block0> a,
+   const_Vector<T1, Block1> b)
+ {
+   assert( a.size(0) == 3 );
+   assert( a.size(1) == 3 );
+   typedef typename Promotion<T0, T1>::type return_type;
+ 
+   Vector<return_type> r(a.size(0));
+ 
+   impl::generic_prod(a, b, r);
+ 
+   return r;
+ }
+ 
+ 
+ 
+ /// [4x4] Matrix-matrix product dispatch.
+ 
+ template <typename T0,
+ 	  typename T1,
+ 	  typename Block0,
+ 	  typename Block1>
+ const_Matrix<typename Promotion<T0, T1>::type>
+ prod4(
+   const_Matrix<T0, Block0> a,
+   const_Matrix<T1, Block1> b)
+ {
+   assert( a.size(0) == 4 );
+   assert( a.size(1) == 4 );
+   typedef typename Promotion<T0, T1>::type return_type;
+ 
+   Matrix<return_type> r(a.size(0), b.size(1));
+ 
+   impl::generic_prod(a, b, r);
+ 
+   return r;
+ }
+ 
+ 
+ 
+ /// [4x4] Matrix-vector product dispatch.
+ 
+ template <typename T0,
+ 	  typename T1,
+ 	  typename Block0,
+ 	  typename Block1>
+ const_Vector<typename Promotion<T0, T1>::type>
+ prod4(
+   const_Matrix<T0, Block0> a,
+   const_Vector<T1, Block1> b)
+ {
+   assert( a.size(0) == 4 );
+   assert( a.size(1) == 4 );
+   typedef typename Promotion<T0, T1>::type return_type;
+ 
+   Vector<return_type> r(a.size(0));
+ 
+   impl::generic_prod(a, b, r);
+ 
+   return r;
+ }
+ 
+ 
+ 
+ /// Matrix-Matrix (with hermitian) product dispatch.
+ 
+ template <typename T0,
+           typename T1,
+           typename Block0,
+           typename Block1>
+ const_Matrix<typename Promotion<complex<T0>, complex<T1> >::type>
+ prodh(
+   const_Matrix<complex<T0>, Block0> m0,
+   const_Matrix<complex<T1>, Block1> m1) 
+     VSIP_NOTHROW
+ {
+   typedef typename Promotion<complex<T0>, complex<T1> >::type return_type;
+ 
+   Matrix<return_type> r(m0.size(0), m1.size(1));
+ 
+   impl::generic_prod(m0, herm(m1), r);
+ 
+   return r;
+ }
+ 
+ 
+ /// Matrix-Matrix (with complex conjugate) product dispatch.
+ 
+ template <typename T0,
+           typename T1,
+           typename Block0,
+           typename Block1>
+ const_Matrix<typename Promotion<complex<T0>, complex<T1> >::type>
+ prodj(
+   const_Matrix<complex<T0>, Block0> m0,
+   const_Matrix<complex<T1>, Block1> m1)
+     VSIP_NOTHROW
+ {
+   typedef typename Promotion<T0, T1>::type return_type;
+ 
+   Matrix<return_type> r(m0.size(0), m1.size(1));
+ 
+   impl::generic_prod(m0, conj(m1), r);
+ 
+   return r;
+ }
+ 
+ 
+ /// Matrix-Matrix (with transpose) product dispatch.
+ 
+ template <typename T0,
+           typename T1,
+           typename Block0,
+           typename Block1>
+ const_Matrix<typename Promotion<T0, T1>::type>
+ prodt(
+   const_Matrix<T0, Block0> m0,
+   const_Matrix<T1, Block1> m1)
+     VSIP_NOTHROW
+ {
+   typedef typename Promotion<T0, T1>::type return_type;
+ 
+   Matrix<return_type> r(m0.size(0), m1.size(1));
+ 
+   impl::generic_prod(m0, trans(m1), r);
+ 
+   return r;
+ }
+ 
+ 
  } // namespace vsip
  
  #endif // VSIP_IMPL_MATVEC_PROD_HPP
Index: tests/matvec-prod.cpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/matvec-prod.cpp,v
retrieving revision 1.1
diff -c -p -r1.1 matvec-prod.cpp
*** tests/matvec-prod.cpp	13 Sep 2005 16:39:45 -0000	1.1
--- tests/matvec-prod.cpp	28 Sep 2005 17:30:50 -0000
*************** ref_prod(
*** 103,108 ****
--- 103,133 ----
  }
  
  
+ template <typename T0,
+ 	  typename T1,
+           typename T2>
+ void
+ check_prod( Matrix<T0> test, Matrix<T1> chk, Matrix<T2> gauge )
+ {
+   typedef typename Promotion<T0, T1>::type return_type;
+   typedef typename vsip::impl::Scalar_of<return_type>::type scalar_type;
+ 
+   Index<2> idx;
+   scalar_type err = maxval(((mag(chk - test)
+ 			     / Precision_traits<scalar_type>::eps)
+ 			    / gauge),
+ 			   idx);
+ 
+ #if VERBOSE
+   cout << "test  =\n" << test;
+   cout << "chk   =\n" << chk;
+   cout << "gauge =\n" << gauge;
+   cout << "err = " << err << endl;
+ #endif
+ 
+   assert(err < 10.0);
+ }
+ 
  
  /***********************************************************************
    Reference Definitions
*************** test_prod_rand(length_type m, length_typ
*** 148,181 ****
        if (!(gauge(i, j) > scalar_type()))
  	gauge(i, j) = scalar_type(1);
  
!   Index<2> idx;
!   scalar_type err1 = maxval(((mag(chk - res1)
! 			     / Precision_traits<scalar_type>::eps)
! 			    / gauge),
! 			   idx);
!   scalar_type err2 = maxval(((mag(chk - res2)
! 			     / Precision_traits<scalar_type>::eps)
! 			    / gauge),
! 			   idx);
!   scalar_type err3 = maxval(((mag(chk - res3)
! 			     / Precision_traits<scalar_type>::eps)
! 			    / gauge),
! 			   idx);
  
  #if VERBOSE
    cout << "a     =\n" << a;
    cout << "b     =\n" << b;
-   cout << "res   =\n" << res;
-   cout << "chk   =\n" << chk;
-   cout << "gauge =\n" << gauge;
-   cout << "err1 = " << err1 << "  "
-        << "err2 = " << err2 << "  "
-        << "err3 = " << err3 << endl;
  #endif
  
!   assert(err1 < 10.0);
!   assert(err2 < 10.0);
!   assert(err3 < 10.0);
  }
  
  
--- 173,393 ----
        if (!(gauge(i, j) > scalar_type()))
  	gauge(i, j) = scalar_type(1);
  
! #if VERBOSE
!   cout << "a     =\n" << a;
!   cout << "b     =\n" << b;
! #endif
! 
!   check_prod( res1, chk, gauge );
!   check_prod( res2, chk, gauge );
!   check_prod( res3, chk, gauge );
! }
! 
! 
! 
! template <typename T0,
! 	  typename T1>
! void
! test_prod3_rand()
! {
!   typedef typename Promotion<T0, T1>::type return_type;
!   typedef typename vsip::impl::Scalar_of<return_type>::type scalar_type;
!   const length_type m = 3;
!   const length_type n = 3;
!   const length_type k = 3;
! 
!   Matrix<T0> a(m, n);
!   Matrix<T1> b(n, k);
!   Matrix<return_type> res1(m, k);
!   Matrix<return_type> res2(m, k);
!   Matrix<return_type> chk(m, k);
!   Matrix<scalar_type> gauge(m, k);
! 
!   randm(a);
!   randm(b);
! 
!   // Test matrix-matrix prod
!   res1   = prod3(a, b);
! 
!   // Test matrix-vector prod
!   for (index_type i=0; i<k; ++i)
!     res2.col(i) = prod3(a, b.col(i));
! 
!   chk   = ref_prod(a, b);
!   gauge = ref_prod(mag(a), mag(b));
! 
!   for (index_type i=0; i<gauge.size(0); ++i)
!     for (index_type j=0; j<gauge.size(1); ++j)
!       if (!(gauge(i, j) > scalar_type()))
! 	gauge(i, j) = scalar_type(1);
! 
! #if VERBOSE
!   cout << "a     =\n" << a;
!   cout << "b     =\n" << b;
! #endif
! 
!   check_prod( res1, chk, gauge );
!   check_prod( res2, chk, gauge );
! }
! 
! 
! template <typename T0,
! 	  typename T1>
! void
! test_prod4_rand()
! {
!   typedef typename Promotion<T0, T1>::type return_type;
!   typedef typename vsip::impl::Scalar_of<return_type>::type scalar_type;
!   const length_type m = 4;
!   const length_type n = 4;
!   const length_type k = 4;
! 
!   Matrix<T0> a(m, n);
!   Matrix<T1> b(n, k);
!   Matrix<return_type> res1(m, k);
!   Matrix<return_type> res2(m, k);
!   Matrix<return_type> chk(m, k);
!   Matrix<scalar_type> gauge(m, k);
! 
!   randm(a);
!   randm(b);
! 
!   // Test matrix-matrix prod
!   res1   = prod4(a, b);
! 
!   // Test matrix-vector prod
!   for (index_type i=0; i<k; ++i)
!     res2.col(i) = prod4(a, b.col(i));
! 
!   chk   = ref_prod(a, b);
!   gauge = ref_prod(mag(a), mag(b));
! 
!   for (index_type i=0; i<gauge.size(0); ++i)
!     for (index_type j=0; j<gauge.size(1); ++j)
!       if (!(gauge(i, j) > scalar_type()))
! 	gauge(i, j) = scalar_type(1);
  
  #if VERBOSE
    cout << "a     =\n" << a;
    cout << "b     =\n" << b;
  #endif
  
!   check_prod( res1, chk, gauge );
!   check_prod( res2, chk, gauge );
! }
! 
! 
! 
! template <typename T0,
! 	  typename T1>
! void
! test_prodh_rand(length_type m, length_type n, length_type k)
! {
!   typedef typename Promotion<T0, T1>::type return_type;
!   typedef typename vsip::impl::Scalar_of<return_type>::type scalar_type;
! 
!   Matrix<T0> a(m, n);
!   Matrix<T1> b(k, n);
!   Matrix<return_type> res1(m, k);
!   Matrix<return_type> chk(m, k);
!   Matrix<scalar_type> gauge(m, k);
! 
!   randm(a);
!   randm(b);
! 
!   // Test matrix-matrix prod for hermitian
!   res1   = prod(a, herm(b));
! 
!   chk   = ref_prod(a, herm(b));
!   gauge = ref_prod(mag(a), mag(herm(b)));
! 
!   for (index_type i=0; i<gauge.size(0); ++i)
!     for (index_type j=0; j<gauge.size(1); ++j)
!       if (!(gauge(i, j) > scalar_type()))
! 	gauge(i, j) = scalar_type(1);
! 
! #if VERBOSE
!   cout << "a     =\n" << a;
!   cout << "b     =\n" << b;
! #endif
! 
!   check_prod( res1, chk, gauge );
! }
! 
! 
! 
! template <typename T0,
! 	  typename T1>
! void
! test_prodj_rand(length_type m, length_type n, length_type k)
! {
!   typedef typename Promotion<T0, T1>::type return_type;
!   typedef typename vsip::impl::Scalar_of<return_type>::type scalar_type;
! 
!   Matrix<T0> a(m, n);
!   Matrix<T1> b(n, k);
!   Matrix<return_type> res1(m, k);
!   Matrix<return_type> chk(m, k);
!   Matrix<scalar_type> gauge(m, k);
! 
!   randm(a);
!   randm(b);
! 
!   // Test matrix-matrix prod for hermitian
!   res1   = prod(a, conj(b));
! 
!   chk   = ref_prod(a, conj(b));
!   gauge = ref_prod(mag(a), mag(conj(b)));
! 
!   for (index_type i=0; i<gauge.size(0); ++i)
!     for (index_type j=0; j<gauge.size(1); ++j)
!       if (!(gauge(i, j) > scalar_type()))
! 	gauge(i, j) = scalar_type(1);
! 
! #if VERBOSE
!   cout << "a     =\n" << a;
!   cout << "b     =\n" << b;
! #endif
! 
!   check_prod( res1, chk, gauge );
! }
! 
! 
! 
! template <typename T0,
! 	  typename T1>
! void
! test_prodt_rand(length_type m, length_type n, length_type k)
! {
!   typedef typename Promotion<T0, T1>::type return_type;
!   typedef typename vsip::impl::Scalar_of<return_type>::type scalar_type;
! 
!   Matrix<T0> a(m, n);
!   Matrix<T1> b(k, n);
!   Matrix<return_type> res1(m, k);
!   Matrix<return_type> chk(m, k);
!   Matrix<scalar_type> gauge(m, k);
! 
!   randm(a);
!   randm(b);
! 
!   // Test matrix-matrix prod for transpose
!   res1   = prod(a, trans(b));
! 
!   chk   = ref_prod(a, trans(b));
!   gauge = ref_prod(mag(a), mag(trans(b)));
! 
!   for (index_type i=0; i<gauge.size(0); ++i)
!     for (index_type j=0; j<gauge.size(1); ++j)
!       if (!(gauge(i, j) > scalar_type()))
! 	gauge(i, j) = scalar_type(1);
! 
! #if VERBOSE
!   cout << "a     =\n" << a;
!   cout << "b     =\n" << b;
! #endif
! 
!   check_prod( res1, chk, gauge );
  }
  
  
*************** prod_cases()
*** 189,197 ****
--- 401,433 ----
    test_prod_rand<T0, T1>(5, 7, 9);
    test_prod_rand<T0, T1>(9, 5, 7);
    test_prod_rand<T0, T1>(9, 7, 5);
+ 
+   test_prod3_rand<T0, T1>();
+   test_prod4_rand<T0, T1>();
+ 
+   test_prodt_rand<T0, T1>(5, 5, 5);
+   test_prodt_rand<T0, T1>(5, 7, 9);
+   test_prodt_rand<T0, T1>(9, 5, 7);
+   test_prodt_rand<T0, T1>(9, 7, 5);
  }
  
  
+ template <typename T0,
+ 	  typename T1>
+ void
+ prod_cases_complex_only()
+ {
+   test_prodh_rand<T0, T1>(5, 5, 5);
+   test_prodh_rand<T0, T1>(5, 7, 9);
+   test_prodh_rand<T0, T1>(9, 5, 7);
+   test_prodh_rand<T0, T1>(9, 7, 5);
+ 
+   test_prodj_rand<T0, T1>(5, 5, 5);
+   test_prodj_rand<T0, T1>(5, 7, 9);
+   test_prodj_rand<T0, T1>(9, 5, 7);
+   test_prodj_rand<T0, T1>(9, 7, 5);
+ }
+ 
  
  void
  prod_types()
*************** prod_types()
*** 204,213 ****
--- 440,453 ----
    prod_cases<complex<float>, complex<float> >();
    prod_cases<float,          complex<float> >();
    prod_cases<complex<float>, float >();
+ 
+   prod_cases_complex_only<complex<float>, complex<float> >();
  }
  
  
  
+ 
+ 
  /***********************************************************************
    Main
  ***********************************************************************/