Inverse_of_Cholesky_decomposition_intermediate_hdf5
C++ Function Reference
1 Signature
void BigDataStatMeth::Inverse_of_Cholesky_decomposition_intermediate_hdf5(BigDataStatMeth::hdf5Dataset *InOutDataset, int idim0, int idim1, long dElementsBlock, Rcpp::Nullable< int > threads)2 Description
Computes inverse of Cholesky factor in-place.
3 Parameters
InOutDataset(BigDataStatMeth::hdf5Dataset *): Dataset containing Cholesky factor L (will be overwritten with L^-1)idim0(int): Number of rowsidim1(int): Number of columnsdElementsBlock(long): Block size for processingthreads(Rcpp::Nullable< int >): Number of threads for parallel processing (optional)
4 Details
InOutDatasetDataset containing Cholesky factor L (will be overwritten with L^-1) idim0Number of rows idim1Number of columns dElementsBlockBlock size for processing threadsNumber of threads for parallel processing (optional) Implements block-wise inversion of lower triangular matrix:Processes matrix in blocks for memory efficiencyOverwrites input with its inverseUses parallel processing for independent computations
5 Call Graph
6 Source Code
NoteImplementation
File: inst/include/hdf5Algebra/matrixInvCholesky.hpp • Lines 711-873
inline void Inverse_of_Cholesky_decomposition_intermediate_hdf5( BigDataStatMeth::hdf5Dataset* InOutDataset,
int idim0, int idim1, long dElementsBlock,
Rcpp::Nullable<int> threads = R_NilValue)
{
try{
int dimensionSize = idim0,
readedCols = 0,
colstoRead,
minimumBlockSize;
std::vector<hsize_t> offset = {0,0},
count = {1,1},
stride = {1,1},
block = {1,1};
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky_decomposition_intermediate_hdf5 -> Inici";
// Set minimum elements in block (mandatory : minimum = 2 * longest line)
if( dElementsBlock < dimensionSize * 2 ) {
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> if(dElementsBlock < dimensionSize * 2)";
minimumBlockSize = dimensionSize * 2;
} else {
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> else";
minimumBlockSize = dElementsBlock;
}
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> getDiagonalfromMatrix(InOutDataset)";
Eigen::VectorXd Diagonal = Rcpp::as<Eigen::VectorXd>(getDiagonalfromMatrix(InOutDataset));
// Set the new diagonal in the result matrix
setDiagonalMatrix( InOutDataset, Rcpp::wrap(Diagonal.cwiseInverse()) );
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> Despres setDiagonal";
if( idim0 == idim1) {
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> idim0 == idim1";
while ( readedCols < dimensionSize ) {
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> readedCols < dimensionSize";
colstoRead = ( -2 * readedCols - 1 + std::sqrt( pow(2*readedCols, 2) - 4 * readedCols + 8 * minimumBlockSize + 1) ) / 2;
if (colstoRead <= 0) {
colstoRead = minimumBlockSize; // Minimum progress to avoid infinite loop
}
if( readedCols + colstoRead > idim0) { // Max size bigger than data to read ?
colstoRead = idim0 - readedCols;
}
size_t potential_elements = static_cast<size_t>(dimensionSize - readedCols) * static_cast<size_t>(readedCols + colstoRead);
size_t optimalSize = getOptimalBlockElements();
if (potential_elements > optimalSize) {
int max_cols = static_cast<int>(MAXCHOLBLOCKSIZE / (dimensionSize - readedCols));
colstoRead = max_cols - readedCols;
if (colstoRead < 1) colstoRead = 1;
}
offset[0] = readedCols;
count[0] = dimensionSize - offset[0];
count[1] = readedCols + colstoRead;
Eigen::MatrixXd verticalData;
std::vector<double> vverticalData( count[0] * count[1] );
InOutDataset->readDatasetBlock( {offset[0], offset[1]}, {count[0], count[1]}, stride, block, vverticalData.data() );
verticalData = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> (vverticalData.data(), count[0], count[1] );
//. 20251121 .// for (int j = 1; j < dimensionSize - offset[0]; j++)
//!!!!!!!!!!!! for (int j = 1; j < dimensionSize - offset[0] && j < static_cast<int>(count[0]); j++)
for (int j = 1; j < dimensionSize - static_cast<int>(offset[0]); j++)
{
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> j < dimensionSize - static_cast<int>(offset[0])";
Eigen::VectorXd vR = Eigen::VectorXd::Zero(j+static_cast<int>(offset[0]));
Eigen::ArrayXd ar_j;
int size_j;
if( j < colstoRead) {
size_j = j;
} else {
size_j = colstoRead;
}
ar_j = verticalData.block( j, static_cast<int>(offset[0]), 1, size_j).transpose().array();
vR = Eigen::VectorXd::Zero(static_cast<int>(offset[0]) + ar_j.size());
#pragma omp parallel for num_threads( get_number_threads(threads, R_NilValue) ) shared (ar_j, j, verticalData, offset, colstoRead, vR) schedule(dynamic)
for (int i = 0; i < static_cast<int>(offset[0]) + ar_j.size() ; i++) {
Eigen::ArrayXd ar_i = verticalData.block( 0, i, size_j, 1).array();
if( static_cast<int>(offset[0]) > 0 ) {
if( j <= colstoRead ) {
if( i < static_cast<int>(offset[0]) ){
vR(i) = (( verticalData.coeff(j, i) + ((ar_j.transpose() * ar_i.transpose()).sum())) * (-1)) / Diagonal[j+static_cast<int>(offset[0])];
} else {
vR(i) = ((ar_j.transpose() * ar_i.transpose()) * (-1)).sum() / Diagonal[j+static_cast<int>(offset[0])];
}
} else {
if( i < static_cast<int>(offset[0]) ){
vR(i) = ( verticalData.coeff(j, i) + ((ar_j.transpose() * ar_i.transpose()).sum()));
} else {
vR(i) = (ar_j.transpose() * ar_i.transpose()).sum();
}
}
} else {
if( j <= colstoRead ) {
vR(i) = ((ar_j.transpose() * ar_i.transpose()) * (-1)).sum() / Diagonal[j+static_cast<int>(offset[0])];
} else {
vR(i) = (ar_j.transpose() * ar_i.transpose()).sum();
}
}
}
verticalData.block(j, 0, 1, vR.size()) = vR.transpose();
}
InOutDataset->writeDatasetBlock( Rcpp::wrap(verticalData), offset, count, stride, block, false);
readedCols = readedCols + colstoRead; // Ho preparem perquè desprès necessitarem llegir a partir de la línea anterior
}
} else {
throw std::range_error("non-conformable arguments");
}
// Rcpp::Rcout<<"\nDins Inverse_of_Cholesky -> Bye Bye...";
} catch( H5::FileIException& error ) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr<<"c++ exception Inverse_of_Cholesky_decomposition_intermediate_hdf5 (File IException)";
return void();
} catch( H5::GroupIException & error ) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr << "c++ exception Inverse_of_Cholesky_decomposition_intermediate_hdf5 (Group IException)";
return void();
} catch( H5::DataSetIException& error ) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr << "c++ exception Inverse_of_Cholesky_decomposition_intermediate_hdf5 (DataSet IException)";
return void();
} catch(std::exception& ex) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr << "c++ exception Inverse_of_Cholesky_decomposition_intermediate_hdf5" << ex.what();
return void();
} catch (...) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr<<"\nC++ exception Inverse_of_Cholesky_decomposition_intermediate_hdf5 (unknown reason)";
return void();
}
return void();
}7 Usage Example
#include "BigDataStatMeth.hpp"
// Example usage
auto result = Inverse_of_Cholesky_decomposition_intermediate_hdf5(...);