RcppQRHdf5

C++ Function Reference

1 Signature

void BigDataStatMeth::RcppQRHdf5(BigDataStatMeth::hdf5Dataset *dsA, BigDataStatMeth::hdf5Dataset *dsQ, BigDataStatMeth::hdf5Dataset *dsR, bool bthin, Rcpp::Nullable< int > block_size, Rcpp::Nullable< int > threads)

2 Description

QR decomposition for HDF5 matrices (Householder / LAPACK backend)

3 Parameters

  • dsA (BigDataStatMeth::hdf5Dataset *): Input matrix dataset
  • dsQ (BigDataStatMeth::hdf5Dataset *): Output Q matrix dataset
  • dsR (BigDataStatMeth::hdf5Dataset *): Output R matrix dataset
  • bthin (bool): Whether to compute thin QR
  • block_size (Rcpp::Nullable< int >): Block size for processing (optional, currently unused for LAPACK)
  • threads (Rcpp::Nullable< int >): Number of Eigen internal threads (-1 = auto)

4 Details

Computes the QR decomposition of a matrix stored in HDF5 format. Supports both full and thin decomposition with optional Eigen thread control.

5 Call Graph

Function dependencies

6 Source Code

File: inst/include/hdf5Algebra/matrixQR.hpp • Lines 233-407

inline void RcppQRHdf5( BigDataStatMeth::hdf5Dataset* dsA, 
                        BigDataStatMeth::hdf5Dataset* dsQ, 
                        BigDataStatMeth::hdf5Dataset* dsR, 
                        bool bthin,
                        Rcpp::Nullable<int> block_size = R_NilValue,
                        Rcpp::Nullable<int> threads    = R_NilValue )
{
    
    try {
        // ── Optionally set Eigen BLAS thread count ────────────────────────────
        // Eigen::setNbThreads uses its own internal pool (separate from OpenMP),
        // so it is safe to call here; it respects the getDTthreads limit via the
        // caller capping threads before passing them in.
        if (threads.isNotNull()) {
            int nth = std::max(1, Rcpp::as<int>(threads));
            Eigen::setNbThreads(nth);
        }

        // ── Read A from HDF5 (row-major → transpose to get R column-major) ───
        std::vector<hsize_t> offset = {0,0},
            count  = {dsA->nrows(), dsA->ncols()},
            stride = {1,1},
            block  = {1,1};
        
        std::vector<double> vdA( count[0] * count[1] );
        dsA->readDatasetBlock( {offset[0], offset[1]}, {count[0], count[1]},
                               stride, block, vdA.data() );
        
        // count[0]=ncols_R, count[1]=nrows_R  (HDF5 is transposed vs R)
        Eigen::MatrixXd A = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> (vdA.data(), count[0], count[1]);
        A.transposeInPlace();   // now A is m×n  (R dimensions)
        
        const int m = static_cast<int>(A.rows());
        const int n = static_cast<int>(A.cols());
        const int k = std::min(m, n);   // number of significant components
        
        // ── QR decomposition ─────────────────────────────────────────────────
        Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);
        
        if (!bthin) {
            // ── Full QR: Q is m×m, R is m×n (upper trapezoidal) ─────────────
            // Q
            Eigen::MatrixXd Qfull = qr.householderQ();
            dsQ->inheritCompressionLevel(dsA->getCompressionLevel());
            dsQ->createDataset( Qfull.rows(), Qfull.cols(), "real" );
            dsQ->writeDataset( Rcpp::wrap(Qfull) );
            
            // R: m×n upper trapezoidal (zeros below the first k rows diagonal)
            Eigen::MatrixXd Rfull = Eigen::MatrixXd::Zero(m, n);
            Rfull.topRows(k) = qr.matrixQR().topRows(k)
                 .template triangularView<Eigen::Upper>();
            dsR->inheritCompressionLevel(dsA->getCompressionLevel());
            dsR->createDataset( Rfull.rows(), Rfull.cols(), "real" );
            dsR->writeDataset( Rcpp::wrap(Rfull) );
            
        } else {
            // ── Thin QR: Q is m×k, R is k×n ─────────────────────────────────
            // Q thin
            Eigen::MatrixXd Qthin = qr.householderQ()
            * Eigen::MatrixXd::Identity(m, k);
            dsQ->inheritCompressionLevel(dsA->getCompressionLevel());
            dsQ->createDataset( Qthin.rows(), Qthin.cols(), "real" );
            dsQ->writeDataset( Rcpp::wrap(Qthin) );
            
            // R thin: k×n upper triangular
            Eigen::MatrixXd Rthin = qr.matrixQR().topRows(k)
                                      .template triangularView<Eigen::Upper>();
            dsR->inheritCompressionLevel(dsA->getCompressionLevel());
            dsR->createDataset( Rthin.rows(), Rthin.cols(), "real" );
            dsR->writeDataset( Rcpp::wrap(Rthin) );
        }
        
    } catch( H5::FileIException& error ) {
        throw std::runtime_error("c++ exception RcppQRHdf5 (File IException)");
    } catch( H5::GroupIException & error ) {
        throw std::runtime_error("c++ exception RcppQRHdf5 (Group IException)");
    } catch( H5::DataSetIException& error ) {
        throw std::runtime_error("c++ exception RcppQRHdf5 (DataSet IException)");
    } catch(std::exception& ex) {
        throw std::runtime_error(std::string("c++ exception RcppQRHdf5: ") + ex.what());
    }
    /**
    
    try {

        int irank; //,
            // iblockfactor = 1;
        std::vector<hsize_t> offset = {0,0},
            count = {dsA->nrows(), dsA->ncols()},
            stride = {1,1},
            block = {1,1};
        Eigen::MatrixXd Q;
        
        
        std::vector<double> vdA( count[0] * count[1] ); 
        dsA->readDatasetBlock( {offset[0], offset[1]}, {count[0], count[1]}, stride, block, vdA.data() );
        Eigen::MatrixXd A = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> (vdA.data(), count[0], count[1] );
        A.transposeInPlace();
        
        
        Eigen::FullPivLU<Eigen::MatrixXd>lu_decomp(A);
        Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);
        
        qr.compute(A);
        irank = lu_decomp.rank();
        
        
        //..//if (irank == count[0] + 1 || irank == count[1] + 1 )
        if (static_cast<unsigned long long>(irank) == count[0] + 1ULL ||
            static_cast<unsigned long long>(irank) == count[1] + 1ULL)
        {
            Eigen::MatrixXd R = qr.matrixQR().template triangularView<Eigen::Upper>();
            dsR->inheritCompressionLevel(dsA->getCompressionLevel());
            dsR->createDataset( R.rows(), R.cols(), "real" );
            dsR->writeDataset(Rcpp::wrap(R));
        } else {
            Eigen::MatrixXd R = qr.matrixQR().topLeftCorner(irank, irank).template triangularView<Eigen::Upper>();
            dsR->inheritCompressionLevel(dsA->getCompressionLevel());
            dsR->createDataset( R.rows(), R.cols(), "real" );
            dsR->writeDataset(Rcpp::wrap(R));
        }
        
        
        if (bthin == false)
        {
            Eigen::MatrixXd Q = qr.householderQ();
            dsQ->inheritCompressionLevel(dsA->getCompressionLevel());
            dsQ->createDataset( Q.rows(), Q.cols(), "real" );
            dsQ->writeDataset( Rcpp::wrap(Q) );
        } else {
            if (bthin) {
                Eigen::MatrixXd R = qr.matrixQR()
                                      .topLeftCorner(irank, A.cols())
                                      .triangularView<Eigen::Upper>();
                dsR->inheritCompressionLevel(dsA->getCompressionLevel());
                dsR->createDataset( R.rows(), R.cols(), "real" );
                dsR->writeDataset(Rcpp::wrap(R));
            } else {
                // Full QR: R is m×n upper trapezoidal
                Eigen::MatrixXd R = Eigen::MatrixXd::Zero(A.rows(), A.cols());
                R.topRows(irank) = qr.matrixQR()
                 .topRows(irank)
                 .triangularView<Eigen::Upper>();
                dsR->inheritCompressionLevel(dsA->getCompressionLevel());
                dsR->createDataset( R.rows(), R.cols(), "real" );
                dsR->writeDataset(Rcpp::wrap(R));
            }
        }
        
        //. 2026/03/02 .//  else {
            
            //. 2026/03/02 .//      //. 2025/01/15 error with thin calculus.// Eigen::MatrixXd Qthin = qr.householderQ() * Eigen::MatrixXd::Identity(count[0], count[1]);
            //. 2026/03/02 .//  //. 2025/01/15 .// dsQ->createDataset( Qthin.rows(), Qthin.cols(), "real" );
            //. 2026/03/02 .//  //. 2025/01/15 .// dsQ->writeDataset( Rcpp::wrap(Qthin));
            
            //. 2026/03/02 .//  Eigen::MatrixXd Qthin = qr.householderQ() * Eigen::MatrixXd::Identity(count[1], count[0]);
            //. 2026/03/02 .//  dsQ->createDataset( Qthin.rows(), Qthin.cols(), "real" );
            //. 2026/03/02 .//  dsQ->writeDataset( Rcpp::wrap(Qthin));
            
        //. 2026/03/02 .//  }
        
    } catch( H5::FileIException& error ) { // catch failure caused by the H5File operations
        throw std::runtime_error("c++ exception RcppQRHdf5 (File IException)");
    } catch( H5::GroupIException & error ) { // catch failure caused by the DataSet operations
        throw std::runtime_error("c++ exception RcppQRHdf5 (Group IException)");
    } catch( H5::DataSetIException& error ) { // catch failure caused by the DataSet operations
        throw std::runtime_error("c++ exception RcppQRHdf5 (DataSet IException)");
    } catch(std::exception& ex) {
        throw std::runtime_error(std::string("c++ exception RcppQRHdf5: ") + ex.what());
    }
     **/
    
    return void();
    
}

7 Usage Example

#include "BigDataStatMeth.hpp"

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