Inverse_of_Cholesky_decomposition_outofcore_hdf5

C++ Function Reference

1 Signature

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

2 Description

Out-of-core inverse of Cholesky factor using tiled back-substitution.

3 Parameters

  • InOutDataset (BigDataStatMeth::hdf5Dataset *): Dataset containing Cholesky factor L (overwritten with L^-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 Cholesky factor L (overwritten with L^-1) idim0Number of rows idim1Number of columns dElementsBlockBlock size (not used, tiles are fixed) threadsNumber of threads for parallel processing (optional) Inverts lower triangular matrix using fixed tiles:Processes backward from last tile to firstEach tile inversion updates dependent tilesMemory usage: ~800 MB constantInverts lower triangular matrix L in-place using tiles

5 Call Graph

Function dependencies

6 Source Code

File: inst/include/hdf5Algebra/matrixInvCholesky.hpp • Lines 881-944

inline void Inverse_of_Cholesky_decomposition_outofcore_hdf5(BigDataStatMeth::hdf5Dataset* InOutDataset, 
                                                             int idim0, int idim1, long dElementsBlock, 
                                                             Rcpp::Nullable<int> threads = R_NilValue)
{
    
    try {
        int dimensionSize = idim0;
        int tileSize = 10000;
        int numTiles = (dimensionSize + tileSize - 1) / tileSize;
        
        std::vector<hsize_t> stride = {1,1}, block = {1,1};
        
        // Process tiles backward for inverse: last diagonal to first
        for (int k = numTiles - 1; k >= 0; k--) {
            
            hsize_t kStart = k * tileSize;
            hsize_t kSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - kStart));
            
            // 1. Invert diagonal tile L(k,k) in-place
            Eigen::MatrixXd Lkk(kSize, kSize);
            std::vector<double> vLkk(kSize * kSize);
            InOutDataset->readDatasetBlock({kStart, kStart}, {kSize, kSize}, stride, block, vLkk.data());
            Lkk = Eigen::Map<Eigen::MatrixXd>(vLkk.data(), kSize, kSize);
            
            // Invert lower triangular tile
            Lkk = Lkk.triangularView<Eigen::Lower>().solve(Eigen::MatrixXd::Identity(kSize, kSize));
            
            // Write inverted diagonal tile
            InOutDataset->writeDatasetBlock(Rcpp::wrap(Lkk), {kStart, kStart}, {kSize, kSize}, stride, block, false);
            
            // 2. Update column tiles below diagonal: L(i,k) = -L(i,i) * L_old(i,k) * L(k,k)
            #pragma omp parallel for num_threads(get_number_threads(threads, R_NilValue)) schedule(dynamic)
            for (int i = k + 1; i < numTiles; i++) {
                
                hsize_t iStart = i * tileSize;
                hsize_t iSize = std::min(static_cast<hsize_t>(tileSize), static_cast<hsize_t>(dimensionSize - iStart));
                
                // Load L(i,k) - old value
                Eigen::MatrixXd Lik_old(iSize, kSize);
                std::vector<double> vLik(iSize * kSize);
                InOutDataset->readDatasetBlock({iStart, kStart}, {iSize, kSize}, stride, block, vLik.data());
                Lik_old = Eigen::Map<Eigen::MatrixXd>(vLik.data(), iSize, kSize);
                
                // Load L(i,i) - already inverted
                Eigen::MatrixXd Lii(iSize, iSize);
                std::vector<double> vLii(iSize * iSize);
                InOutDataset->readDatasetBlock({iStart, iStart}, {iSize, iSize}, stride, block, vLii.data());
                Lii = Eigen::Map<Eigen::MatrixXd>(vLii.data(), iSize, iSize);
                
                // Update: L(i,k) = -L(i,i) * L_old(i,k) * L(k,k)
                Eigen::MatrixXd Lik_new = -Lii * Lik_old * Lkk;
                
                // Write updated column tile
                InOutDataset->writeDatasetBlock(Rcpp::wrap(Lik_new), {iStart, kStart}, {iSize, kSize}, stride, block, false);
            }
        }
        
    } catch(std::exception& ex) {
        Rcpp::Rcerr << "c++ exception Inverse_of_Cholesky_decomposition_outofcore_hdf5: " << ex.what();
        return void();
    }
    
    return void();
}

7 Usage Example

#include "BigDataStatMeth.hpp"

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