Rcpp_block_matrix_substract_hdf5

C++ Function Reference

1 Signature

BigDataStatMeth::hdf5Dataset * BigDataStatMeth::Rcpp_block_matrix_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 subtraction for HDF5 matrices.

3 Parameters

  • dsA (BigDataStatMeth::hdf5Dataset *): First input matrix dataset
  • dsB (BigDataStatMeth::hdf5Dataset *): Second 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 subtraction C = A - B where A, B, and C are HDF5 datasets. Optimized for large matrices with parallel processing.

6 Call Graph

Function dependencies

7 Source Code

File: inst/include/hdf5Algebra/matrixSubstract.hpp • Lines 69-188

inline BigDataStatMeth::hdf5Dataset*  Rcpp_block_matrix_substract_hdf5( 
            BigDataStatMeth::hdf5Dataset* dsA, BigDataStatMeth::hdf5Dataset* dsB, BigDataStatMeth::hdf5Dataset* dsC,
            hsize_t hdf5_block, bool bparal, Rcpp::Nullable<int> threads  = R_NilValue)
    {
        
        try {
            
            hsize_t K = dsA->nrows();
            hsize_t N = dsA->ncols();
            
            if( K == dsB->nrows() && N == dsB->ncols())
            {
                
                // Parallellization and Block variables 
                // unsigned int ithreads;
                std::vector<hsize_t> stride = {1, 1},
                                     block = {1, 1},
                                     vstart, vsizetoRead;
                // hsize_t isize = hdf5_block + 1;
                
                dsC->createDataset( N, K, "real"); 
                
                // ithreads = get_threads(bparal, threads);
                
                if( K<=N ) {
                    
                    getBlockPositionsSizes( N, 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)
                        for (hsize_t ii = 0; ii < vstart.size(); ii ++)
                        {
                            
                            std::vector<double> vdA( K * vsizetoRead[ii] ); 
                            #pragma omp critical (accessFile)
                            {
                            dsA->readDatasetBlock( {0, vstart[ii]}, { K, vsizetoRead[ii]}, stride, block, vdA.data() );
                            }
                            
                            std::vector<double> vdB( K * vsizetoRead[ii] ); 
                            #pragma omp critical (accessFile)
                            {
                            dsB->readDatasetBlock( {0, vstart[ii]}, {K, vsizetoRead[ii]}, stride, block, vdB.data() );
                            }
                            std::transform (vdA.begin(), vdA.end(),
                                            vdB.begin(), vdA.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(vdA, offset, count, stride, block);
                            }
                        }
                    }
    
                } else {
                    
                    getBlockPositionsSizes( K, hdf5_block, vstart, vsizetoRead );
                    // int chunks = vstart.size()/ithreads;
                    
                    #pragma omp parallel num_threads( get_threads(bparal, threads) ) shared(dsA, dsB, dsC)
                    {
                        #pragma omp for schedule (dynamic)
                        for (hsize_t ii = 0; ii < vstart.size(); ii++)
                        {
                            std::vector<double> vdA( vsizetoRead[ii] * N ); 
                            #pragma omp critical (accessFile)
                            {
                                dsA->readDatasetBlock( {vstart[ii], 0}, { vsizetoRead[ii], N}, stride, block, vdA.data() );
                            }
                            
                            std::vector<double> vdB( vsizetoRead[ii] * N); 
                            #pragma omp critical (accessFile)
                            {
                                dsB->readDatasetBlock( {vstart[ii], 0}, {vsizetoRead[ii], N}, stride, block, vdB.data() );
                            }
                            
                            std::transform (vdA.begin(), vdA.end(),
                                            vdB.begin(), vdA.begin(), std::minus<double>());
                            
                            std::vector<hsize_t> offset = { vstart[ii], 0 };
                            std::vector<hsize_t> count = { vsizetoRead[ii], N };
                            #pragma omp critical (accessFile)
                            {
                                dsC->writeDatasetBlock(vdA, offset, count, stride, block);
                            }
                        }
                    }
                }
            } else {
                Rcpp::Rcout<<"matrix substract error: non-conformable arguments\n";
            }
            
        } catch( H5::FileIException& error ) { 
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_substract_hdf5 (File IException)";
            return(dsC);
        } catch( H5::GroupIException & error ) { 
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_substract_hdf5 (Group IException)";
            return(dsC);
        } catch( H5::DataSetIException& error ) { 
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_substract_hdf5 (DataSet IException)";
            return(dsC);
        } catch(std::exception& ex) {
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nc++ exception Rcpp_block_matrix_substract_hdf5" << ex.what();
            return(dsC);
        } catch (...) {
            checkClose_file(dsA, dsB, dsC);
            Rcpp::Rcerr<<"\nC++ exception Rcpp_block_matrix_substract_hdf5 (unknown reason)";
            return(dsC);
        }
        
        return(dsC);
    }

8 Usage Example

#include "BigDataStatMeth.hpp"

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