Inverse_Matrix_Cholesky_outofcore_hdf5

C++ Function Reference

1 Signature

void BigDataStatMeth::Inverse_Matrix_Cholesky_outofcore_hdf5(BigDataStatMeth::hdf5Dataset *InOutDataset, int idim0, int idim1, long dElementsBlock, Rcpp::Nullable< int > threads)

2 Description

Out-of-core computation of A^-1 = L^-T L^-1 using tiles.

3 Parameters

InOutDataset (BigDataStatMeth::hdf5Dataset *): Dataset containing L^-1 (overwritten with A^-1)
idim0 (int): Number of rows
idim1 (int): Number of columns
dElementsBlock (long): Block size (not used, tiles are fixed)
threads (Rcpp::Nullable< int >): Number of threads for parallel processing (optional)

4 Details

InOutDatasetDataset containing L^-1 (overwritten with A^-1) idim0Number of rows idim1Number of columns dElementsBlockBlock size (not used, tiles are fixed) threadsNumber of threads for parallel processing (optional) Computes final inverse from inverted Cholesky factor:Accumulates triple-nested tile productsExploits symmetry (only lower triangle computed)Memory usage: ~1.6 GB constantComputes final inverse from inverted Cholesky factor Formula: A^-1(i,j) = sum_k L^-1(k,i) * L^-1(k,j) Reads column tiles vertically to compute dot products

5 Call Graph

6 Source Code

Implementation

File: inst/include/hdf5Algebra/matrixInvCholesky.hpp • Lines 1137-1237

inline void Inverse_Matrix_Cholesky_outofcore_hdf5(BigDataStatMeth::hdf5Dataset* InOutDataset, 
                                                   int idim0, int idim1, long dElementsBlock, 
                                                   Rcpp::Nullable<int> threads = R_NilValue)
{
    
    //  BigDataStatMeth::hdf5Dataset* dsA = nullptr;
    
    try {
        int dimensionSize = idim0;
        int tileSize = 10000; // Fixed tile size for memory control
        int numTiles = (dimensionSize + tileSize - 1) / tileSize;
        
        std::vector<hsize_t> stride = {1,1}, block = {1,1};
        
        // Create temporary dataset to avoid overwriting L^-1 during computation
        std::string tempDatasetPath = InOutDataset->getGroup() + "/.tmp_inverse_" + InOutDataset->getDatasetName();
        BigDataStatMeth::hdf5DatasetInternal tempDataset(InOutDataset->getFileName(), tempDatasetPath, true);
        tempDataset.createDataset(idim0, idim1, "real");
        
        // Parallel computation over result tiles (i,j)
        // Result is symmetric, only compute lower triangle
#pragma omp parallel for num_threads(get_number_threads(threads, R_NilValue)) schedule(dynamic)
        for (int i = 0; i < numTiles; i++) {
            
            hsize_t iStart = i * tileSize;
            hsize_t iSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - iStart));
            
            // Only compute j <= i (lower triangle)
            for (int j = 0; j <= i; j++) {
                
                hsize_t jStart = j * tileSize;
                hsize_t jSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - jStart));
                
                // Initialize result tile
                Eigen::MatrixXd Rij = Eigen::MatrixXd::Zero(iSize, jSize);
                
                // Sum over k: A^-1(i,j) = sum_k L^-1(k,i) * L^-1(k,j)
                // L^-1 is lower triangular, so L^-1(k,i) = 0 for k < i
                // Start from max(i,j) to skip zero products
                int k_start = (i > j) ? i : j;
                
                for (int k = k_start; k < numTiles; k++) {
                    
                    hsize_t kStart = k * tileSize;
                    hsize_t kSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - kStart));
                    
                    // Read column tile i: rows k, columns i
                    // This gives L^-1(k,i) which we need for the product
                    Eigen::MatrixXd Lki(kSize, iSize);
                    std::vector<double> vLki(kSize * iSize);
                    InOutDataset->readDatasetBlock({kStart, iStart}, {kSize, iSize}, stride, block, vLki.data());
                    Lki = Eigen::Map<Eigen::MatrixXd>(vLki.data(), kSize, iSize);
                    
                    // Read column tile j: rows k, columns j  
                    // This gives L^-1(k,j)
                    Eigen::MatrixXd Lkj(kSize, jSize);
                    std::vector<double> vLkj(kSize * jSize);
                    InOutDataset->readDatasetBlock({kStart, jStart}, {kSize, jSize}, stride, block, vLkj.data());
                    Lkj = Eigen::Map<Eigen::MatrixXd>(vLkj.data(), kSize, jSize);
                    
                    // Accumulate: R(i,j) += L^-1(k,i)^T * L^-1(k,j)
                    // Lki has rows k, cols i → transpose gives rows i, cols k
                    // Result: (i × k) × (k × j) = (i × j)
                    Rij += Lki.transpose() * Lkj;
                }
                
                // Write result tile to temporary dataset
                tempDataset.writeDatasetBlock(Rcpp::wrap(Rij), {iStart, jStart}, {iSize, jSize}, stride, block, false);
            }
        }
        
        // Copy result from temporary to original dataset
        for (int i = 0; i < numTiles; i++) {
            hsize_t iStart = i * tileSize;
            hsize_t iSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - iStart));
            
            for (int j = 0; j <= i; j++) {
                hsize_t jStart = j * tileSize;
                hsize_t jSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - jStart));
                
                // Read tile from temporary dataset
                Eigen::MatrixXd Rij(iSize, jSize);
                std::vector<double> vRij(iSize * jSize);
                tempDataset.readDatasetBlock({iStart, jStart}, {iSize, jSize}, stride, block, vRij.data());
                Rij = Eigen::Map<Eigen::MatrixXd>(vRij.data(), iSize, jSize);
                
                // Write to original dataset
                InOutDataset->writeDatasetBlock(Rcpp::wrap(Rij), {iStart, jStart}, {iSize, jSize}, stride, block, false);
            }
        }
        
        // Clean up temporary dataset
        tempDataset.remove();
        
    } catch(std::exception& ex) {
        Rcpp::Rcerr << "c++ exception Inverse_Matrix_Cholesky_outofcore_hdf5: " << ex.what();
        return void();
    }
    
    return void();
}

7 Usage Example

#include "BigDataStatMeth.hpp"

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

--- title: "Inverse_Matrix_Cholesky_outofcore_hdf5" subtitle: "C++ Function Reference" --- ## Signature ```cpp void BigDataStatMeth::Inverse_Matrix_Cholesky_outofcore_hdf5(BigDataStatMeth::hdf5Dataset *InOutDataset, int idim0, int idim1, long dElementsBlock, Rcpp::Nullable< int > threads) ``` ## Description Out-of-core computation of A^-1 = L^-T L^-1 using tiles. ## Parameters - `InOutDataset` (`BigDataStatMeth::hdf5Dataset *`): Dataset containing L^-1 (overwritten with A^-1) - `idim0` (`int`): Number of rows - `idim1` (`int`): Number of columns - `dElementsBlock` (`long`): Block size (not used, tiles are fixed) - `threads` (`Rcpp::Nullable< int >`): Number of threads for parallel processing (optional) ## Details InOutDatasetDataset containing L^-1 (overwritten with A^-1) idim0Number of rows idim1Number of columns dElementsBlockBlock size (not used, tiles are fixed) threadsNumber of threads for parallel processing (optional) Computes final inverse from inverted Cholesky factor:Accumulates triple-nested tile productsExploits symmetry (only lower triangle computed)Memory usage: ~1.6 GB constantComputes final inverse from inverted Cholesky factor Formula: A^-1(i,j) = sum_k L^-1(k,i) * L^-1(k,j) Reads column tiles vertically to compute dot products ## Call Graph ::: {.diagram-section} <object type="image/svg+xml" data="images/namespace_big_data_stat_meth_a6cf11002c5ce11473daebd0aefb52444_cgraph.svg" class="class-diagram"> <img src="images/namespace_big_data_stat_meth_a6cf11002c5ce11473daebd0aefb52444_cgraph.svg" alt="Function dependencies"/> </object> ::: ## Source Code ::: {.callout-note collapse="true"} ### Implementation **File:** `inst/include/hdf5Algebra/matrixInvCholesky.hpp` • **Lines 1137-1237** ```cpp inline void Inverse_Matrix_Cholesky_outofcore_hdf5(BigDataStatMeth::hdf5Dataset* InOutDataset, int idim0, int idim1, long dElementsBlock, Rcpp::Nullable<int> threads = R_NilValue) { // BigDataStatMeth::hdf5Dataset* dsA = nullptr; try { int dimensionSize = idim0; int tileSize = 10000; // Fixed tile size for memory control int numTiles = (dimensionSize + tileSize - 1) / tileSize; std::vector<hsize_t> stride = {1,1}, block = {1,1}; // Create temporary dataset to avoid overwriting L^-1 during computation std::string tempDatasetPath = InOutDataset->getGroup() + "/.tmp_inverse_" + InOutDataset->getDatasetName(); BigDataStatMeth::hdf5DatasetInternal tempDataset(InOutDataset->getFileName(), tempDatasetPath, true); tempDataset.createDataset(idim0, idim1, "real"); // Parallel computation over result tiles (i,j) // Result is symmetric, only compute lower triangle #pragma omp parallel for num_threads(get_number_threads(threads, R_NilValue)) schedule(dynamic) for (int i = 0; i < numTiles; i++) { hsize_t iStart = i * tileSize; hsize_t iSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - iStart)); // Only compute j <= i (lower triangle) for (int j = 0; j <= i; j++) { hsize_t jStart = j * tileSize; hsize_t jSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - jStart)); // Initialize result tile Eigen::MatrixXd Rij = Eigen::MatrixXd::Zero(iSize, jSize); // Sum over k: A^-1(i,j) = sum_k L^-1(k,i) * L^-1(k,j) // L^-1 is lower triangular, so L^-1(k,i) = 0 for k < i // Start from max(i,j) to skip zero products int k_start = (i > j) ? i : j; for (int k = k_start; k < numTiles; k++) { hsize_t kStart = k * tileSize; hsize_t kSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - kStart)); // Read column tile i: rows k, columns i // This gives L^-1(k,i) which we need for the product Eigen::MatrixXd Lki(kSize, iSize); std::vector<double> vLki(kSize * iSize); InOutDataset->readDatasetBlock({kStart, iStart}, {kSize, iSize}, stride, block, vLki.data()); Lki = Eigen::Map<Eigen::MatrixXd>(vLki.data(), kSize, iSize); // Read column tile j: rows k, columns j // This gives L^-1(k,j) Eigen::MatrixXd Lkj(kSize, jSize); std::vector<double> vLkj(kSize * jSize); InOutDataset->readDatasetBlock({kStart, jStart}, {kSize, jSize}, stride, block, vLkj.data()); Lkj = Eigen::Map<Eigen::MatrixXd>(vLkj.data(), kSize, jSize); // Accumulate: R(i,j) += L^-1(k,i)^T * L^-1(k,j) // Lki has rows k, cols i → transpose gives rows i, cols k // Result: (i × k) × (k × j) = (i × j) Rij += Lki.transpose() * Lkj; } // Write result tile to temporary dataset tempDataset.writeDatasetBlock(Rcpp::wrap(Rij), {iStart, jStart}, {iSize, jSize}, stride, block, false); } } // Copy result from temporary to original dataset for (int i = 0; i < numTiles; i++) { hsize_t iStart = i * tileSize; hsize_t iSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - iStart)); for (int j = 0; j <= i; j++) { hsize_t jStart = j * tileSize; hsize_t jSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - jStart)); // Read tile from temporary dataset Eigen::MatrixXd Rij(iSize, jSize); std::vector<double> vRij(iSize * jSize); tempDataset.readDatasetBlock({iStart, jStart}, {iSize, jSize}, stride, block, vRij.data()); Rij = Eigen::Map<Eigen::MatrixXd>(vRij.data(), iSize, jSize); // Write to original dataset InOutDataset->writeDatasetBlock(Rcpp::wrap(Rij), {iStart, jStart}, {iSize, jSize}, stride, block, false); } } // Clean up temporary dataset tempDataset.remove(); } catch(std::exception& ex) { Rcpp::Rcerr << "c++ exception Inverse_Matrix_Cholesky_outofcore_hdf5: " << ex.what(); return void(); } return void(); } ``` ::: ## Usage Example ```cpp #include "BigDataStatMeth.hpp" // Example usage auto result = Inverse_Matrix_Cholesky_outofcore_hdf5(...); ```