multiplication
C++ Function Reference
1 Signature
void BigDataStatMeth::multiplication(BigDataStatMeth::hdf5Dataset *dsA, BigDataStatMeth::hdf5Dataset *dsB, BigDataStatMeth::hdf5Dataset *dsC, bool transpose_A, bool transpose_B, Rcpp::Nullable< bool > bparal, Rcpp::Nullable< int > hdf5_block, Rcpp::Nullable< int > threads)2 Description
Main matrix multiplication function for HDF5 matrices.
3 Parameters
dsA(BigDataStatMeth::hdf5Dataset *): First input matrix datasetdsB(BigDataStatMeth::hdf5Dataset *): Second input matrix datasetdsC(BigDataStatMeth::hdf5Dataset *): Output matrix datasettranspose_A(bool): Whether to transpose matrix Atranspose_B(bool): Whether to transpose matrix Bbparal(Rcpp::Nullable< bool >): Whether to use parallel processinghdf5_block(Rcpp::Nullable< int >): Block size for HDF5 I/O operationsthreads(Rcpp::Nullable< int >): Number of threads for parallel processing
4 Details
Performs matrix multiplication C = A * B where A, B, and C are HDF5 datasets. Supports parallel processing and block-based computation for memory efficiency.
5 Call Graph
6 Source Code
NoteImplementation
File: inst/include/hdf5Algebra/multiplication.hpp • Lines 242-455
inline void multiplication( BigDataStatMeth::hdf5Dataset* dsA, BigDataStatMeth::hdf5Dataset* dsB, BigDataStatMeth::hdf5Dataset* dsC,
bool transpose_A, bool transpose_B, Rcpp::Nullable<bool> bparal, Rcpp::Nullable<int> hdf5_block, Rcpp::Nullable<int> threads = R_NilValue)
{
try {
// int ihdf5_block;
// hsize_t K = dsA->nrows();
// hsize_t N = dsA->ncols();
// hsize_t L = dsB->ncols();
// hsize_t M = dsB->nrows();
//
hsize_t K, N, L, M;
if (transpose_A) {
K = dsA->ncols(); // t(A): filas de A se vuelven columnas
N = dsA->nrows(); // t(A): columnas de A se vuelven filas
} else {
K = dsA->nrows(); // A normal
N = dsA->ncols();
}
if (transpose_B) {
L = dsB->nrows(); // t(B): columnas de B se vuelven filas
M = dsB->ncols(); // t(B): filas de B se vuelven columnas
} else {
L = dsB->ncols(); // B normal
M = dsB->nrows();
}
// if( hdf5_block.isNotNull()) {
// ihdf5_block = Rcpp::as<int>(hdf5_block);
// } else {
// ihdf5_block = MAXBLOCKSIZE/3;
// }
int ihdf5_block_N, ihdf5_block_M, ihdf5_block_K;
if( hdf5_block.isNotNull()) {
ihdf5_block_N = ihdf5_block_M = ihdf5_block_K = Rcpp::as<int>(hdf5_block);
} else {
BlockSizes blocks = calculate_multiplication_blocks(N, M, K);
ihdf5_block_N = ihdf5_block_M = blocks.output_block;
ihdf5_block_K = blocks.inner_block;
}
if( K == L )
{
std::vector<hsize_t> stride = {1, 1},
block = {1, 1},
vsizetoRead, vstart,
vsizetoReadM, vstartM,
vsizetoReadK, vstartK;
dsC->createDataset( N, M, "real");
if( dsC->getDatasetptr() != nullptr) {
getBlockPositionsSizes_hdf5( N, ihdf5_block_N, vstart, vsizetoRead );
getBlockPositionsSizes_hdf5( M, ihdf5_block_M, vstartM, vsizetoReadM );
getBlockPositionsSizes_hdf5( K, ihdf5_block_K, vstartK, vsizetoReadK );
// getBlockPositionsSizes_hdf5( N, ihdf5_block, vstart, vsizetoRead );
// getBlockPositionsSizes_hdf5( M, ihdf5_block, vstartM, vsizetoReadM );
// getBlockPositionsSizes_hdf5( K, ihdf5_block, vstartK, vsizetoReadK );
// int ithreads = get_number_threads(threads, R_NilValue);
// int chunks = vstart.size()/ithreads;
#pragma omp parallel num_threads( get_number_threads(threads, R_NilValue) ) shared(dsA, dsB, dsC, vstart, vsizetoRead) // chunks
{
#pragma omp for schedule(dynamic) nowait
for (hsize_t ii = 0; ii < vstart.size(); ii++)
{
for (hsize_t jj = 0; jj < vstartM.size(); jj++)
{
Eigen::MatrixXd C_accumulator = Eigen::MatrixXd::Zero(vsizetoReadM[jj], vsizetoRead[ii]);
for (hsize_t kk = 0; kk < vstartK.size(); kk++)
{
hsize_t iColsA = vsizetoReadK[kk],
iRowsA = vsizetoRead[ii],
iColsB = vsizetoReadM[jj],
iRowsB = vsizetoReadK[kk];
std::vector<double> vdA( iRowsA * iColsA );
#pragma omp critical(accessFile)
{
dsA->readDatasetBlock( {vstartK[kk], vstart[ii]}, {vsizetoReadK[kk], vsizetoRead[ii]}, stride, block, vdA.data() );
}
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> A (vdA.data(), vsizetoReadK[kk], vsizetoRead[ii] );
std::vector<double> vdB( iRowsB * iColsB );
#pragma omp critical(accessFile)
{
dsB->readDatasetBlock( {vstartM[jj], vstartK[kk]}, {vsizetoReadM[jj], vsizetoReadK[kk]}, stride, block, vdB.data() );
}
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> B (vdB.data(), vsizetoReadM[jj], vsizetoReadK[kk] );
// C_accumulator += B * A;
// Operación según flags de transposición
if (!transpose_A && !transpose_B) {
C_accumulator += B * A; // A * B
} else if (transpose_A && !transpose_B) {
C_accumulator += B * A.transpose(); // t(A) * B
} else if (!transpose_A && transpose_B) {
C_accumulator += B.transpose() * A; // A * t(B)
} else {
C_accumulator += B.transpose() * A.transpose(); // t(A) * t(B)
}
}
std::vector<double> vdC_final(vsizetoReadM[jj] * vsizetoRead[ii]);
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> C_final_map(vdC_final.data(), vsizetoReadM[jj], vsizetoRead[ii]);
C_final_map = C_accumulator;
std::vector<hsize_t> offset = {vstartM[jj], vstart[ii]};
std::vector<hsize_t> count = {vsizetoReadM[jj], vsizetoRead[ii]};
#pragma omp critical(accessFile)
{
dsC->writeDatasetBlock(vdC_final, offset, count, stride, block);
}
}
}
// #pragma omp for schedule(dynamic) nowait
// for (hsize_t ii = 0; ii < vstart.size(); ii++)
// {
//
// for (hsize_t jj = 0; jj < vstartM.size(); jj++)
// {
//
// for (hsize_t kk = 0; kk < vstartK.size(); kk++)
// {
//
// hsize_t iColsA = vsizetoReadK[kk],
// iRowsA = vsizetoRead[ii],
// iColsB = vsizetoReadM[jj],
// iRowsB = vsizetoReadK[kk];
//
// std::vector<double> vdA( iRowsA * iColsA );
// #pragma omp critical(accessFile)
// {
// dsA->readDatasetBlock( {vstartK[kk], vstart[ii]}, {vsizetoReadK[kk], vsizetoRead[ii]}, stride, block, vdA.data() );
// }
//
// Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> A (vdA.data(), vsizetoReadK[kk], vsizetoRead[ii] );
//
// std::vector<double> vdB( iRowsB * iColsB );
// #pragma omp critical(accessFile)
// {
// dsB->readDatasetBlock( {vstartM[jj], vstartK[kk]}, {vsizetoReadM[jj], vsizetoReadK[kk]}, stride, block, vdB.data() );
// }
//
// Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> B (vdB.data(), vsizetoReadM[jj], vsizetoReadK[kk] );
//
// std::vector<double> vdC( vsizetoReadM[jj] * vsizetoRead[ii] );
// #pragma omp critical(accessFile)
// {
// dsC->readDatasetBlock( {vstartM[jj], vstart[ii]}, {vsizetoReadM[jj], vsizetoRead[ii]}, stride, block, vdC.data() );
// }
//
// Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> C (vdC.data(), B.rows(), A.cols() );
//
// C = C + B * A;
//
// std::vector<hsize_t> offset = {vstartM[jj], vstart[ii]};
// std::vector<hsize_t> count = {vsizetoReadM[jj], vsizetoRead[ii] };
//
// #pragma omp critical(accessFile)
// {
// dsC->writeDatasetBlock(vdC, offset, count, stride, block);
// }
// }
// }
// }
}
}
} else {
throw std::range_error("multiplication error: non-conformable arguments");
}
} catch( H5::FileIException& error ) { // catch failure caused by the H5File operations
checkClose_file(dsA, dsB, dsC);
Rcpp::Rcerr<<"\nc++ c++ exception multiplication (File IException)\n";
// return void();
} catch( H5::DataSetIException& error ) { // catch failure caused by the DataSet operations
checkClose_file(dsA, dsB, dsC);
Rcpp::Rcerr<<"\nc++ exception multiplication (DataSet IException)\n";
// return void();
} catch(std::exception &ex) {
checkClose_file(dsA, dsB, dsC);
Rcpp::Rcerr<<"\nc++ exception multiplication " << ex.what();
// return void();
} catch (...) {
checkClose_file(dsA, dsB, dsC);
Rcpp::Rcerr<<"\nC++ exception multiplication (unknown reason)";
// return void();
}
return void();
}7 Usage Example
#include "BigDataStatMeth.hpp"
// Example usage
auto result = multiplication(...);