Chris@49: // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2008-2011 Conrad Sanderson
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup gemm_mixed
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! \brief
Chris@49: //! Matrix multplication where the matrices have differing element types.
Chris@49: //! Uses caching for speedup.
Chris@49: //! Matrix 'C' is assumed to have been set to the correct size (i.e. taking into account transposes)
Chris@49: 
Chris@49: template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
Chris@49: class gemm_mixed_large
Chris@49:   {
Chris@49:   public:
Chris@49:   
Chris@49:   template<typename out_eT, typename in_eT1, typename in_eT2>
Chris@49:   arma_hot
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   apply
Chris@49:     (
Chris@49:           Mat<out_eT>& C,
Chris@49:     const Mat<in_eT1>& A,
Chris@49:     const Mat<in_eT2>& B,
Chris@49:     const out_eT alpha = out_eT(1),
Chris@49:     const out_eT beta  = out_eT(0)
Chris@49:     )
Chris@49:     {
Chris@49:     arma_extra_debug_sigprint();
Chris@49:     
Chris@49:     const uword A_n_rows = A.n_rows;
Chris@49:     const uword A_n_cols = A.n_cols;
Chris@49:     
Chris@49:     const uword B_n_rows = B.n_rows;
Chris@49:     const uword B_n_cols = B.n_cols;
Chris@49:     
Chris@49:     if( (do_trans_A == false) && (do_trans_B == false) )
Chris@49:       {
Chris@49:       podarray<in_eT1> tmp(A_n_cols);
Chris@49:       in_eT1* A_rowdata = tmp.memptr();
Chris@49:       
Chris@49:       for(uword row_A=0; row_A < A_n_rows; ++row_A)
Chris@49:         {
Chris@49:         tmp.copy_row(A, row_A);
Chris@49:         
Chris@49:         for(uword col_B=0; col_B < B_n_cols; ++col_B)
Chris@49:           {
Chris@49:           const in_eT2* B_coldata = B.colptr(col_B);
Chris@49:           
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i=0; i < B_n_rows; ++i)
Chris@49:             {
Chris@49:             acc += upgrade_val<in_eT1,in_eT2>::apply(A_rowdata[i]) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
Chris@49:             }
Chris@49:         
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = acc + beta*C.at(row_A,col_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = alpha*acc + beta*C.at(row_A,col_B);
Chris@49:             }
Chris@49:           
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:     if( (do_trans_A == true) && (do_trans_B == false) )
Chris@49:       {
Chris@49:       for(uword col_A=0; col_A < A_n_cols; ++col_A)
Chris@49:         {
Chris@49:         // col_A is interpreted as row_A when storing the results in matrix C
Chris@49:         
Chris@49:         const in_eT1* A_coldata = A.colptr(col_A);
Chris@49:         
Chris@49:         for(uword col_B=0; col_B < B_n_cols; ++col_B)
Chris@49:           {
Chris@49:           const in_eT2* B_coldata = B.colptr(col_B);
Chris@49:           
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i=0; i < B_n_rows; ++i)
Chris@49:             {
Chris@49:             acc += upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
Chris@49:             }
Chris@49:         
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = acc + beta*C.at(col_A,col_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = alpha*acc + beta*C.at(col_A,col_B);
Chris@49:             }
Chris@49:           
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:     if( (do_trans_A == false) && (do_trans_B == true) )
Chris@49:       {
Chris@49:       Mat<in_eT2> B_tmp;
Chris@49:       
Chris@49:       op_strans::apply_noalias(B_tmp, B);
Chris@49:       
Chris@49:       gemm_mixed_large<false, false, use_alpha, use_beta>::apply(C, A, B_tmp, alpha, beta);
Chris@49:       }
Chris@49:     else
Chris@49:     if( (do_trans_A == true) && (do_trans_B == true) )
Chris@49:       {
Chris@49:       // mat B_tmp = trans(B);
Chris@49:       // dgemm_arma<true, false,  use_alpha, use_beta>::apply(C, A, B_tmp, alpha, beta);
Chris@49:       
Chris@49:       
Chris@49:       // By using the trans(A)*trans(B) = trans(B*A) equivalency,
Chris@49:       // transpose operations are not needed
Chris@49:       
Chris@49:       podarray<in_eT2> tmp(B_n_cols);
Chris@49:       in_eT2* B_rowdata = tmp.memptr();
Chris@49:       
Chris@49:       for(uword row_B=0; row_B < B_n_rows; ++row_B)
Chris@49:         {
Chris@49:         tmp.copy_row(B, row_B);
Chris@49:         
Chris@49:         for(uword col_A=0; col_A < A_n_cols; ++col_A)
Chris@49:           {
Chris@49:           const in_eT1* A_coldata = A.colptr(col_A);
Chris@49:           
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i=0; i < A_n_rows; ++i)
Chris@49:             {
Chris@49:             acc += upgrade_val<in_eT1,in_eT2>::apply(B_rowdata[i]) * upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]);
Chris@49:             }
Chris@49:         
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = acc + beta*C.at(col_A,row_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = alpha*acc + beta*C.at(col_A,row_B);
Chris@49:             }
Chris@49:           
Chris@49:           }
Chris@49:         }
Chris@49:       
Chris@49:       }
Chris@49:     }
Chris@49:     
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Matrix multplication where the matrices have different element types.
Chris@49: //! Simple version (no caching).
Chris@49: //! Matrix 'C' is assumed to have been set to the correct size (i.e. taking into account transposes)
Chris@49: template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
Chris@49: class gemm_mixed_small
Chris@49:   {
Chris@49:   public:
Chris@49:   
Chris@49:   template<typename out_eT, typename in_eT1, typename in_eT2>
Chris@49:   arma_hot
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   apply
Chris@49:     (
Chris@49:           Mat<out_eT>& C,
Chris@49:     const Mat<in_eT1>& A,
Chris@49:     const Mat<in_eT2>& B,
Chris@49:     const out_eT alpha = out_eT(1),
Chris@49:     const out_eT beta  = out_eT(0)
Chris@49:     )
Chris@49:     {
Chris@49:     arma_extra_debug_sigprint();
Chris@49:     
Chris@49:     const uword A_n_rows = A.n_rows;
Chris@49:     const uword A_n_cols = A.n_cols;
Chris@49:     
Chris@49:     const uword B_n_rows = B.n_rows;
Chris@49:     const uword B_n_cols = B.n_cols;
Chris@49:     
Chris@49:     if( (do_trans_A == false) && (do_trans_B == false) )
Chris@49:       {
Chris@49:       for(uword row_A = 0; row_A < A_n_rows; ++row_A)
Chris@49:         {
Chris@49:         for(uword col_B = 0; col_B < B_n_cols; ++col_B)
Chris@49:           {
Chris@49:           const in_eT2* B_coldata = B.colptr(col_B);
Chris@49:           
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i = 0; i < B_n_rows; ++i)
Chris@49:             {
Chris@49:             const out_eT val1 = upgrade_val<in_eT1,in_eT2>::apply(A.at(row_A,i));
Chris@49:             const out_eT val2 = upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
Chris@49:             acc += val1 * val2;
Chris@49:             //acc += upgrade_val<in_eT1,in_eT2>::apply(A.at(row_A,i)) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
Chris@49:             }
Chris@49:           
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = acc + beta*C.at(row_A,col_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(row_A,col_B) = alpha*acc + beta*C.at(row_A,col_B);
Chris@49:             }
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:     if( (do_trans_A == true) && (do_trans_B == false) )
Chris@49:       {
Chris@49:       for(uword col_A=0; col_A < A_n_cols; ++col_A)
Chris@49:         {
Chris@49:         // col_A is interpreted as row_A when storing the results in matrix C
Chris@49:         
Chris@49:         const in_eT1* A_coldata = A.colptr(col_A);
Chris@49:         
Chris@49:         for(uword col_B=0; col_B < B_n_cols; ++col_B)
Chris@49:           {
Chris@49:           const in_eT2* B_coldata = B.colptr(col_B);
Chris@49:           
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i=0; i < B_n_rows; ++i)
Chris@49:             {
Chris@49:             acc += upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
Chris@49:             }
Chris@49:         
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = acc + beta*C.at(col_A,col_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,col_B) = alpha*acc + beta*C.at(col_A,col_B);
Chris@49:             }
Chris@49:           
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:     if( (do_trans_A == false) && (do_trans_B == true) )
Chris@49:       {
Chris@49:       for(uword row_A = 0; row_A < A_n_rows; ++row_A)
Chris@49:         {
Chris@49:         for(uword row_B = 0; row_B < B_n_rows; ++row_B)
Chris@49:           {
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i = 0; i < B_n_cols; ++i)
Chris@49:             {
Chris@49:             acc += upgrade_val<in_eT1,in_eT2>::apply(A.at(row_A,i)) * upgrade_val<in_eT1,in_eT2>::apply(B.at(row_B,i));
Chris@49:             }
Chris@49:           
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(row_A,row_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(row_A,row_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(row_A,row_B) = acc + beta*C.at(row_A,row_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(row_A,row_B) = alpha*acc + beta*C.at(row_A,row_B);
Chris@49:             }
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:     if( (do_trans_A == true) && (do_trans_B == true) )
Chris@49:       {
Chris@49:       for(uword row_B=0; row_B < B_n_rows; ++row_B)
Chris@49:         {
Chris@49:         
Chris@49:         for(uword col_A=0; col_A < A_n_cols; ++col_A)
Chris@49:           {
Chris@49:           const in_eT1* A_coldata = A.colptr(col_A);
Chris@49:           
Chris@49:           out_eT acc = out_eT(0);
Chris@49:           for(uword i=0; i < A_n_rows; ++i)
Chris@49:             {
Chris@49:             acc += upgrade_val<in_eT1,in_eT2>::apply(B.at(row_B,i)) * upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]);
Chris@49:             }
Chris@49:         
Chris@49:           if( (use_alpha == false) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == false) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = alpha * acc;
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == false) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = acc + beta*C.at(col_A,row_B);
Chris@49:             }
Chris@49:           else
Chris@49:           if( (use_alpha == true) && (use_beta == true) )
Chris@49:             {
Chris@49:             C.at(col_A,row_B) = alpha*acc + beta*C.at(col_A,row_B);
Chris@49:             }
Chris@49:           
Chris@49:           }
Chris@49:         }
Chris@49:       
Chris@49:       }
Chris@49:     }
Chris@49:     
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! \brief
Chris@49: //! Matrix multplication where the matrices have differing element types.
Chris@49: 
Chris@49: template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
Chris@49: class gemm_mixed
Chris@49:   {
Chris@49:   public:
Chris@49:   
Chris@49:   //! immediate multiplication of matrices A and B, storing the result in C
Chris@49:   template<typename out_eT, typename in_eT1, typename in_eT2>
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   apply
Chris@49:     (
Chris@49:           Mat<out_eT>& C,
Chris@49:     const Mat<in_eT1>& A,
Chris@49:     const Mat<in_eT2>& B,
Chris@49:     const out_eT alpha = out_eT(1),
Chris@49:     const out_eT beta  = out_eT(0)
Chris@49:     )
Chris@49:     {
Chris@49:     arma_extra_debug_sigprint();
Chris@49:     
Chris@49:     Mat<in_eT1> tmp_A;
Chris@49:     Mat<in_eT2> tmp_B;
Chris@49:     
Chris@49:     const bool predo_trans_A = ( (do_trans_A == true) && (is_complex<in_eT1>::value == true) );
Chris@49:     const bool predo_trans_B = ( (do_trans_B == true) && (is_complex<in_eT2>::value == true) );
Chris@49:     
Chris@49:     if(do_trans_A)
Chris@49:       {
Chris@49:       op_htrans::apply_noalias(tmp_A, A);
Chris@49:       }
Chris@49:     
Chris@49:     if(do_trans_B)
Chris@49:       {
Chris@49:       op_htrans::apply_noalias(tmp_B, B);
Chris@49:       }
Chris@49:      
Chris@49:     const Mat<in_eT1>& AA = (predo_trans_A == false) ? A : tmp_A;
Chris@49:     const Mat<in_eT2>& BB = (predo_trans_B == false) ? B : tmp_B;
Chris@49:     
Chris@49:     if( (AA.n_elem <= 64u) && (BB.n_elem <= 64u) )
Chris@49:       {
Chris@49:       gemm_mixed_small<((predo_trans_A) ? false : do_trans_A), ((predo_trans_B) ? false : do_trans_B), use_alpha, use_beta>::apply(C, AA, BB, alpha, beta);
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       gemm_mixed_large<((predo_trans_A) ? false : do_trans_A), ((predo_trans_B) ? false : do_trans_B), use_alpha, use_beta>::apply(C, AA, BB, alpha, beta);
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}