RcppbdEigen_hdf5

C++ Function Reference

1 Signature

void BigDataStatMeth::RcppbdEigen_hdf5(std::string filename, std::string strsubgroup, std::string strdataset, int k=0, const std::string &which="LM", int ncv=0, bool bcenter=false, bool bscale=false, double tolerance=1e-10, int max_iter=1000, bool compute_vectors=true, bool bforce=false, Rcpp::Nullable< int > threads=R_NilValue)

2 Description

Main eigenvalue decomposition function for HDF5 matrices using Spectra.

3 Parameters

  • filename (std::string): Path to HDF5 file containing the input matrix
  • strsubgroup (std::string): Group path within the HDF5 file
  • strdataset (std::string): Dataset name containing the matrix
  • k (int): Number of eigenvalues to compute (0 = auto-select based on matrix size)
  • which (const std::string &): Which eigenvalues to compute (LM, SM, LR, SR, LI, SI, LA, SA)
  • ncv (int): Number of Arnoldi vectors (0 = auto-select)
  • bcenter (bool): Whether to center the data before decomposition
  • bscale (bool): Whether to scale the data before decomposition
  • tolerance (double): Convergence tolerance for Spectra algorithms
  • max_iter (int): Maximum iterations for Spectra algorithms
  • compute_vectors (bool): Whether to compute eigenvectors or just eigenvalues
  • bforce (bool): Whether to overwrite existing results
  • threads (Rcpp::Nullable< int >): Number of threads for parallel computation

4 Details

This is the main interface for eigenvalue decomposition of matrices stored in HDF5. It uses Spectra library for consistent results with RSpectra and automatically handles both symmetric and non-symmetric matrices.

5 Call Graph

Function dependencies

6 Source Code

File: inst/include/hdf5Algebra/matrixEigenDecomposition.hpp • Lines 456-568

inline void RcppbdEigen_hdf5(std::string filename, std::string strsubgroup, std::string strdataset,
                                 int k = 0, const std::string& which = "LM", int ncv = 0,
                                 bool bcenter = false, bool bscale = false,
                                 double tolerance = 1e-10, int max_iter = 1000, 
                                 bool compute_vectors = true, bool bforce = false,
                                 Rcpp::Nullable<int> threads = R_NilValue) {
        
        BigDataStatMeth::hdf5Dataset* dsA = nullptr;
        BigDataStatMeth::hdf5Dataset* dsd = nullptr;
        BigDataStatMeth::hdf5Dataset* dsu = nullptr;
        
        try {
            
            std::vector<hsize_t> stride = {1, 1}, block = {1, 1}, offset = {0, 0}, count = {0, 0};
            
            dsA = new BigDataStatMeth::hdf5Dataset(filename, strsubgroup, strdataset, false);
            dsA->openDataset();
            
            if (dsA->getDatasetptr() != nullptr) {
                
                // Create results folder following SVD pattern
                std::string stroutgroup = "EIGEN/" + strdataset;
                
                std::vector<hsize_t> dims_out = {dsA->nrows(), dsA->ncols()};
                count = {dims_out[0], dims_out[1]};
                
                // Check if matrix is square
                if (dims_out[0] != dims_out[1]) {
                    Rf_error("Matrix must be square for eigendecomposition");
                    return;
                }
                
                hsize_t n = dims_out[0];
                
                // Use parameter validation function
                std::tie(k, ncv) = validateSpectraParams((int)n, k, ncv);
                
                // Create output datasets
                dsd = new BigDataStatMeth::hdf5Dataset(filename, stroutgroup, "values", bforce);
                if (compute_vectors) {
                    dsu = new BigDataStatMeth::hdf5Dataset(filename, stroutgroup, "vectors", bforce);
                }
                
                // Small matrices => Direct eigendecomposition with Spectra
                if (n * n < MAXELEMSINBLOCK / 20) {
                    
                    Eigen::MatrixXd X;
                    eigdecomp reteig;
                    
                    std::vector<double> vdA(count[0] * count[1]);
                    dsA->readDatasetBlock({offset[0], offset[1]}, {count[0], count[1]}, stride, block, vdA.data());
                    
                    X = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
                        vdA.data(), count[0], count[1]);
                    
                    reteig = RcppbdEigen_spectra(X, k, which, ncv, bcenter, bscale, tolerance, max_iter);
                    
                    if (!reteig.bconv) {
                        Rf_error("Eigendecomposition failed to converge");
                        return;
                    }
                    
                    // Write real eigenvalues
                    dsd->createDataset(1, reteig.eigenvalues_real.size(), "real");
                    dsd->writeDataset(Rcpp::wrap(reteig.eigenvalues_real));
                    
                    // Write real eigenvectors if computed
                    if (compute_vectors && reteig.bcomputevectors) {
                        dsu->createDataset(reteig.eigenvectors_real.rows(), reteig.eigenvectors_real.cols(), "real");
                        dsu->writeDataset(Rcpp::wrap(reteig.eigenvectors_real));
                    }
                    
                    // Write imaginary parts if non-zero (for non-symmetric matrices)
                    if (!reteig.is_symmetric && reteig.eigenvalues_imag.cwiseAbs().maxCoeff() > 1e-14) {
                        
                        BigDataStatMeth::hdf5Dataset* dsd_imag = new BigDataStatMeth::hdf5Dataset(filename, stroutgroup, "values_imag", bforce);
                        dsd_imag->createDataset(1, reteig.eigenvalues_imag.size(), "real");
                        dsd_imag->writeDataset(Rcpp::wrap(reteig.eigenvalues_imag));
                        delete dsd_imag;
                        
                        if (compute_vectors && reteig.bcomputevectors && reteig.eigenvectors_imag.cwiseAbs().maxCoeff() > 1e-14) {
                            BigDataStatMeth::hdf5Dataset* dsu_imag = new BigDataStatMeth::hdf5Dataset(filename, stroutgroup, "vectors_imag", bforce);
                            dsu_imag->createDataset(reteig.eigenvectors_imag.rows(), reteig.eigenvectors_imag.cols(), "real");
                            dsu_imag->writeDataset(Rcpp::wrap(reteig.eigenvectors_imag));
                            delete dsu_imag;
                        }
                    }
                    
                } else {
                    // Large matrices => Block-based approach with Spectra  
                    if (dsA->getDatasetptr() != nullptr) {
                        RcppbdEigen_hdf5_Block(dsA, dsd, dsu, k, which, ncv, bcenter, bscale, tolerance, max_iter, compute_vectors, threads);
                    }
                }
            }
            delete dsd; dsd = nullptr;
            if (dsu) { delete dsu; dsu = nullptr; }
            delete dsA; dsA = nullptr;
            
        } catch(H5::FileIException& error) {
            checkClose_file(dsA, dsd, dsu);
            Rcpp::Rcerr << "\nc++ exception RcppbdEigen_hdf5 (File IException)\n";
        } catch(H5::DataSetIException& error) {
            checkClose_file(dsA, dsd, dsu);
            Rcpp::Rcerr << "\nc++ exception RcppbdEigen_hdf5 (DataSet IException) "<<error.getDetailMsg() <<" \n";
        } catch(std::exception &ex) {
            checkClose_file(dsA, dsd, dsu);
            Rcpp::Rcerr << "c++ exception RcppbdEigen_hdf5: " << ex.what();
        } catch (...) {
            checkClose_file(dsA, dsd, dsu);
            Rcpp::Rcerr << "\nC++ exception RcppbdEigen_hdf5 (unknown reason)";
        }
    }

7 Usage Example

#include "BigDataStatMeth.hpp"

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