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 gemm_mixed
max@0: //! @{
max@0: 
max@0: 
max@0: 
max@0: //! \brief
max@0: //! Matrix multplication where the matrices have differing element types.
max@0: //! Uses caching for speedup.
max@0: //! Matrix 'C' is assumed to have been set to the correct size (i.e. taking into account transposes)
max@0: 
max@0: template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
max@0: class gemm_mixed_large
max@0:   {
max@0:   public:
max@0:   
max@0:   template<typename out_eT, typename in_eT1, typename in_eT2>
max@0:   arma_hot
max@0:   inline
max@0:   static
max@0:   void
max@0:   apply
max@0:     (
max@0:           Mat<out_eT>& C,
max@0:     const Mat<in_eT1>& A,
max@0:     const Mat<in_eT2>& B,
max@0:     const out_eT alpha = out_eT(1),
max@0:     const out_eT beta  = out_eT(0)
max@0:     )
max@0:     {
max@0:     arma_extra_debug_sigprint();
max@0:     
max@0:     const uword A_n_rows = A.n_rows;
max@0:     const uword A_n_cols = A.n_cols;
max@0:     
max@0:     const uword B_n_rows = B.n_rows;
max@0:     const uword B_n_cols = B.n_cols;
max@0:     
max@0:     if( (do_trans_A == false) && (do_trans_B == false) )
max@0:       {
max@0:       podarray<in_eT1> tmp(A_n_cols);
max@0:       in_eT1* A_rowdata = tmp.memptr();
max@0:       
max@0:       for(uword row_A=0; row_A < A_n_rows; ++row_A)
max@0:         {
max@0:         tmp.copy_row(A, row_A);
max@0:         
max@0:         for(uword col_B=0; col_B < B_n_cols; ++col_B)
max@0:           {
max@0:           const in_eT2* B_coldata = B.colptr(col_B);
max@0:           
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i=0; i < B_n_rows; ++i)
max@0:             {
max@0:             acc += upgrade_val<in_eT1,in_eT2>::apply(A_rowdata[i]) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
max@0:             }
max@0:         
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(row_A,col_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(row_A,col_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(row_A,col_B) = acc + beta*C.at(row_A,col_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(row_A,col_B) = alpha*acc + beta*C.at(row_A,col_B);
max@0:             }
max@0:           
max@0:           }
max@0:         }
max@0:       }
max@0:     else
max@0:     if( (do_trans_A == true) && (do_trans_B == false) )
max@0:       {
max@0:       for(uword col_A=0; col_A < A_n_cols; ++col_A)
max@0:         {
max@0:         // col_A is interpreted as row_A when storing the results in matrix C
max@0:         
max@0:         const in_eT1* A_coldata = A.colptr(col_A);
max@0:         
max@0:         for(uword col_B=0; col_B < B_n_cols; ++col_B)
max@0:           {
max@0:           const in_eT2* B_coldata = B.colptr(col_B);
max@0:           
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i=0; i < B_n_rows; ++i)
max@0:             {
max@0:             acc += upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
max@0:             }
max@0:         
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,col_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,col_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,col_B) = acc + beta*C.at(col_A,col_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,col_B) = alpha*acc + beta*C.at(col_A,col_B);
max@0:             }
max@0:           
max@0:           }
max@0:         }
max@0:       }
max@0:     else
max@0:     if( (do_trans_A == false) && (do_trans_B == true) )
max@0:       {
max@0:       Mat<in_eT2> B_tmp;
max@0:       
max@0:       op_strans::apply_noalias(B_tmp, B);
max@0:       
max@0:       gemm_mixed_large<false, false, use_alpha, use_beta>::apply(C, A, B_tmp, alpha, beta);
max@0:       }
max@0:     else
max@0:     if( (do_trans_A == true) && (do_trans_B == true) )
max@0:       {
max@0:       // mat B_tmp = trans(B);
max@0:       // dgemm_arma<true, false,  use_alpha, use_beta>::apply(C, A, B_tmp, alpha, beta);
max@0:       
max@0:       
max@0:       // By using the trans(A)*trans(B) = trans(B*A) equivalency,
max@0:       // transpose operations are not needed
max@0:       
max@0:       podarray<in_eT2> tmp(B_n_cols);
max@0:       in_eT2* B_rowdata = tmp.memptr();
max@0:       
max@0:       for(uword row_B=0; row_B < B_n_rows; ++row_B)
max@0:         {
max@0:         tmp.copy_row(B, row_B);
max@0:         
max@0:         for(uword col_A=0; col_A < A_n_cols; ++col_A)
max@0:           {
max@0:           const in_eT1* A_coldata = A.colptr(col_A);
max@0:           
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i=0; i < A_n_rows; ++i)
max@0:             {
max@0:             acc += upgrade_val<in_eT1,in_eT2>::apply(B_rowdata[i]) * upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]);
max@0:             }
max@0:         
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,row_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,row_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,row_B) = acc + beta*C.at(col_A,row_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,row_B) = alpha*acc + beta*C.at(col_A,row_B);
max@0:             }
max@0:           
max@0:           }
max@0:         }
max@0:       
max@0:       }
max@0:     }
max@0:     
max@0:   };
max@0: 
max@0: 
max@0: 
max@0: //! Matrix multplication where the matrices have different element types.
max@0: //! Simple version (no caching).
max@0: //! Matrix 'C' is assumed to have been set to the correct size (i.e. taking into account transposes)
max@0: template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
max@0: class gemm_mixed_small
max@0:   {
max@0:   public:
max@0:   
max@0:   template<typename out_eT, typename in_eT1, typename in_eT2>
max@0:   arma_hot
max@0:   inline
max@0:   static
max@0:   void
max@0:   apply
max@0:     (
max@0:           Mat<out_eT>& C,
max@0:     const Mat<in_eT1>& A,
max@0:     const Mat<in_eT2>& B,
max@0:     const out_eT alpha = out_eT(1),
max@0:     const out_eT beta  = out_eT(0)
max@0:     )
max@0:     {
max@0:     arma_extra_debug_sigprint();
max@0:     
max@0:     const uword A_n_rows = A.n_rows;
max@0:     const uword A_n_cols = A.n_cols;
max@0:     
max@0:     const uword B_n_rows = B.n_rows;
max@0:     const uword B_n_cols = B.n_cols;
max@0:     
max@0:     if( (do_trans_A == false) && (do_trans_B == false) )
max@0:       {
max@0:       for(uword row_A = 0; row_A < A_n_rows; ++row_A)
max@0:         {
max@0:         for(uword col_B = 0; col_B < B_n_cols; ++col_B)
max@0:           {
max@0:           const in_eT2* B_coldata = B.colptr(col_B);
max@0:           
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i = 0; i < B_n_rows; ++i)
max@0:             {
max@0:             const out_eT val1 = upgrade_val<in_eT1,in_eT2>::apply(A.at(row_A,i));
max@0:             const out_eT val2 = upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
max@0:             acc += val1 * val2;
max@0:             //acc += upgrade_val<in_eT1,in_eT2>::apply(A.at(row_A,i)) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
max@0:             }
max@0:           
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(row_A,col_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(row_A,col_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(row_A,col_B) = acc + beta*C.at(row_A,col_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(row_A,col_B) = alpha*acc + beta*C.at(row_A,col_B);
max@0:             }
max@0:           }
max@0:         }
max@0:       }
max@0:     else
max@0:     if( (do_trans_A == true) && (do_trans_B == false) )
max@0:       {
max@0:       for(uword col_A=0; col_A < A_n_cols; ++col_A)
max@0:         {
max@0:         // col_A is interpreted as row_A when storing the results in matrix C
max@0:         
max@0:         const in_eT1* A_coldata = A.colptr(col_A);
max@0:         
max@0:         for(uword col_B=0; col_B < B_n_cols; ++col_B)
max@0:           {
max@0:           const in_eT2* B_coldata = B.colptr(col_B);
max@0:           
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i=0; i < B_n_rows; ++i)
max@0:             {
max@0:             acc += upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]) * upgrade_val<in_eT1,in_eT2>::apply(B_coldata[i]);
max@0:             }
max@0:         
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,col_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,col_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,col_B) = acc + beta*C.at(col_A,col_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,col_B) = alpha*acc + beta*C.at(col_A,col_B);
max@0:             }
max@0:           
max@0:           }
max@0:         }
max@0:       }
max@0:     else
max@0:     if( (do_trans_A == false) && (do_trans_B == true) )
max@0:       {
max@0:       for(uword row_A = 0; row_A < A_n_rows; ++row_A)
max@0:         {
max@0:         for(uword row_B = 0; row_B < B_n_rows; ++row_B)
max@0:           {
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i = 0; i < B_n_cols; ++i)
max@0:             {
max@0:             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));
max@0:             }
max@0:           
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(row_A,row_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(row_A,row_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(row_A,row_B) = acc + beta*C.at(row_A,row_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(row_A,row_B) = alpha*acc + beta*C.at(row_A,row_B);
max@0:             }
max@0:           }
max@0:         }
max@0:       }
max@0:     else
max@0:     if( (do_trans_A == true) && (do_trans_B == true) )
max@0:       {
max@0:       for(uword row_B=0; row_B < B_n_rows; ++row_B)
max@0:         {
max@0:         
max@0:         for(uword col_A=0; col_A < A_n_cols; ++col_A)
max@0:           {
max@0:           const in_eT1* A_coldata = A.colptr(col_A);
max@0:           
max@0:           out_eT acc = out_eT(0);
max@0:           for(uword i=0; i < A_n_rows; ++i)
max@0:             {
max@0:             acc += upgrade_val<in_eT1,in_eT2>::apply(B.at(row_B,i)) * upgrade_val<in_eT1,in_eT2>::apply(A_coldata[i]);
max@0:             }
max@0:         
max@0:           if( (use_alpha == false) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,row_B) = acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == false) )
max@0:             {
max@0:             C.at(col_A,row_B) = alpha * acc;
max@0:             }
max@0:           else
max@0:           if( (use_alpha == false) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,row_B) = acc + beta*C.at(col_A,row_B);
max@0:             }
max@0:           else
max@0:           if( (use_alpha == true) && (use_beta == true) )
max@0:             {
max@0:             C.at(col_A,row_B) = alpha*acc + beta*C.at(col_A,row_B);
max@0:             }
max@0:           
max@0:           }
max@0:         }
max@0:       
max@0:       }
max@0:     }
max@0:     
max@0:   };
max@0: 
max@0: 
max@0: 
max@0: 
max@0: 
max@0: //! \brief
max@0: //! Matrix multplication where the matrices have differing element types.
max@0: 
max@0: template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
max@0: class gemm_mixed
max@0:   {
max@0:   public:
max@0:   
max@0:   //! immediate multiplication of matrices A and B, storing the result in C
max@0:   template<typename out_eT, typename in_eT1, typename in_eT2>
max@0:   inline
max@0:   static
max@0:   void
max@0:   apply
max@0:     (
max@0:           Mat<out_eT>& C,
max@0:     const Mat<in_eT1>& A,
max@0:     const Mat<in_eT2>& B,
max@0:     const out_eT alpha = out_eT(1),
max@0:     const out_eT beta  = out_eT(0)
max@0:     )
max@0:     {
max@0:     arma_extra_debug_sigprint();
max@0:     
max@0:     Mat<in_eT1> tmp_A;
max@0:     Mat<in_eT2> tmp_B;
max@0:     
max@0:     const bool predo_trans_A = ( (do_trans_A == true) && (is_complex<in_eT1>::value == true) );
max@0:     const bool predo_trans_B = ( (do_trans_B == true) && (is_complex<in_eT2>::value == true) );
max@0:     
max@0:     if(do_trans_A)
max@0:       {
max@0:       op_htrans::apply_noalias(tmp_A, A);
max@0:       }
max@0:     
max@0:     if(do_trans_B)
max@0:       {
max@0:       op_htrans::apply_noalias(tmp_B, B);
max@0:       }
max@0:      
max@0:     const Mat<in_eT1>& AA = (predo_trans_A == false) ? A : tmp_A;
max@0:     const Mat<in_eT2>& BB = (predo_trans_B == false) ? B : tmp_B;
max@0:     
max@0:     if( (AA.n_elem <= 64u) && (BB.n_elem <= 64u) )
max@0:       {
max@0:       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);
max@0:       }
max@0:     else
max@0:       {
max@0:       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);
max@0:       }
max@0:     }
max@0:   
max@0:   
max@0:   };
max@0: 
max@0: 
max@0: 
max@0: //! @}