max@0: // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
max@0: // Copyright (C) 2008-2011 Conrad Sanderson
max@0: // 
max@0: // This file is part of the Armadillo C++ library.
max@0: // It is provided without any warranty of fitness
max@0: // for any purpose. You can redistribute this file
max@0: // and/or modify it under the terms of the GNU
max@0: // Lesser General Public License (LGPL) as published
max@0: // by the Free Software Foundation, either version 3
max@0: // of the License or (at your option) any later version.
max@0: // (see http://www.opensource.org/licenses for more info)
max@0: 
max@0: 
max@0: //! \addtogroup op_dot
max@0: //! @{
max@0: 
max@0: 
max@0: 
max@0: 
max@0: //! for two arrays, generic version
max@0: template<typename eT>
max@0: arma_hot
max@0: arma_pure
max@0: inline
max@0: eT
max@0: op_dot::direct_dot_arma(const uword n_elem, const eT* const A, const eT* const B)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   eT val1 = eT(0);
max@0:   eT val2 = eT(0);
max@0:   
max@0:   uword i, j;
max@0:   
max@0:   for(i=0, j=1; j<n_elem; i+=2, j+=2)
max@0:     {
max@0:     val1 += A[i] * B[i];
max@0:     val2 += A[j] * B[j];
max@0:     }
max@0:   
max@0:   if(i < n_elem)
max@0:     {
max@0:     val1 += A[i] * B[i];
max@0:     }
max@0:   
max@0:   return val1 + val2;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! for two arrays, float and double version
max@0: template<typename eT>
max@0: arma_hot
max@0: arma_pure
max@0: inline
max@0: typename arma_float_only<eT>::result
max@0: op_dot::direct_dot(const uword n_elem, const eT* const A, const eT* const B)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   if( n_elem <= (128/sizeof(eT)) )
max@0:     {
max@0:     return op_dot::direct_dot_arma(n_elem, A, B);
max@0:     }
max@0:   else
max@0:     {
max@0:     #if defined(ARMA_USE_ATLAS)
max@0:       {
max@0:       arma_extra_debug_print("atlas::cblas_dot()");
max@0:       
max@0:       return atlas::cblas_dot(n_elem, A, B);
max@0:       }
max@0:     #elif defined(ARMA_USE_BLAS)
max@0:       {
max@0:       arma_extra_debug_print("blas::dot()");
max@0:       
max@0:       return blas::dot(n_elem, A, B);
max@0:       }
max@0:     #else
max@0:       {
max@0:       return op_dot::direct_dot_arma(n_elem, A, B);
max@0:       }
max@0:     #endif
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! for two arrays, complex version
max@0: template<typename eT>
max@0: inline
max@0: arma_hot
max@0: arma_pure
max@0: typename arma_cx_only<eT>::result
max@0: op_dot::direct_dot(const uword n_elem, const eT* const A, const eT* const B)
max@0:   {
max@0:   #if defined(ARMA_USE_ATLAS)
max@0:     {
max@0:     arma_extra_debug_print("atlas::cx_cblas_dot()");
max@0:     
max@0:     return atlas::cx_cblas_dot(n_elem, A, B);
max@0:     }
max@0:   #elif defined(ARMA_USE_BLAS)
max@0:     {
max@0:     // TODO: work out the mess with zdotu() and zdotu_sub() in BLAS
max@0:     return op_dot::direct_dot_arma(n_elem, A, B);
max@0:     }
max@0:   #else
max@0:     {
max@0:     return op_dot::direct_dot_arma(n_elem, A, B);
max@0:     }
max@0:   #endif
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! for two arrays, integral version
max@0: template<typename eT>
max@0: arma_hot
max@0: arma_pure
max@0: inline
max@0: typename arma_integral_only<eT>::result
max@0: op_dot::direct_dot(const uword n_elem, const eT* const A, const eT* const B)
max@0:   {
max@0:   return op_dot::direct_dot_arma(n_elem, A, B);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: 
max@0: //! for three arrays
max@0: template<typename eT>
max@0: arma_hot
max@0: arma_pure
max@0: inline
max@0: eT
max@0: op_dot::direct_dot(const uword n_elem, const eT* const A, const eT* const B, const eT* C)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   eT val = eT(0);
max@0:   
max@0:   for(uword i=0; i<n_elem; ++i)
max@0:     {
max@0:     val += A[i] * B[i] * C[i];
max@0:     }
max@0: 
max@0:   return val;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1, typename T2>
max@0: arma_hot
max@0: arma_inline
max@0: typename T1::elem_type
max@0: op_dot::apply(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   if( (is_Mat<T1>::value == true) && (is_Mat<T2>::value == true) )
max@0:     {
max@0:     return op_dot::apply_unwrap(X,Y);
max@0:     }
max@0:   else
max@0:     {
max@0:     return op_dot::apply_proxy(X,Y);
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1, typename T2>
max@0: arma_hot
max@0: arma_inline
max@0: typename T1::elem_type
max@0: op_dot::apply_unwrap(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   
max@0:   const unwrap<T1> tmp1(X.get_ref());
max@0:   const unwrap<T2> tmp2(Y.get_ref());
max@0:   
max@0:   const Mat<eT>& A = tmp1.M;
max@0:   const Mat<eT>& B = tmp2.M;
max@0:   
max@0:   arma_debug_check( (A.n_elem != B.n_elem), "dot(): objects must have the same number of elements" );
max@0:   
max@0:   return op_dot::direct_dot(A.n_elem, A.mem, B.mem);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1, typename T2>
max@0: arma_hot
max@0: inline
max@0: typename T1::elem_type
max@0: op_dot::apply_proxy(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type      eT;
max@0:   typedef typename Proxy<T1>::ea_type ea_type1;
max@0:   typedef typename Proxy<T2>::ea_type ea_type2;
max@0:   
max@0:   const Proxy<T1> A(X.get_ref());
max@0:   const Proxy<T2> B(Y.get_ref());
max@0:   
max@0:   const bool prefer_at_accessor = (Proxy<T1>::prefer_at_accessor) && (Proxy<T2>::prefer_at_accessor);
max@0:   
max@0:   if(prefer_at_accessor == false)
max@0:     {
max@0:     arma_debug_check( (A.get_n_elem() != B.get_n_elem()), "dot(): objects must have the same number of elements" );
max@0:   
max@0:     const uword    N  = A.get_n_elem();
max@0:           ea_type1 PA = A.get_ea();
max@0:           ea_type2 PB = B.get_ea();
max@0:     
max@0:     eT val1 = eT(0);
max@0:     eT val2 = eT(0);
max@0:     
max@0:     uword i,j;
max@0:     
max@0:     for(i=0, j=1; j<N; i+=2, j+=2)
max@0:       {
max@0:       val1 += PA[i] * PB[i];
max@0:       val2 += PA[j] * PB[j];
max@0:       }
max@0:     
max@0:     if(i < N)
max@0:       {
max@0:       val1 += PA[i] * PB[i];
max@0:       }
max@0:     
max@0:     return val1 + val2;
max@0:     }
max@0:   else
max@0:     {
max@0:     return op_dot::apply_unwrap(A.Q, B.Q);
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // op_norm_dot
max@0: 
max@0: 
max@0: 
max@0: template<typename T1, typename T2>
max@0: arma_hot
max@0: inline
max@0: typename T1::elem_type
max@0: op_norm_dot::apply(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type      eT;
max@0:   typedef typename Proxy<T1>::ea_type ea_type1;
max@0:   typedef typename Proxy<T2>::ea_type ea_type2;
max@0:   
max@0:   const bool prefer_at_accessor = (Proxy<T1>::prefer_at_accessor) && (Proxy<T2>::prefer_at_accessor);
max@0:   
max@0:   if(prefer_at_accessor == false)
max@0:     {
max@0:     const Proxy<T1> A(X.get_ref());
max@0:     const Proxy<T2> B(Y.get_ref());
max@0:     
max@0:     arma_debug_check( (A.get_n_elem() != B.get_n_elem()), "norm_dot(): objects must have the same number of elements" );
max@0:     
max@0:     const uword    N  = A.get_n_elem();
max@0:           ea_type1 PA = A.get_ea();
max@0:           ea_type2 PB = B.get_ea();
max@0:     
max@0:     eT acc1 = eT(0);
max@0:     eT acc2 = eT(0);
max@0:     eT acc3 = eT(0);
max@0:     
max@0:     for(uword i=0; i<N; ++i)
max@0:       {
max@0:       const eT tmpA = PA[i];
max@0:       const eT tmpB = PB[i];
max@0:       
max@0:       acc1 += tmpA * tmpA;
max@0:       acc2 += tmpB * tmpB;
max@0:       acc3 += tmpA * tmpB;
max@0:       }
max@0:       
max@0:     return acc3 / ( std::sqrt(acc1 * acc2) );
max@0:     }
max@0:   else
max@0:     {
max@0:     return op_norm_dot::apply_unwrap(X, Y);
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1, typename T2>
max@0: arma_hot
max@0: inline
max@0: typename T1::elem_type
max@0: op_norm_dot::apply_unwrap(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   
max@0:   const unwrap<T1> tmp1(X.get_ref());
max@0:   const unwrap<T2> tmp2(Y.get_ref());
max@0:   
max@0:   const Mat<eT>& A = tmp1.M;
max@0:   const Mat<eT>& B = tmp2.M;
max@0:   
max@0:   
max@0:   arma_debug_check( (A.n_elem != B.n_elem), "norm_dot(): objects must have the same number of elements" );
max@0:   
max@0:   const uword N = A.n_elem;
max@0:   
max@0:   const eT* A_mem = A.memptr();
max@0:   const eT* B_mem = B.memptr();
max@0:   
max@0:   eT acc1 = eT(0);
max@0:   eT acc2 = eT(0);
max@0:   eT acc3 = eT(0);
max@0:   
max@0:   for(uword i=0; i<N; ++i)
max@0:     {
max@0:     const eT tmpA = A_mem[i];
max@0:     const eT tmpB = B_mem[i];
max@0:     
max@0:     acc1 += tmpA * tmpA;
max@0:     acc2 += tmpB * tmpB;
max@0:     acc3 += tmpA * tmpB;
max@0:     }
max@0:     
max@0:   return acc3 / ( std::sqrt(acc1 * acc2) );
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //
max@0: // op_cdot
max@0: 
max@0: 
max@0: 
max@0: template<typename T1, typename T2>
max@0: arma_hot
max@0: arma_inline
max@0: typename T1::elem_type
max@0: op_cdot::apply(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type      eT;
max@0:   typedef typename Proxy<T1>::ea_type ea_type1;
max@0:   typedef typename Proxy<T2>::ea_type ea_type2;
max@0:   
max@0:   const Proxy<T1> A(X.get_ref());
max@0:   const Proxy<T2> B(Y.get_ref());
max@0:   
max@0:   arma_debug_check( (A.get_n_elem() != B.get_n_elem()), "cdot(): objects must have the same number of elements" );
max@0:   
max@0:   const uword    N  = A.get_n_elem();
max@0:         ea_type1 PA = A.get_ea();
max@0:         ea_type2 PB = B.get_ea();
max@0:   
max@0:   eT val1 = eT(0);
max@0:   eT val2 = eT(0);
max@0:   
max@0:   uword i,j;
max@0:   for(i=0, j=1; j<N; i+=2, j+=2)
max@0:     {
max@0:     val1 += std::conj(PA[i]) * PB[i];
max@0:     val2 += std::conj(PA[j]) * PB[j];
max@0:     }
max@0:   
max@0:   if(i < N)
max@0:     {
max@0:     val1 += std::conj(PA[i]) * PB[i];
max@0:     }
max@0:   
max@0:   return val1 + val2;
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! @}