Rcpp_block_matrix_mul

C++ Function Reference

1 Signature

Eigen::MatrixXd BigDataStatMeth::Rcpp_block_matrix_mul(T X, U Y, Rcpp::Nullable< int > iblock_size)

2 Description

Block-based matrix multiplication.

3 Parameters

  • X (T): First input matrix
  • Y (U): Second input matrix
  • iblock_size (Rcpp::Nullable< int >): Block size for computation

4 Returns

Result of matrix multiplication

5 Details

Implements efficient block-based matrix multiplication using cache-friendly algorithms.

6 Call Graph

Function dependencies

7 Source Code

File: inst/include/memAlgebra/memMultiplication.hpp • Lines 109-188

inline Eigen::MatrixXd Rcpp_block_matrix_mul( T X, U Y, Rcpp::Nullable<int>  iblock_size)
   {
       
       Eigen::MatrixXd C;
       
       try{

           static_assert(std::is_same<T, Eigen::MatrixXd >::value || 
                         std::is_same<T, Eigen::Map< Eigen::MatrixXd >>::value || 
                         std::is_same<T, Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> >::value || 
                         std::is_same<T, Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>> >::value,
                         "Error - type not allowed");
           
           static_assert(std::is_same<U, Eigen::MatrixXd >::value || 
                         std::is_same<U, Eigen::Map< Eigen::MatrixXd >>::value || 
                         std::is_same<U, Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> >::value || 
                         std::is_same<U, Eigen::Map< Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>> >::value,
                         "Error - type not allowed");
           
           Eigen::MatrixXd A = X,
                           B = Y;
           
           // if(transpX == true){ A = X.transpose(); }
           // if(transpY == true){ B = Y.transpose(); }
           
           // int chunks;//, tid;
           hsize_t block_size;
           
           std::vector<hsize_t> vsizetoReadN, vstartN,
                                vsizetoReadM, vstartM,
                                vsizetoReadK, vstartK;
           
           // unsigned int ithreads;
           hsize_t M = A.rows();
           hsize_t K = A.cols();
           hsize_t N = B.cols();
           
           //. 20251121 .// if( K == B.rows()) {
           if( static_cast<Eigen::Index>(K) == B.rows()) {
               if( iblock_size.isNotNull()) {
                   block_size =  Rcpp::as<int>(iblock_size);  
               } else {
                   block_size =  MAXBLOCKSIZE/3;  
               }
               
               C = Eigen::MatrixXd::Zero(M,N) ;
               if(block_size > std::min( N, std::min(M,K)) )
                   block_size = std::min( N, std::min(M,K)); 
               
               getBlockPositionsSizes( N, block_size, vstartN, vsizetoReadN );
               getBlockPositionsSizes( M, block_size, vstartM, vsizetoReadM );
               getBlockPositionsSizes( K, block_size, vstartK, vsizetoReadK );
               
               for (size_t ii = 0; ii < vstartM.size(); ii++)
               {
                   for (size_t jj = 0; jj < vstartN.size(); jj++)
                   {
                       for(size_t kk = 0; kk < vstartK.size(); kk++)
                       {
                           C.block(vstartM[ii], vstartN[jj], vsizetoReadM[ii], vsizetoReadN[jj]) =  C.block(vstartM[ii], vstartN[jj], vsizetoReadM[ii], vsizetoReadN[jj]) + 
                               (A.block(vstartM[ii], vstartK[kk], vsizetoReadM[ii], vsizetoReadK[kk]) * B.block(vstartK[kk], vstartN[jj], vsizetoReadK[kk], vsizetoReadN[jj]));
                       }
                   }
               }
               
           } else {
               throw std::range_error("non-conformable arguments");
               }
           
           
       } catch(std::exception &ex) {
           Rcpp::Rcout<<"c++ error : Bblock_matrix_mul : " <<ex.what();
           return(C);
           
       } catch(...) { 
           Rf_error("c++ exception in Bblock_matrix_mul (unknown reason)"); 
       }
       
       return(C);
   }

8 Usage Example

#include "BigDataStatMeth.hpp"

// Example usage
auto result = Rcpp_block_matrix_mul(...);