Inverse_Matrix_Cholesky_intermediate_hdf5
C++ Function Reference
1 Signature
void BigDataStatMeth::Inverse_Matrix_Cholesky_intermediate_hdf5(BigDataStatMeth::hdf5Dataset *InOutDataset, int idim0, int idim1, long dElementsBlock, Rcpp::Nullable< int > threads)2 Description
Computes final matrix inverse using inverted Cholesky factors.
3 Parameters
InOutDataset(BigDataStatMeth::hdf5Dataset *): Dataset containing L^-1 (will be overwritten with A^-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 L^-1 (will be overwritten with A^-1) idim0Number of rows idim1Number of columns dElementsBlockBlock size for processing threadsNumber of threads for parallel processing (optional) Computes A^-1 = L^-T L^-1 where L^-1 is the inverse of Cholesky factor:Uses block matrix multiplication for memory efficiencyImplements parallel processing for block operationsResult is symmetric but only lower triangular part is computed
5 Call Graph
6 Source Code
NoteImplementation
File: inst/include/hdf5Algebra/matrixInvCholesky.hpp • Lines 962-1128
inline void Inverse_Matrix_Cholesky_intermediate_hdf5( BigDataStatMeth::hdf5Dataset* InOutDataset,
int idim0, int idim1, long dElementsBlock,
Rcpp::Nullable<int> threads = R_NilValue)
{
try {
int dimensionSize = idim0,
readedCols = 0,
colstoRead;
//.. 20251121 ..// minimumBlockSize;
Eigen::VectorXd newDiag(idim0);
std::vector<hsize_t> offset = {0,0},
count = {1,1},
stride = {1,1},
block = {1,1};
// Rcpp::Rcout<<"\nDins Inverse_Matrix_Cholesky_intermediate_hdf5 -> Hi...";
// Set minimum elements in block (mandatory : minimum = 2 * longest line)
//.. 20251121 ..// if( dElementsBlock < dimensionSize * 2 ) {
//.. 20251121 ..// minimumBlockSize = dimensionSize * 2;
//.. 20251121 ..// } else {
//.. 20251121 ..// minimumBlockSize = dElementsBlock;
//.. 20251121 ..// }
if( idim0 == idim1 )
{
// Rcpp::Rcout<<"\nDins Inverse_Matrix_Cholesky_parallel -> idim0 == idim1";
while ( readedCols < dimensionSize ) {
// Rcpp::Rcout<<"\nDins Inverse_Matrix_Cholesky_parallel -> readedCols < dimensionSize";
// Set minimum elements in block - prevent overflow and optimize for any matrix size
//.. 20251121 ..// minimumBlockSize = std::min(static_cast<long>(dElementsBlock), static_cast<long>(std::max(2000L, dimensionSize * 2L)));
// Adaptive block size calculation for any matrix dimensions
// Target: ~50MB blocks for optimal cache usage and memory efficiency
// const int targetElements = 6250000; // ~50MB / 8 bytes per double
const int targetElements = (dimensionSize > 10000) ? 25000000 : 6250000;
const int minCols = std::max(1, std::min(16, dimensionSize / 100)); // 1-16 or 1% of matrix
const int maxCols = std::min(dimensionSize, static_cast<int>(std::sqrt(targetElements / dimensionSize))); // Cache-friendly limit
colstoRead = std::max(minCols, std::min(maxCols, (dimensionSize - readedCols + 2) / 3));
// Ensure we don't exceed remaining columns
colstoRead = std::min(colstoRead, dimensionSize - readedCols);
// Final safety check
if (colstoRead <= 0) colstoRead = 1;
// colstoRead = ( -2 * readedCols - 1 + std::sqrt( pow(2*readedCols, 2) - 4 * readedCols + 8 * minimumBlockSize + 1) ) / 2;
if(colstoRead ==1) {
colstoRead = 2;
}
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 < 2) colstoRead = 2; // Esta función necesita mínimo 2
}
offset[0] = readedCols;
count[0] = dimensionSize - readedCols;
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] );
// Rcpp::Rcout<<"\nDins Inverse_Matrix_Cholesky_parallel -> Anem a per el for parallel";
//!!!!! for ( int i = 0; i < colstoRead + offset[0]; i++) // Columns
// #pragma omp parallel for num_threads(ithreads) shared (verticalData, colstoRead, offset) schedule(dynamic)
// int max_i = std::min(static_cast<int>(colstoRead + offset[0]), static_cast<int>(count[1]));
if (offset[0] > static_cast<hsize_t>(std::numeric_limits<Eigen::Index>::max())) {
Rcpp::stop("offset[0] is too large to convert to Eigen::Index");
}
#pragma omp parallel for num_threads( get_number_threads(threads, R_NilValue) ) shared (verticalData, colstoRead, offset) schedule(static) ordered
for ( int i = 0; i < colstoRead + static_cast<int>(offset[0]); i++) // Columns
// for ( int i = 0; i < max_i; i++) // Columns
{
int init;
if(static_cast<int>(offset[0]) == 0) {
init = i + 1;
newDiag(i) = verticalData.block(i, i, idim0-i, 1 ).array().pow(2).sum();
} else {
if( i < static_cast<int>(offset[0])) {
init = 0;
} else {
newDiag(i) = verticalData.block( i-static_cast<int>(offset[0]), i, idim0-i, 1 ).array().pow(2).sum();
if( i < static_cast<int>(offset[0]) + colstoRead - 1) {
init = i - static_cast<int>(offset[0]) + 1;
} else {
init = colstoRead; // force end
}
}
}
#pragma omp ordered
{
for ( int j = init; j < colstoRead ; j++) { // Rows
if( static_cast<int>(offset[0]) + j < verticalData.cols()) {
verticalData(j,i) = (verticalData.block( j, i , verticalData.rows() - j, 1).array() * verticalData.block( j, j + static_cast<int>(offset[0]), verticalData.rows() - j, 1).array()).sum();
}
}
}
}
// Rcpp::Rcout<<"\nDins Inverse_Matrix_Cholesky_parallel -> Fi volta, anem a escriure";
InOutDataset->writeDatasetBlock( Rcpp::wrap(verticalData), offset, count, stride, block, false);
readedCols = readedCols + colstoRead;
}
setDiagonalMatrix( InOutDataset, Rcpp::wrap(newDiag) );
} else {
throw std::range_error("non-conformable arguments");
}
// Rcpp::Rcout<<"\nDins Inverse_Matrix_Cholesky_parallel -> Bye Bye...";
} catch( H5::FileIException& error ) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr<<"c++ exception Inverse_Matrix_Cholesky_intermediate_hdf5 (File IException)";
return void();
} catch( H5::GroupIException & error ) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr << "c++ exception Inverse_Matrix_Cholesky_intermediate_hdf5 (Group IException)";
return void();
} catch( H5::DataSetIException& error ) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr << "c++ exception Inverse_Matrix_Cholesky_intermediate_hdf5 (DataSet IException)";
return void();
} catch(std::exception& ex) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr << "c++ exception Inverse_Matrix_Cholesky_intermediate_hdf5: " << ex.what();
return void();
} catch (...) {
checkClose_file(InOutDataset);
InOutDataset = nullptr;
Rcpp::Rcerr<<"\nC++ exception Inverse_Matrix_Cholesky_intermediate_hdf5 (unknown reason)";
return void();
}
return void();
}7 Usage Example
#include "BigDataStatMeth.hpp"
// Example usage
auto result = Inverse_Matrix_Cholesky_intermediate_hdf5(...);