Rcpp_block_matrix_vector_substract_hdf5

C++ Function Reference

1 Signature

BigDataStatMeth::hdf5Dataset * BigDataStatMeth::Rcpp_block_matrix_vector_substract_hdf5(BigDataStatMeth::hdf5Dataset *dsA, BigDataStatMeth::hdf5Dataset *dsB, BigDataStatMeth::hdf5Dataset *dsC, hsize_t hdf5_block, bool bparal, Rcpp::Nullable< int > threads=R_NilValue)

2 Description

Block-based matrix-vector subtraction for HDF5 matrices.

3 Parameters

  • dsA (BigDataStatMeth::hdf5Dataset *): Input vector dataset
  • dsB (BigDataStatMeth::hdf5Dataset *): Input matrix dataset
  • dsC (BigDataStatMeth::hdf5Dataset *): Output matrix dataset
  • hdf5_block (hsize_t): Block size for HDF5 I/O operations
  • bparal (bool): Whether to use parallel processing
  • threads (Rcpp::Nullable< int >): Number of threads for parallel processing (optional)

4 Returns

Pointer to result dataset

5 Details

Performs block-based matrix-vector subtraction where one operand is a vector and the other is a matrix. Supports both row and column vectors.

6 Call Graph

Function dependencies

7 Source Code

File: inst/include/hdf5Algebra/matrixSubstract.hpp • Lines 206-350

inline BigDataStatMeth::hdf5Dataset* Rcpp_block_matrix_vector_substract_hdf5( 
            BigDataStatMeth::hdf5Dataset* dsA, BigDataStatMeth::hdf5Dataset* dsB, BigDataStatMeth::hdf5Dataset* dsC,
            hsize_t hdf5_block, bool bparal, Rcpp::Nullable<int> threads  = R_NilValue)
    {
        
        try {
            
            // Vector
            hsize_t K = dsA->nrows();
            hsize_t N = dsA->ncols();
            
            // Matrix
            hsize_t M = dsB->nrows();
            hsize_t L = dsB->ncols();
            
            std::vector<hsize_t> vstart, vsizetoRead;
            // unsigned int ithreads;
            
            if(hdf5_block == 1) {
                hdf5_block = ceil(MAXBLOCKSIZE/(K*N));
            }
            
            // hsize_t isize = hdf5_block + 1;
            std::vector<hsize_t> stride = {1, 1},
                                 block = {1, 1};
            
            dsC->createDataset( L, M, "real");
            
            std::vector<double> vdA( K * N );
            dsA->readDatasetBlock( {0, 0}, {K, N}, stride, block, vdA.data() );
            
            // ithreads = get_threads(bparal, threads);
            
            if(  K == M )
            { // Sum vector to every col
                
                getBlockPositionsSizes( L, hdf5_block, vstart, vsizetoRead );
                // int chunks = vstart.size()/ithreads;
                
                #pragma omp parallel num_threads( get_threads(bparal, threads) ) shared(dsA, dsB, dsC) //, chunks)
                {
                    #pragma omp for schedule (dynamic) // collapse(2)
                    for (hsize_t ii = 0; ii < vstart.size(); ii ++)
                    {
                        std::vector<double> vdB( K * vsizetoRead[ii] );
                        #pragma omp critical (accessFile)
                        {
                            dsB->readDatasetBlock( {0, vstart[ii]}, {K, vsizetoRead[ii]}, stride, block, vdB.data() );
                        }
                        
                        // Duplicate vector
                        std::size_t const no_of_duplicates = (K * vsizetoRead[ii]) / vdA.size();
                        
                        std::vector<double> v = vdA; 
                        v.reserve(vdA.size() * no_of_duplicates);
                        auto end = std::end(v);
                        
                        for(std::size_t i = 1; i < no_of_duplicates; ++i)
                            v.insert(std::end(v), std::begin(v), end);
                        
                        // Sum vector to matrix by columns / rows
                        std::transform (vdB.begin(), vdB.end(),
                                        v.begin(), vdB.begin(), std::minus<double>());
                        
                        
                        std::vector<hsize_t> offset = { 0, vstart[ii]};
                        std::vector<hsize_t> count = { K, vsizetoRead[ii]};
                        
                        #pragma omp critical (accessFile)
                        {
                            dsC->writeDatasetBlock( vdB, offset, count, stride, block);
                        }
                    }
                }
    
            } else if(  K == L ) { // Sum vector to every row
                
                getBlockPositionsSizes( M, hdf5_block, vstart, vsizetoRead );
                // int chunks = vstart.size()/ithreads;
                
                #pragma omp parallel num_threads( get_threads(bparal, threads) ) shared(dsA, dsB, dsC) //, chunks)
                {
                    #pragma omp for schedule (dynamic) // collapse(2)
                    for (hsize_t ii = 0; ii < vstart.size(); ii ++)
                    {
                        std::vector<double> vdB( L * vsizetoRead[ii] );
                        #pragma omp critical (accessFile)
                        {
                            dsB->readDatasetBlock( {vstart[ii], 0}, {vsizetoRead[ii], K}, stride, block, vdB.data() );
                        }
                        Rcpp::NumericMatrix B (vsizetoRead[ii], K, vdB.begin());
                        
                        Rcpp::transpose(B);
                        
                        // Duplicate vector
                        std::size_t const no_of_duplicates = (vsizetoRead[ii] * K) / vdA.size();
                        
                        std::vector<double> v = vdA; 
                        v.reserve(v.size() * no_of_duplicates);
                        auto end = std::end(v);
                        
                        for(std::size_t i = 1; i < no_of_duplicates; ++i)
                            v.insert(std::end(v), std::begin(v), end);
                        
                        // Sum vector to matrix by columns / rows
                        std::transform (B.begin(), B.end(),
                                        v.begin(), B.begin() , std::minus<double>());
                        
                        std::vector<hsize_t> offset = { vstart[ii], 0};
                        std::vector<hsize_t> count = { vsizetoRead[ii], K};
                        
                        #pragma omp critical (accessFile)
                        {
                            dsC->writeDatasetBlock( vdB, offset, count, stride, block);
                        }
                    }
                }
            } else {
                Rcpp::Rcout<< "vector substract error: non-conformable arguments\n";
            }
            
        } catch( H5::FileIException& error ) { 
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_vector_substract_hdf5 (File IException)";
            return(dsC);
        } catch( H5::GroupIException & error ) { 
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_vector_substract_hdf5 (Group IException)";
            return(dsC);
        } catch( H5::DataSetIException& error ) { 
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_vector_substract_hdf5 (DataSet IException)";
            return(dsC);
        } catch(std::exception& ex) {
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_vector_substract_hdf5 " << ex.what();
            return(dsC);
        } catch (...) {
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nC++ exception Rcpp_block_matrix_vector_substract_hdf5 (unknown reason)";
            return(dsC);
        }
        
        return(dsC);
    }

8 Usage Example

#include "BigDataStatMeth.hpp"

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